Mercurial > repos > greg > insect_phenology_model
diff insect_phenology_model.R @ 96:f0c65c9f1078 draft
Uploaded
| author | greg |
|---|---|
| date | Mon, 04 Dec 2017 10:57:59 -0500 |
| parents | 89492de287a6 |
| children | 58bc1a2ca936 |
line wrap: on
line diff
--- a/insect_phenology_model.R Mon Dec 04 10:57:52 2017 -0500 +++ b/insect_phenology_model.R Mon Dec 04 10:57:59 2017 -0500 @@ -111,7 +111,7 @@ dh[i] <- T[i] - 8.4 } } - else if (i > settime) { + else if (i > settime) { n <- i - settime T[i] = dnp + (ts - dnp) * exp( - b * n / z) if (T[i] < 8.4) { @@ -201,9 +201,23 @@ cat("Number of days: ", opt$num_days, "\n") -# Initialize matrix for results from all replications. -S0.rep <- S1.rep <- S2.rep <- S3.rep <- S4.rep <- S5.rep <- matrix(rep(0, opt$num_days * opt$replications), ncol = opt$replications) -newborn.rep <- death.rep <- adult.rep <- pop.rep <- g0.rep <- g1.rep <- g2.rep <- g0a.rep <- g1a.rep <- g2a.rep <- matrix(rep(0, opt$num_days * opt$replications), ncol=opt$replications) +# Initialize matrices for results from all replications. +S0.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) +S1.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) +S2.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) +S3.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) +S4.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) +S5.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) +newborn.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) +death.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) +adult.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) +pop.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) +g0.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) +g1.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) +g2.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) +g0a.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) +g1a.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) +g2a.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) # Loop through replications. for (N.rep in 1:opt$replications) { @@ -274,7 +288,7 @@ u.d <- runif(1) if (u.d < death.prob) { death.vec <- c(death.vec, i) - } + } else { # Aggregrate index of dead bug. # Event 1 end of diapause. @@ -412,7 +426,7 @@ vec.ind[4] <- vec.ind[4] + 1 } vec.mat[i,] <- vec.ind - } + } # Event 3.3 old nymph to adult: previttelogenic or diapausing? if (vec.ind[2] == 2) { # Old nymph stage. @@ -454,22 +468,22 @@ num_insects.newborn <- length(birth.vec[,1]) vec.mat <- rbind(vec.mat, birth.vec) # Update population size for the next day. - num_insects <- num_insects - num_insects.death + num_insects.newborn + num_insects <- num_insects - num_insects.death + num_insects.newborn # Aggregate results by day. - tot.pop <- c(tot.pop, num_insects) + tot.pop <- c(tot.pop, num_insects) # Egg. - s0 <- sum(vec.mat[,2] == 0) + sum_eggs <- sum(vec.mat[,2] == 0) # Young nymph. - s1 <- sum(vec.mat[,2] == 1) + sum_young_nymphs <- sum(vec.mat[,2] == 1) # Old nymph. - s2 <- sum(vec.mat[,2] == 2) + sum_old_nymphs <- sum(vec.mat[,2] == 2) # Previtellogenic. - s3 <- sum(vec.mat[,2] == 3) + sum_previtellogenic <- sum(vec.mat[,2] == 3) # Vitellogenic. - s4 <- sum(vec.mat[,2] == 4) + sum_vitellogenic <- sum(vec.mat[,2] == 4) # Diapausing. - s5 <- sum(vec.mat[,2] == 5) + sum_diapausing <- sum(vec.mat[,2] == 5) # Overwintering adult. gen0 <- sum(vec.mat[,1] == 0) # First generation. @@ -484,12 +498,12 @@ gen1.pop[row] <- gen1 gen2.pop[row] <- gen2 - S0[row] <- s0 - S1[row] <- s1 - S2[row] <- s2 - S3[row] <- s3 - S4[row] <- s4 - S5[row] <- s5 + S0[row] <- sum_eggs + S1[row] <- sum_young_nymphs + S2[row] <- sum_old_nymphs + S3[row] <- sum_previtellogenic + S4[row] <- sum_vitellogenic + S5[row] <- sum_diapausing g0.adult[row] <- sum(vec.mat[,1] == 0) g1.adult[row] <- sum((vec.mat[,1] == 1 & vec.mat[,2] == 3) | (vec.mat[,1] == 1 & vec.mat[,2] == 4) | (vec.mat[,1] == 1 & vec.mat[,2] == 5)) @@ -521,42 +535,42 @@ g2a.rep[,N.rep] <- g2.adult } -# Mean value for adults +# Mean value for adults. mean_value_adult <- apply((S3.rep + S4.rep + S5.rep), 1, mean) -# Mean value for nymphs +# Mean value for nymphs. mean_value_nymph <- apply((S1.rep + S2.rep), 1, mean) -# Mean value for eggs +# Mean value for eggs. mean_value_egg <- apply(S0.rep, 1, mean) -# Mean value for P +# Mean value for P. g0 <- apply(g0.rep, 1, mean) -# Mean value for F1 +# Mean value for F1. g1 <- apply(g1.rep, 1, mean) -# Mean value for F2 +# Mean value for F2. g2 <- apply(g2.rep, 1, mean) -# Mean value for P adult +# Mean value for P adult. g0a <- apply(g0a.rep, 1, mean) -# Mean value for F1 adult +# Mean value for F1 adult. g1a <- apply(g1a.rep, 1, mean) -# Mean value for F2 adult +# Mean value for F2 adult. g2a <- apply(g2a.rep, 1, mean) -# Standard error for adults +# Standard error for adults. mean_value_adult.std_error <- apply((S3.rep + S4.rep + S5.rep), 1, sd) / sqrt(opt$replications) -# Standard error for nymphs +# Standard error for nymphs. mean_value_nymph.std_error <- apply((S1.rep + S2.rep) / sqrt(opt$replications), 1, sd) -# Standard error for eggs +# Standard error for eggs. mean_value_egg.std_error <- apply(S0.rep, 1, sd) / sqrt(opt$replications) -# Standard error value for P +# Standard error value for P. g0.std_error <- apply(g0.rep, 1, sd) / sqrt(opt$replications) -# Standard error for F1 +# Standard error for F1. g1.std_error <- apply(g1.rep, 1, sd) / sqrt(opt$replications) -# Standard error for F2 +# Standard error for F2. g2.std_error <- apply(g2.rep, 1, sd) / sqrt(opt$replications) -# Standard error for P adult +# Standard error for P adult. g0a.std_error <- apply(g0a.rep, 1, sd) / sqrt(opt$replications) -# Standard error for F1 adult +# Standard error for F1 adult. g1a.std_error <- apply(g1a.rep, 1, sd) / sqrt(opt$replications) -# Standard error for F2 adult +# Standard error for F2 adult. g2a.std_error <- apply(g2a.rep, 1, sd) / sqrt(opt$replications) dev.new(width=20, height=30) @@ -583,11 +597,11 @@ axis(2, cex.axis=3) leg.text <- c("Egg", "Nymph", "Adult") legend("topleft", leg.text, lty=c(1, 1, 1), col=c(4, 2, 1), cex=3) -if (opt$se_plot == 1) { +if (opt$std_error_plot == 1) { # Add Standard error lines to plot # Standard error for adults lines (day.all, mean_value_adult+mean_value_adult.std_error, lty=2) - lines (day.all, mean_value_adult-mean_value_adult.std_error, lty=2) + lines (day.all, mean_value_adult-mean_value_adult.std_error, lty=2) # Standard error for nymphs lines (day.all, mean_value_nymph+mean_value_nymph.std_error, col=2, lty=2) lines (day.all, mean_value_nymph-mean_value_nymph.std_error, col=2, lty=2) @@ -605,11 +619,11 @@ axis(2, cex.axis=3) leg.text <- c("P", "F1", "F2") legend("topleft", leg.text, lty=c(1, 1, 1), col=c(1, 2, 4), cex=3) -if (opt$se_plot == 1) { +if (opt$std_error_plot == 1) { # Add Standard error lines to plot # Standard error for adults lines (day.all, g0+g0.std_error, lty=2) - lines (day.all, g0-g0.std_error, lty=2) + lines (day.all, g0-g0.std_error, lty=2) # Standard error for nymphs lines (day.all, g1+g1.std_error, col=2, lty=2) lines (day.all, g1-g1.std_error, col=2, lty=2) @@ -631,7 +645,7 @@ # Add Standard error lines to plot # Standard error for adults lines (day.all, g0a+g0a.std_error, lty=2) - lines (day.all, g0a-g0a.std_error, lty=2) + lines (day.all, g0a-g0a.std_error, lty=2) # Standard error for nymphs lines (day.all, g1a+g1a.std_error, col=2, lty=2) lines (day.all, g1a-g1a.std_error, col=2, lty=2)