Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 98:d33a401c0153 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Mon, 04 Dec 2017 11:54:10 -0500
parents 58bc1a2ca936
children 51045fde125e
comparison
equal deleted inserted replaced
97:58bc1a2ca936 98:d33a401c0153
233 # Time series of population size. 233 # Time series of population size.
234 tot.pop <- NULL 234 tot.pop <- NULL
235 gen0.pop <- rep(0, opt$num_days) 235 gen0.pop <- rep(0, opt$num_days)
236 gen1.pop <- rep(0, opt$num_days) 236 gen1.pop <- rep(0, opt$num_days)
237 gen2.pop <- rep(0, opt$num_days) 237 gen2.pop <- rep(0, opt$num_days)
238 S0 <- S1 <- S2 <- S3 <- S4 <- S5 <- rep(0, opt$num_days) 238 S0 <- rep(0, opt$num_days)
239 g0.adult <- g1.adult <- g2.adult <- rep(0, opt$num_days) 239 S1 <- rep(0, opt$num_days)
240 N.newborn <- N.death <- N.adult <- rep(0, opt$num_days) 240 S2 <- rep(0, opt$num_days)
241 S3 <- rep(0, opt$num_days)
242 S4 <- rep(0, opt$num_days)
243 S5 <- rep(0, opt$num_days)
244 g0.adult <- rep(0, opt$num_days)
245 g1.adult <- rep(0, opt$num_days)
246 g2.adult <- rep(0, opt$num_days)
247 N.newborn <- rep(0, opt$num_days)
248 N.death <- rep(0, opt$num_days)
249 N.adult <- rep(0, opt$num_days)
241 degree_days.day <- rep(0, opt$num_days) 250 degree_days.day <- rep(0, opt$num_days)
242 # All the days included in the input temperature dataset. 251 # All the days included in the input temperature dataset.
243 for (row in 1:opt$num_days) { 252 for (row in 1:opt$num_days) {
244 # Get the integer day of the year for the current row. 253 # Get the integer day of the year for the current row.
245 doy <- temperature_data_frame$DOY[row] 254 doy <- temperature_data_frame$DOY[row]
483 # Vitellogenic. 492 # Vitellogenic.
484 sum_vitellogenic <- sum(vector.matrix[,2] == 4) 493 sum_vitellogenic <- sum(vector.matrix[,2] == 4)
485 # Diapausing. 494 # Diapausing.
486 sum_diapausing <- sum(vector.matrix[,2] == 5) 495 sum_diapausing <- sum(vector.matrix[,2] == 5)
487 # Overwintering adult. 496 # Overwintering adult.
488 gen0 <- sum(vector.matrix[,1] == 0) 497 sum_overwintering_adults <- sum(vector.matrix[,1] == 0)
489 # First generation. 498 # First generation.
490 gen1 <- sum(vector.matrix[,1] == 1) 499 sum_first_generation <- sum(vector.matrix[,1] == 1)
491 # Second generation. 500 # Second generation.
492 gen2 <- sum(vector.matrix[,1] == 2) 501 sum_second_generation <- sum(vector.matrix[,1] == 2)
493 # Sum of all adults. 502 # Sum of all adults.
494 num_insects.adult <- sum(vector.matrix[,2] == 3) + sum(vector.matrix[,2] == 4) + sum(vector.matrix[,2] == 5) 503 num_insects.adult <- sum(vector.matrix[,2] == 3) + sum(vector.matrix[,2] == 4) + sum(vector.matrix[,2] == 5)
495 504
496 # Generation 0 pop size. 505 # Population sizes.
497 gen0.pop[row] <- gen0 506 gen0.pop[row] <- sum_overwintering_adults
498 gen1.pop[row] <- gen1 507 gen1.pop[row] <- sum_first_generation
499 gen2.pop[row] <- gen2 508 gen2.pop[row] <- sum_second_generation
500 509
501 S0[row] <- sum_eggs 510 S0[row] <- sum_eggs
502 S1[row] <- sum_young_nymphs 511 S1[row] <- sum_young_nymphs
503 S2[row] <- sum_old_nymphs 512 S2[row] <- sum_old_nymphs
504 S3[row] <- sum_previtellogenic 513 S3[row] <- sum_previtellogenic
550 # Mean value for P adult. 559 # Mean value for P adult.
551 g0a <- apply(g0a.replications, 1, mean) 560 g0a <- apply(g0a.replications, 1, mean)
552 # Mean value for F1 adult. 561 # Mean value for F1 adult.
553 g1a <- apply(g1a.replications, 1, mean) 562 g1a <- apply(g1a.replications, 1, mean)
554 # Mean value for F2 adult. 563 # Mean value for F2 adult.
555 g2a <- apply(g2a.replications, 1, mean) 564 F2_adult_mean <- apply(g2a.replications, 1, mean)
556 565
557 # Standard error for adults. 566 # Standard error for adults.
558 adult_mean.std_error <- apply((S3.replications + S4.replications + S5.replications), 1, sd) / sqrt(opt$replications) 567 adult_mean.std_error <- apply((S3.replications + S4.replications + S5.replications), 1, sd) / sqrt(opt$replications)
559 # Standard error for nymphs. 568 # Standard error for nymphs.
560 nymph_mean.std_error <- apply((S1.replications + S2.replications) / sqrt(opt$replications), 1, sd) 569 nymph_mean.std_error <- apply((S1.replications + S2.replications) / sqrt(opt$replications), 1, sd)
630 lines (day.all, g2-g2.std_error, col=4, lty=2) 639 lines (day.all, g2-g2.std_error, col=4, lty=2)
631 } 640 }
632 641
633 # Subfigure 3: adult population size by generation 642 # Subfigure 3: adult population size by generation
634 title <- paste(opt$insect, ": Adult pop. by generation :", opt$location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ") 643 title <- paste(opt$insect, ": Adult pop. by generation :", opt$location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ")
635 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) 644 plot(day.all, g0a, ylim=c(0, max(F2_adult_mean) + 100), main=title, type="l", axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3)
636 lines(day.all, g1a, lwd = 2, lty = 1, col=2) 645 lines(day.all, g1a, lwd = 2, lty = 1, col=2)
637 lines(day.all, g2a, lwd = 2, lty = 1, col=4) 646 lines(day.all, F2_adult_mean, lwd = 2, lty = 1, col=4)
638 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")) 647 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"))
639 axis(2, cex.axis=3) 648 axis(2, cex.axis=3)
640 legend("topleft", c("P", "F1", "F2"), lty=c(1, 1, 1), col=c(1, 2, 4), cex=3) 649 legend("topleft", c("P", "F1", "F2"), lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)
641 if (opt$std_error_plot == 1) { 650 if (opt$std_error_plot == 1) {
642 # Add Standard error lines to plot 651 # Add Standard error lines to plot
645 lines (day.all, g0a-g0a.std_error, lty=2) 654 lines (day.all, g0a-g0a.std_error, lty=2)
646 # Standard error for nymphs 655 # Standard error for nymphs
647 lines (day.all, g1a+g1a.std_error, col=2, lty=2) 656 lines (day.all, g1a+g1a.std_error, col=2, lty=2)
648 lines (day.all, g1a-g1a.std_error, col=2, lty=2) 657 lines (day.all, g1a-g1a.std_error, col=2, lty=2)
649 # Standard error for eggs 658 # Standard error for eggs
650 lines (day.all, g2a+g2a.std_error, col=4, lty=2) 659 lines (day.all, F2_adult_mean+g2a.std_error, col=4, lty=2)
651 lines (day.all, g2a-g2a.std_error, col=4, lty=2) 660 lines (day.all, F2_adult_mean-g2a.std_error, col=4, lty=2)
652 } 661 }
653 662
654 # Turn off device driver to flush output. 663 # Turn off device driver to flush output.
655 dev.off() 664 dev.off()