Mercurial > repos > greg > insect_phenology_model
changeset 103:4d9e40d8af0f draft
Uploaded
| author | greg |
|---|---|
| date | Tue, 05 Dec 2017 11:01:23 -0500 |
| parents | ce996e90f5a7 |
| children | d6e37828b094 |
| files | insect_phenology_model.R |
| diffstat | 1 files changed, 86 insertions(+), 67 deletions(-) [+] |
line wrap: on
line diff
--- a/insect_phenology_model.R Tue Dec 05 10:32:04 2017 -0500 +++ b/insect_phenology_model.R Tue Dec 05 11:01:23 2017 -0500 @@ -83,9 +83,9 @@ # Mean temperature for current row. curr_mean_temp <- 0.5 * (curr_min_temp + curr_max_temp) # Initialize degree day accumulation - degree_days <- 0 + averages <- 0 if (curr_max_temp < threshold) { - degree_days <- 0 + averages <- 0 } else { # Initialize hourly temperature. @@ -138,9 +138,9 @@ } } } - degree_days <- sum(dh) / 24 + averages <- sum(dh) / 24 } - return(c(curr_mean_temp, degree_days)) + return(c(curr_mean_temp, averages)) } mortality.adult = function(temperature) { @@ -197,15 +197,19 @@ group1, group2, group3, group1_std_error, group2_std_error, group3_std_error) { if (chart_type == "pop_size_by_life_stage") { title <- paste(insect, ": Total pop. by life stage :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ") - legend("topleft", c("Egg", "Nymph", "Adult"), lty=c(1, 1, 1), col=c(4, 2, 1), cex=3) + legend_text <- c("Egg", "Nymph", "Adult") + columns <- c(4, 2, 1) } else if (chart_type == "pop_size_by_generation") { title <- paste(insect, ": Total pop. by generation :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ") - legend("topleft", c("P", "F1", "F2"), lty=c(1, 1, 1), col=c(1, 2, 4), cex=3) + legend_text <- c("P", "F1", "F2") + columns <- col=c(1, 2, 4) } else if (chart_type == "adult_pop_size_by_generation") { title <- paste(insect, ": Adult pop. by generation :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ") - legend("topleft", c("P", "F1", "F2"), lty=c(1, 1, 1), col=c(1, 2, 4), cex=3) + legend_text <- c("P", "F1", "F2") + columns <- col=c(1, 2, 4) } plot(days, group1, main=title, type="l", ylim=c(0, maxval), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3) + legend("topleft", legend_text, lty=c(1, 1, 1), columns, cex=3) lines(days, group2, lwd=2, lty=1, col=2) lines(days, group3, lwd=2, lty=1, col=4) axis(1, at=c(1:12) * 30 - 15, cex.axis=3, labels=c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")) @@ -237,17 +241,20 @@ Previtellogenic.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) Vitellogenic.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) Diapausing.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) + newborn.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) +adult.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) death.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) -adult.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) -population.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) + P.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) +P_adults.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) F1.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) +F1_adults.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) F2.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) -P_adults.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) -F1_adults.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) F2_adults.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) +population.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) + # Process replications. for (N.replications in 1:opt$replications) { # During each replication start with 1000 individuals. @@ -260,23 +267,28 @@ # Complete matrix for the population. vector.matrix <- base::t(matrix(vector.matrix, nrow=5)) # Time series of population size. - total.population <- NULL - overwintering_adult.population <- rep(0, opt$num_days) - first_generation.population <- rep(0, opt$num_days) - second_generation.population <- rep(0, opt$num_days) Eggs <- rep(0, opt$num_days) YoungNymphs <- rep(0, opt$num_days) OldNymphs <- rep(0, opt$num_days) Previtellogenic <- rep(0, opt$num_days) Vitellogenic <- rep(0, opt$num_days) Diapausing <- rep(0, opt$num_days) + + N.newborn <- rep(0, opt$num_days) + N.adult <- rep(0, opt$num_days) + N.death <- rep(0, opt$num_days) + + overwintering_adult.population <- rep(0, opt$num_days) + first_generation.population <- rep(0, opt$num_days) + second_generation.population <- rep(0, opt$num_days) + P.adult <- rep(0, opt$num_days) F1.adult <- rep(0, opt$num_days) F2.adult <- rep(0, opt$num_days) - N.newborn <- rep(0, opt$num_days) - N.death <- rep(0, opt$num_days) - N.adult <- rep(0, opt$num_days) - degree_days.day <- rep(0, opt$num_days) + + total.population <- NULL + + averages.day <- rep(0, opt$num_days) # All the days included in the input temperature dataset. for (row in 1:opt$num_days) { # Get the integer day of the year for the current row. @@ -285,18 +297,17 @@ photoperiod <- temperature_data_frame$DAYLEN[row] temp.profile <- get_temperature_at_hour(latitude, temperature_data_frame, row, opt$num_days) mean.temp <- temp.profile[1] - degree_days.temp <- temp.profile[2] - degree_days.day[row] <- degree_days.temp + averages.temp <- temp.profile[2] + averages.day[row] <- averages.temp # Trash bin for death. death.vector <- NULL # Newborn. birth.vector <- NULL # All individuals. for (i in 1:num_insects) { - # Find individual record. + # Individual record. vector.individual <- vector.matrix[i,] - # First of all, still alive? - # Adjustment for late season mortality rate. + # Adjustment for late season mortality rate (still alive?). if (latitude < 40.0) { post.mortality <- 1 day.kill <- 300 @@ -313,7 +324,7 @@ death.probability = opt$nymph_mortality * mortality.nymph(mean.temp) } else if (vector.individual[2] == 3 | vector.individual[2] == 4 | vector.individual[2] == 5) { - # For adult. + # Adult. if (doy < day.kill) { death.probability = opt$adult_mortality * mortality.adult(mean.temp) } @@ -322,14 +333,14 @@ death.probability = opt$adult_mortality * post.mortality * mortality.adult(mean.temp) } } - # (or dependent on temperature and life stage?) + # Dependent on temperature and life stage? u.d <- runif(1) if (u.d < death.probability) { death.vector <- c(death.vector, i) } else { # Aggregrate index of dead bug. - # Event 1 end of diapause. + # End of diapause. if (vector.individual[1] == 0 && vector.individual[2] == 3) { # Overwintering adult (previttelogenic). if (photoperiod > opt$photoperiod && vector.individual[3] > 68 && doy < 180) { @@ -339,8 +350,8 @@ vector.matrix[i,] <- vector.individual } else { - # Add to degree_days. - vector.individual[3] <- vector.individual[3] + degree_days.temp + # Add to # Add average temperature for current day.. + vector.individual[3] <- vector.individual[3] + averages.temp # Add 1 day in current stage. vector.individual[4] <- vector.individual[4] + 1 vector.matrix[i,] <- vector.individual @@ -356,8 +367,8 @@ vector.matrix[i,] <- vector.individual } else { - # Add to degree_days. - vector.individual[3] <- vector.individual[3] + degree_days.temp + # Add average temperature for current day. + vector.individual[3] <- vector.individual[3] + averages.temp # Add 1 day in current stage. vector.individual[4] <- vector.individual[4] + 1 vector.matrix[i,] <- vector.individual @@ -378,8 +389,8 @@ num_insects.birth = round(runif(1, 2, 8)) } } - # Add to degree_days. - vector.individual[3] <- vector.individual[3] + degree_days.temp + # Add average temperature for current day. + vector.individual[3] <- vector.individual[3] + averages.temp # Add 1 day in current stage. vector.individual[4] <- vector.individual[4] + 1 vector.matrix[i,] <- vector.individual @@ -410,8 +421,8 @@ num_insects.birth = round(runif(1, 2, 8)) } } - # Add to degree_days. - vector.individual[3] <- vector.individual[3] + degree_days.temp + # Add average temperature for current day. + vector.individual[3] <- vector.individual[3] + averages.temp # Add 1 day in current stage. vector.individual[4] <- vector.individual[4] + 1 vector.matrix[i,] <- vector.individual @@ -431,8 +442,8 @@ # Event 3.1 egg development to young nymph (vector.individual[2]=0 -> egg). if (vector.individual[2] == 0) { # Egg stage. - # Add to degree_days. - vector.individual[3] <- vector.individual[3] + degree_days.temp + # Add average temperature for current day. + vector.individual[3] <- vector.individual[3] + averages.temp if (vector.individual[3] >= (68+opt$young_nymph_accumulation)) { # From egg to young nymph, degree-days requirement met. current.gen <- vector.individual[1] @@ -448,8 +459,8 @@ # Event 3.2 young nymph to old nymph (vector.individual[2]=1 -> young nymph: determines diapause). if (vector.individual[2] == 1) { # Young nymph stage. - # Add to degree_days. - vector.individual[3] <- vector.individual[3] + degree_days.temp + # Add average temperature for current day. + vector.individual[3] <- vector.individual[3] + averages.temp if (vector.individual[3] >= (250+opt$old_nymph_accumulation)) { # From young to old nymph, degree_days requirement met. current.gen <- vector.individual[1] @@ -468,8 +479,8 @@ # Event 3.3 old nymph to adult: previttelogenic or diapausing? if (vector.individual[2] == 2) { # Old nymph stage. - # Add to degree_days. - vector.individual[3] <- vector.individual[3] + degree_days.temp + # Add average temperature for current day. + vector.individual[3] <- vector.individual[3] + averages.temp if (vector.individual[3] >= (200+opt$adult_accumulation)) { # From old to adult, degree_days requirement met. current.gen <- vector.individual[1] @@ -490,7 +501,7 @@ } # Event 4 growing of diapausing adult (unimportant, but still necessary). if (vector.individual[2] == 5) { - vector.individual[3] <- vector.individual[3] + degree_days.temp + vector.individual[3] <- vector.individual[3] + averages.temp vector.individual[4] <- vector.individual[4] + 1 vector.matrix[i,] <- vector.individual } @@ -544,7 +555,7 @@ N.adult[row] <- num_insects.adult } # End of days specified in the input temperature data. - degree_days.cum <- cumsum(degree_days.day) + averages.cum <- cumsum(averages.day) # Define the output values. Eggs.replications[,N.replications] <- Eggs @@ -565,40 +576,48 @@ F2_adults.replications[,N.replications] <- F2.adult } -# Mean value for adults. -adults <- apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, mean) -# Mean value for nymphs. -nymphs <- apply((YoungNymphs.replications+OldNymphs.replications), 1, mean) # Mean value for eggs. eggs <- apply(Eggs.replications, 1, mean) +# Standard error for eggs. +eggs.std_error <- apply(Eggs.replications, 1, sd) / sqrt(opt$replications) + +# Mean value for nymphs. +nymphs <- apply((YoungNymphs.replications+OldNymphs.replications), 1, mean) +# Standard error for nymphs. +nymphs.std_error <- apply((YoungNymphs.replications+OldNymphs.replications) / sqrt(opt$replications), 1, sd) + +# Mean value for adults. +adults <- apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, mean) +# Standard error for adults. +adults.std_error <- apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, sd) / sqrt(opt$replications) + # Mean value for P. P <- apply(P.replications, 1, mean) +# Standard error for P. +P.std_error <- apply(P.replications, 1, sd) / sqrt(opt$replications) + +# Mean value for P adults. +P_adults <- apply(P_adults.replications, 1, mean) +# Standard error for P_adult. +P_adults.std_error <- apply(P_adults.replications, 1, sd) / sqrt(opt$replications) + # Mean value for F1. F1 <- apply(F1.replications, 1, mean) -# Mean value for F2. -F2 <- apply(F2.replications, 1, mean) -# Mean value for P adults. -P_adults <- apply(P_adults.replications, 1, mean) +# Standard error for F1. +F1.std_error <- apply(F1.replications, 1, sd) / sqrt(opt$replications) + # Mean value for F1 adults. F1_adults <- apply(F1_adults.replications, 1, mean) +# Standard error for F1 adult. +F1_adults.std_error <- apply(F1_adults.replications, 1, sd) / sqrt(opt$replications) + +# Mean value for F2. +F2 <- apply(F2.replications, 1, mean) +# Standard error for F2. +F2.std_error <- apply(F2.replications, 1, sd) / sqrt(opt$replications) + # Mean value for F2 adults. F2_adults <- apply(F2_adults.replications, 1, mean) -# Standard error for adults. -adults.std_error <- apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, sd) / sqrt(opt$replications) -# Standard error for nymphs. -nymphs.std_error <- apply((YoungNymphs.replications+OldNymphs.replications) / sqrt(opt$replications), 1, sd) -# Standard error for eggs. -eggs.std_error <- apply(Eggs.replications, 1, sd) / sqrt(opt$replications) -# Standard error value for P. -P.std_error <- apply(P.replications, 1, sd) / sqrt(opt$replications) -# Standard error for F1. -F1.std_error <- apply(F1.replications, 1, sd) / sqrt(opt$replications) -# Standard error for F2. -F2.std_error <- apply(F2.replications, 1, sd) / sqrt(opt$replications) -# Standard error for P adult. -P_adults.std_error <- apply(P_adults.replications, 1, sd) / sqrt(opt$replications) -# Standard error for F1 adult. -F1_adults.std_error <- apply(F1_adults.replications, 1, sd) / sqrt(opt$replications) # Standard error for F2 adult. F2_adults.std_error <- apply(F2_adults.replications, 1, sd) / sqrt(opt$replications)