insect_phenology_model: insect_phenology

comparison insect_phenology_model.R @ 97:58bc1a2ca936 draft

Uploaded

author	greg
date	Mon, 04 Dec 2017 11:21:25 -0500
parents	f0c65c9f1078
children	d33a401c0153

comparison

equal deleted inserted replaced

-:f0c65c9f1078
+:58bc1a2ca936
 # so get the value from the first row.
 latitude <- temperature_data_frame$LATITUDE[1]
 cat("Number of days: ", opt$num_days, "\n")
-# Initialize matrices for results from all replications.
+# Initialize matrices.
-S0.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+S0.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-S1.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+S1.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-S2.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+S2.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-S3.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+S3.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-S4.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+S4.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-S5.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+S5.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-newborn.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+newborn.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-death.rep <- 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.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+adult.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-pop.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+pop.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-g0.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+g0.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-g1.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+g1.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-g2.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+g2.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-g0a.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+g0a.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-g1a.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+g1a.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-g2a.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+g2a.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
 # Loop through replications.
-for (N.rep in 1:opt$replications) {
+for (N.replications in 1:opt$replications) {
 # During each replication start with 1000 individuals.
 # TODO: user definable as well?
 num_insects <- 1000
 # Generation, Stage, degree-days, T, Diapause.
-vec.ini <- c(0, 3, 0, 0, 0)
+vector.ini <- c(0, 3, 0, 0, 0)
 # Overwintering, previttelogenic, degree-days=0, T=0, no-diapause.
-vec.mat <- rep(vec.ini, num_insects)
+vector.matrix <- rep(vector.ini, num_insects)
 # Complete matrix for the population.
-vec.mat <- base::t(matrix(vec.mat, nrow=5))
+vector.matrix <- base::t(matrix(vector.matrix, nrow=5))
 # Time series of population size.
 tot.pop <- NULL
 gen0.pop <- rep(0, opt$num_days)
 gen1.pop <- rep(0, opt$num_days)
 gen2.pop <- rep(0, opt$num_days)
 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
 # Trash bin for death.
-death.vec <- NULL
+death.vector <- NULL
 # Newborn.
-birth.vec <- NULL
+birth.vector <- NULL
 # All individuals.
 for (i in 1:num_insects) {
 # Find individual record.
-vec.ind <- vec.mat[i,]
+vector.ind <- vector.matrix[i,]
 # First of all, still alive?
 # Adjustment for late season mortality rate.
 if (latitude < 40.0) {
 post.mort <- 1
 day.kill <- 300
 }
 else {
 post.mort <- 2
 day.kill <- 250
 }
-if (vec.ind[2] == 0) {
+if (vector.ind[2] == 0) {
 # Egg.
 death.prob = opt$egg_mort * mortality.egg(mean.temp)
 }
-else if (vec.ind[2] == 1 | vec.ind[2] == 2) {
+else if (vector.ind[2] == 1 | vector.ind[2] == 2) {
 death.prob = opt$nymph_mort * mortality.nymph(mean.temp)
 }
-else if (vec.ind[2] == 3 | vec.ind[2] == 4 | vec.ind[2] == 5) {
+else if (vector.ind[2] == 3 | vector.ind[2] == 4 | vector.ind[2] == 5) {
 # For adult.
 if (doy < day.kill) {
 death.prob = opt$adult_mort * mortality.adult(mean.temp)
 }
 else {
 }
 }
 # (or dependent on temperature and life stage?)
 u.d <- runif(1)
 if (u.d < death.prob) {
-death.vec <- c(death.vec, i)
+death.vector <- c(death.vector, i)
 }
 else {
 # Aggregrate index of dead bug.
 # Event 1 end of diapause.
-if (vec.ind[1] == 0 && vec.ind[2] == 3) {
+if (vector.ind[1] == 0 && vector.ind[2] == 3) {
 # Overwintering adult (previttelogenic).
-if (photoperiod > opt$photoperiod && vec.ind[3] > 68 && doy < 180) {
+if (photoperiod > opt$photoperiod && vector.ind[3] > 68 && doy < 180) {
 # Add 68C to become fully reproductively matured.
 # Transfer to vittelogenic.
-vec.ind <- c(0, 4, 0, 0, 0)
+vector.ind <- c(0, 4, 0, 0, 0)
-vec.mat[i,] <- vec.ind
+vector.matrix[i,] <- vector.ind
 }
 else {
 # Add to degree_days.
-vec.ind[3] <- vec.ind[3] + degree_days.temp
+vector.ind[3] <- vector.ind[3] + degree_days.temp
 # Add 1 day in current stage.
-vec.ind[4] <- vec.ind[4] + 1
+vector.ind[4] <- vector.ind[4] + 1
-vec.mat[i,] <- vec.ind
+vector.matrix[i,] <- vector.ind
 }
 }
-if (vec.ind[1] != 0 && vec.ind[2] == 3) {
+if (vector.ind[1] != 0 && vector.ind[2] == 3) {
 # Not overwintering adult (previttelogenic).
-current.gen <- vec.ind[1]
+current.gen <- vector.ind[1]
-if (vec.ind[3] > 68) {
+if (vector.ind[3] > 68) {
 # Add 68C to become fully reproductively matured.
 # Transfer to vittelogenic.
-vec.ind <- c(current.gen, 4, 0, 0, 0)
+vector.ind <- c(current.gen, 4, 0, 0, 0)
-vec.mat[i,] <- vec.ind
+vector.matrix[i,] <- vector.ind
 }
 else {
 # Add to degree_days.
-vec.ind[3] <- vec.ind[3] + degree_days.temp
+vector.ind[3] <- vector.ind[3] + degree_days.temp
 # Add 1 day in current stage.
-vec.ind[4] <- vec.ind[4] + 1
+vector.ind[4] <- vector.ind[4] + 1
-vec.mat[i,] <- vec.ind
+vector.matrix[i,] <- vector.ind
 }
 }
 # Event 2 oviposition -- where population dynamics comes from.
-if (vec.ind[2] == 4 && vec.ind[1] == 0 && mean.temp > 10) {
+if (vector.ind[2] == 4 && vector.ind[1] == 0 && mean.temp > 10) {
 # Vittelogenic stage, overwintering generation.
-if (vec.ind[4] == 0) {
+if (vector.ind[4] == 0) {
 # Just turned in vittelogenic stage.
 num_insects.birth = round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size))
 }
 else {
 # Daily probability of birth.
 if (u1 < p.birth) {
 num_insects.birth = round(runif(1, 2, 8))
 }
 }
 # Add to degree_days.
-vec.ind[3] <- vec.ind[3] + degree_days.temp
+vector.ind[3] <- vector.ind[3] + degree_days.temp
 # Add 1 day in current stage.
-vec.ind[4] <- vec.ind[4] + 1
+vector.ind[4] <- vector.ind[4] + 1
-vec.mat[i,] <- vec.ind
+vector.matrix[i,] <- vector.ind
 if (num_insects.birth > 0) {
 # Add new birth -- might be in different generations.
-new.gen <- vec.ind[1] + 1
+new.gen <- vector.ind[1] + 1
 # Egg profile.
 new.ind <- c(new.gen, 0, 0, 0, 0)
-new.vec <- rep(new.ind, num_insects.birth)
+new.vector <- rep(new.ind, num_insects.birth)
 # Update batch of egg profile.
-new.vec <- t(matrix(new.vec, nrow=5))
+new.vector <- t(matrix(new.vector, nrow=5))
 # Group with total eggs laid in that day.
-birth.vec <- rbind(birth.vec, new.vec)
+birth.vector <- rbind(birth.vector, new.vector)
 }
 }
 # Event 2 oviposition -- for generation 1.
-if (vec.ind[2] == 4 && vec.ind[1] == 1 && mean.temp > 12.5 && doy < 222) {
+if (vector.ind[2] == 4 && vector.ind[1] == 1 && mean.temp > 12.5 && doy < 222) {
 # Vittelogenic stage, 1st generation
-if (vec.ind[4] == 0) {
+if (vector.ind[4] == 0) {
 # Just turned in vittelogenic stage.
 num_insects.birth=round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size))
 }
 else {
 # Daily probability of birth.
 if (u1 < p.birth) {
 num_insects.birth = round(runif(1, 2, 8))
 }
 }
 # Add to degree_days.
-vec.ind[3] <- vec.ind[3] + degree_days.temp
+vector.ind[3] <- vector.ind[3] + degree_days.temp
 # Add 1 day in current stage.
-vec.ind[4] <- vec.ind[4] + 1
+vector.ind[4] <- vector.ind[4] + 1
-vec.mat[i,] <- vec.ind
+vector.matrix[i,] <- vector.ind
 if (num_insects.birth > 0) {
 # Add new birth -- might be in different generations.
-new.gen <- vec.ind[1] + 1
+new.gen <- vector.ind[1] + 1
 # Egg profile.
 new.ind <- c(new.gen, 0, 0, 0, 0)
-new.vec <- rep(new.ind, num_insects.birth)
+new.vector <- rep(new.ind, num_insects.birth)
 # Update batch of egg profile.
-new.vec <- t(matrix(new.vec, nrow=5))
+new.vector <- t(matrix(new.vector, nrow=5))
 # Group with total eggs laid in that day.
-birth.vec <- rbind(birth.vec, new.vec)
+birth.vector <- rbind(birth.vector, new.vector)
 }
 }
 # Event 3 development (with diapause determination).
-# Event 3.1 egg development to young nymph (vec.ind[2]=0 -> egg).
+# Event 3.1 egg development to young nymph (vector.ind[2]=0 -> egg).
-if (vec.ind[2] == 0) {
+if (vector.ind[2] == 0) {
 # Egg stage.
 # Add to degree_days.
-vec.ind[3] <- vec.ind[3] + degree_days.temp
+vector.ind[3] <- vector.ind[3] + degree_days.temp
-if (vec.ind[3] >= (68 + opt$young_nymph_accum)) {
+if (vector.ind[3] >= (68 + opt$young_nymph_accum)) {
 # From egg to young nymph, degree-days requirement met.
-current.gen <- vec.ind[1]
+current.gen <- vector.ind[1]
 # Transfer to young nymph stage.
-vec.ind <- c(current.gen, 1, 0, 0, 0)
+vector.ind <- c(current.gen, 1, 0, 0, 0)
 }
 else {
 # Add 1 day in current stage.
-vec.ind[4] <- vec.ind[4] + 1
+vector.ind[4] <- vector.ind[4] + 1
 }
-vec.mat[i,] <- vec.ind
+vector.matrix[i,] <- vector.ind
 }
-# Event 3.2 young nymph to old nymph (vec.ind[2]=1 -> young nymph: determines diapause).
+# Event 3.2 young nymph to old nymph (vector.ind[2]=1 -> young nymph: determines diapause).
-if (vec.ind[2] == 1) {
+if (vector.ind[2] == 1) {
 # Young nymph stage.
 # Add to degree_days.
-vec.ind[3] <- vec.ind[3] + degree_days.temp
+vector.ind[3] <- vector.ind[3] + degree_days.temp
-if (vec.ind[3] >= (250 + opt$old_nymph_accum)) {
+if (vector.ind[3] >= (250 + opt$old_nymph_accum)) {
 # From young to old nymph, degree_days requirement met.
-current.gen <- vec.ind[1]
+current.gen <- vector.ind[1]
 # Transfer to old nym stage.
-vec.ind <- c(current.gen, 2, 0, 0, 0)
+vector.ind <- c(current.gen, 2, 0, 0, 0)
 if (photoperiod < opt$photoperiod && doy > 180) {
-vec.ind[5] <- 1
+vector.ind[5] <- 1
 } # Prepare for diapausing.
 }
 else {
 # Add 1 day in current stage.
-vec.ind[4] <- vec.ind[4] + 1
+vector.ind[4] <- vector.ind[4] + 1
 }
-vec.mat[i,] <- vec.ind
+vector.matrix[i,] <- vector.ind
 }
 # Event 3.3 old nymph to adult: previttelogenic or diapausing?
-if (vec.ind[2] == 2) {
+if (vector.ind[2] == 2) {
 # Old nymph stage.
 # Add to degree_days.
-vec.ind[3] <- vec.ind[3] + degree_days.temp
+vector.ind[3] <- vector.ind[3] + degree_days.temp
-if (vec.ind[3] >= (200 + opt$adult_accum)) {
+if (vector.ind[3] >= (200 + opt$adult_accum)) {
 # From old to adult, degree_days requirement met.
-current.gen <- vec.ind[1]
+current.gen <- vector.ind[1]
-if (vec.ind[5] == 0) {
+if (vector.ind[5] == 0) {
 # Non-diapausing adult -- previttelogenic.
-vec.ind <- c(current.gen, 3, 0, 0, 0)
+vector.ind <- c(current.gen, 3, 0, 0, 0)
 }
 else {
 # Diapausing.
-vec.ind <- c(current.gen, 5, 0, 0, 1)
+vector.ind <- c(current.gen, 5, 0, 0, 1)
 }
 }
 else {
 # Add 1 day in current stage.
-vec.ind[4] <- vec.ind[4] + 1
+vector.ind[4] <- vector.ind[4] + 1
 }
-vec.mat[i,] <- vec.ind
+vector.matrix[i,] <- vector.ind
 }
 # Event 4 growing of diapausing adult (unimportant, but still necessary).
-if (vec.ind[2] == 5) {
+if (vector.ind[2] == 5) {
-vec.ind[3] <- vec.ind[3] + degree_days.temp
+vector.ind[3] <- vector.ind[3] + degree_days.temp
-vec.ind[4] <- vec.ind[4] + 1
+vector.ind[4] <- vector.ind[4] + 1
-vec.mat[i,] <- vec.ind
+vector.matrix[i,] <- vector.ind
 }
 } # Else if it is still alive.
 } # End of the individual bug loop.
 # Find how many died.
-num_insects.death <- length(death.vec)
+num_insects.death <- length(death.vector)
 if (num_insects.death > 0) {
-vec.mat <- vec.mat[-death.vec, ]
+vector.matrix <- vector.matrix[-death.vector, ]
 }
 # Remove record of dead.
 # Find how many new born.
-num_insects.newborn <- length(birth.vec[,1])
+num_insects.newborn <- length(birth.vector[,1])
-vec.mat <- rbind(vec.mat, birth.vec)
+vector.matrix <- rbind(vector.matrix, birth.vector)
 # Update population size for the next day.
 num_insects <- num_insects - num_insects.death + num_insects.newborn
 # Aggregate results by day.
 tot.pop <- c(tot.pop, num_insects)
 # Egg.
-sum_eggs <- sum(vec.mat[,2] == 0)
+sum_eggs <- sum(vector.matrix[,2] == 0)
 # Young nymph.
-sum_young_nymphs <- sum(vec.mat[,2] == 1)
+sum_young_nymphs <- sum(vector.matrix[,2] == 1)
 # Old nymph.
-sum_old_nymphs <- sum(vec.mat[,2] == 2)
+sum_old_nymphs <- sum(vector.matrix[,2] == 2)
 # Previtellogenic.
-sum_previtellogenic <- sum(vec.mat[,2] == 3)
+sum_previtellogenic <- sum(vector.matrix[,2] == 3)
 # Vitellogenic.
-sum_vitellogenic <- sum(vec.mat[,2] == 4)
+sum_vitellogenic <- sum(vector.matrix[,2] == 4)
 # Diapausing.
-sum_diapausing <- sum(vec.mat[,2] == 5)
+sum_diapausing <- sum(vector.matrix[,2] == 5)
 # Overwintering adult.
-gen0 <- sum(vec.mat[,1] == 0)
+gen0 <- sum(vector.matrix[,1] == 0)
 # First generation.
-gen1 <- sum(vec.mat[,1] == 1)
+gen1 <- sum(vector.matrix[,1] == 1)
 # Second generation.
-gen2 <- sum(vec.mat[,1] == 2)
+gen2 <- sum(vector.matrix[,1] == 2)
 # Sum of all adults.
-num_insects.adult <- sum(vec.mat[,2] == 3) + sum(vec.mat[,2] == 4) + sum(vec.mat[,2] == 5)
+num_insects.adult <- sum(vector.matrix[,2] == 3) + sum(vector.matrix[,2] == 4) + sum(vector.matrix[,2] == 5)
 # Generation 0 pop size.
 gen0.pop[row] <- gen0
 gen1.pop[row] <- gen1
 gen2.pop[row] <- gen2
 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)
+g0.adult[row] <- sum(vector.matrix[,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))
+g1.adult[row] <- sum((vector.matrix[,1] == 1 & vector.matrix[,2] == 3) | (vector.matrix[,1] == 1 & vector.matrix[,2] == 4) | (vector.matrix[,1] == 1 & vector.matrix[,2] == 5))
-g2.adult[row] <- sum((vec.mat[,1]== 2 & vec.mat[,2] == 3) | (vec.mat[,1] == 2 & vec.mat[,2] == 4) | (vec.mat[,1] == 2 & vec.mat[,2] == 5))
+g2.adult[row] <- sum((vector.matrix[,1]== 2 & vector.matrix[,2] == 3) | (vector.matrix[,1] == 2 & vector.matrix[,2] == 4) | (vector.matrix[,1] == 2 & vector.matrix[,2] == 5))
 N.newborn[row] <- num_insects.newborn
 N.death[row] <- num_insects.death
 N.adult[row] <- num_insects.adult
 }   # end of days specified in the input temperature data
 degree_days.cum <- cumsum(degree_days.day)
 # Collect all the outputs.
-S0.rep[,N.rep] <- S0
+S0.replications[,N.replications] <- S0
-S1.rep[,N.rep] <- S1
+S1.replications[,N.replications] <- S1
-S2.rep[,N.rep] <- S2
+S2.replications[,N.replications] <- S2
-S3.rep[,N.rep] <- S3
+S3.replications[,N.replications] <- S3
-S4.rep[,N.rep] <- S4
+S4.replications[,N.replications] <- S4
-S5.rep[,N.rep] <- S5
+S5.replications[,N.replications] <- S5
-newborn.rep[,N.rep] <- N.newborn
+newborn.replications[,N.replications] <- N.newborn
-death.rep[,N.rep] <- N.death
+death.replications[,N.replications] <- N.death
-adult.rep[,N.rep] <- N.adult
+adult.replications[,N.replications] <- N.adult
-pop.rep[,N.rep] <- tot.pop
+pop.replications[,N.replications] <- tot.pop
-g0.rep[,N.rep] <- gen0.pop
+g0.replications[,N.replications] <- gen0.pop
-g1.rep[,N.rep] <- gen1.pop
+g1.replications[,N.replications] <- gen1.pop
-g2.rep[,N.rep] <- gen2.pop
+g2.replications[,N.replications] <- gen2.pop
-g0a.rep[,N.rep] <- g0.adult
+g0a.replications[,N.replications] <- g0.adult
-g1a.rep[,N.rep] <- g1.adult
+g1a.replications[,N.replications] <- g1.adult
-g2a.rep[,N.rep] <- g2.adult
+g2a.replications[,N.replications] <- g2.adult
 }
 # Mean value for adults.
-mean_value_adult <- apply((S3.rep + S4.rep + S5.rep), 1, mean)
+adult_mean <- apply((S3.replications + S4.replications + S5.replications), 1, mean)
 # Mean value for nymphs.
-mean_value_nymph <- apply((S1.rep + S2.rep), 1, mean)
+nymph_mean <- apply((S1.replications + S2.replications), 1, mean)
 # Mean value for eggs.
-mean_value_egg <- apply(S0.rep, 1, mean)
+egg_mean <- apply(S0.replications, 1, mean)
 # Mean value for P.
-g0 <- apply(g0.rep, 1, mean)
+g0 <- apply(g0.replications, 1, mean)
 # Mean value for F1.
-g1 <- apply(g1.rep, 1, mean)
+g1 <- apply(g1.replications, 1, mean)
 # Mean value for F2.
-g2 <- apply(g2.rep, 1, mean)
+g2 <- apply(g2.replications, 1, mean)
 # Mean value for P adult.
-g0a <- apply(g0a.rep, 1, mean)
+g0a <- apply(g0a.replications, 1, mean)
 # Mean value for F1 adult.
-g1a <- apply(g1a.rep, 1, mean)
+g1a <- apply(g1a.replications, 1, mean)
 # Mean value for F2 adult.
-g2a <- apply(g2a.rep, 1, mean)
+g2a <- apply(g2a.replications, 1, mean)
 # Standard error for adults.
-mean_value_adult.std_error <- apply((S3.rep + S4.rep + S5.rep), 1, sd) / sqrt(opt$replications)
+adult_mean.std_error <- apply((S3.replications + S4.replications + S5.replications), 1, sd) / sqrt(opt$replications)
 # Standard error for nymphs.
-mean_value_nymph.std_error <- apply((S1.rep + S2.rep) / sqrt(opt$replications), 1, sd)
+nymph_mean.std_error <- apply((S1.replications + S2.replications) / sqrt(opt$replications), 1, sd)
 # Standard error for eggs.
-mean_value_egg.std_error <- apply(S0.rep, 1, sd) / sqrt(opt$replications)
+egg_mean.std_error <- apply(S0.replications, 1, sd) / sqrt(opt$replications)
 # Standard error value for P.
-g0.std_error <- apply(g0.rep, 1, sd) / sqrt(opt$replications)
+g0.std_error <- apply(g0.replications, 1, sd) / sqrt(opt$replications)
 # Standard error for F1.
-g1.std_error <- apply(g1.rep, 1, sd) / sqrt(opt$replications)
+g1.std_error <- apply(g1.replications, 1, sd) / sqrt(opt$replications)
 # Standard error for F2.
-g2.std_error <- apply(g2.rep, 1, sd) / sqrt(opt$replications)
+g2.std_error <- apply(g2.replications, 1, sd) / sqrt(opt$replications)
 # Standard error for P adult.
-g0a.std_error <- apply(g0a.rep, 1, sd) / sqrt(opt$replications)
+g0a.std_error <- apply(g0a.replications, 1, sd) / sqrt(opt$replications)
 # Standard error for F1 adult.
-g1a.std_error <- apply(g1a.rep, 1, sd) / sqrt(opt$replications)
+g1a.std_error <- apply(g1a.replications, 1, sd) / sqrt(opt$replications)
 # Standard error for F2 adult.
-g2a.std_error <- apply(g2a.rep, 1, sd) / sqrt(opt$replications)
+g2a.std_error <- apply(g2a.replications, 1, sd) / sqrt(opt$replications)
 dev.new(width=20, height=30)
 # Start PDF device driver to save charts to output.
 pdf(file=opt$output, width=20, height=30, bg="white")
 start_date <- temperature_data_frame$DATE[1]
 end_date <- temperature_data_frame$DATE[opt$num_days]
 # Subfigure 1: population size by life stage
 title <- paste(opt$insect, ": Total pop. by life stage :", opt$location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ")
-plot(day.all, mean_value_adult, main=title, type="l", ylim=c(0, max(mean_value_egg + mean_value_egg.std_error, mean_value_nymph + mean_value_nymph.std_error, mean_value_adult + mean_value_adult.std_error)), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3)
+plot(day.all, adult_mean, main=title, type="l", ylim=c(0, max(egg_mean + egg_mean.std_error, nymph_mean + nymph_mean.std_error, adult_mean + adult_mean.std_error)), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3)
 # Young and old nymphs.
-lines(day.all, mean_value_nymph, lwd=2, lty=1, col=2)
+lines(day.all, nymph_mean, lwd=2, lty=1, col=2)
 # Eggs
-lines(day.all, mean_value_egg, lwd=2, lty=1, col=4)
+lines(day.all, egg_mean, 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"))
 axis(2, cex.axis=3)
-leg.text <- c("Egg", "Nymph", "Adult")
+legend("topleft", c("Egg", "Nymph", "Adult"), lty=c(1, 1, 1), col=c(4, 2, 1), cex=3)
-legend("topleft", leg.text, lty=c(1, 1, 1), col=c(4, 2, 1), cex=3)
 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, adult_mean+adult_mean.std_error, lty=2)
-lines (day.all, mean_value_adult-mean_value_adult.std_error, lty=2)
+lines (day.all, adult_mean-adult_mean.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, nymph_mean+nymph_mean.std_error, col=2, lty=2)
-lines (day.all, mean_value_nymph-mean_value_nymph.std_error, col=2, lty=2)
+lines (day.all, nymph_mean-nymph_mean.std_error, col=2, lty=2)
 # Standard error for eggs
-lines (day.all, mean_value_egg+mean_value_egg.std_error, col=4, lty=2)
+lines (day.all, egg_mean+egg_mean.std_error, col=4, lty=2)
-lines (day.all, mean_value_egg-mean_value_egg.std_error, col=4, lty=2)
+lines (day.all, egg_mean-egg_mean.std_error, col=4, lty=2)
 }
 # Subfigure 2: population size by generation
 title <- paste(opt$insect, ": Total pop. by generation :", opt$location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ")
 plot(day.all, g0, main=title, type="l", ylim=c(0, max(g2)), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3)
 lines(day.all, g1, lwd = 2, lty = 1, col=2)
 lines(day.all, g2, 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"))
 axis(2, cex.axis=3)
-leg.text <- c("P", "F1", "F2")
+legend("topleft", c("P", "F1", "F2"), lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)
-legend("topleft", leg.text, lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)
 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)
 plot(day.all, g0a, ylim=c(0, max(g2a) + 100), main=title, type="l", axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3)
 lines(day.all, g1a, lwd = 2, lty = 1, col=2)
 lines(day.all, g2a, 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"))
 axis(2, cex.axis=3)
-leg.text <- c("P", "F1", "F2")
+legend("topleft", c("P", "F1", "F2"), lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)
-legend("topleft", leg.text, lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)
 if (opt$std_error_plot == 1) {
 # 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)

Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 97:58bc1a2ca936 draft