Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 100:52114847afff 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 13:21:20 -0500
parents 51045fde125e
children 691b677d8f83
comparison
equal deleted inserted replaced
99:51045fde125e 100:52114847afff
210 S5.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 210 S5.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
211 newborn.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 211 newborn.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
212 death.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 212 death.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
213 adult.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 213 adult.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
214 pop.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 214 pop.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
215 g0.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 215 P.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
216 g1.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 216 F1.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
217 g2.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 217 F2.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
218 g0a.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 218 P_adults.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
219 g1a.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 219 F1_adults.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
220 F2_adults.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 220 F2_adults.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
221 221
222 # Loop through replications. 222 # Loop through replications.
223 for (N.replications in 1:opt$replications) { 223 for (N.replications in 1:opt$replications) {
224 # During each replication start with 1000 individuals. 224 # During each replication start with 1000 individuals.
239 S1 <- rep(0, opt$num_days) 239 S1 <- rep(0, opt$num_days)
240 S2 <- rep(0, opt$num_days) 240 S2 <- rep(0, opt$num_days)
241 S3 <- rep(0, opt$num_days) 241 S3 <- rep(0, opt$num_days)
242 S4 <- rep(0, opt$num_days) 242 S4 <- rep(0, opt$num_days)
243 S5 <- rep(0, opt$num_days) 243 S5 <- rep(0, opt$num_days)
244 g0.adult <- rep(0, opt$num_days) 244 P.adult <- rep(0, opt$num_days)
245 g1.adult <- rep(0, opt$num_days) 245 F1.adult <- rep(0, opt$num_days)
246 g2.adult <- rep(0, opt$num_days) 246 F2.adult <- rep(0, opt$num_days)
247 N.newborn <- rep(0, opt$num_days) 247 N.newborn <- rep(0, opt$num_days)
248 N.death <- rep(0, opt$num_days) 248 N.death <- rep(0, opt$num_days)
249 N.adult <- rep(0, opt$num_days) 249 N.adult <- rep(0, opt$num_days)
250 degree_days.day <- rep(0, opt$num_days) 250 degree_days.day <- rep(0, opt$num_days)
251 # All the days included in the input temperature dataset. 251 # All the days included in the input temperature dataset.
503 # Vitellogenic population size. 503 # Vitellogenic population size.
504 S4[row] <- sum(vector.matrix[,2] == 4) 504 S4[row] <- sum(vector.matrix[,2] == 4)
505 # Diapausing population size. 505 # Diapausing population size.
506 S5[row] <- sum(vector.matrix[,2] == 5) 506 S5[row] <- sum(vector.matrix[,2] == 5)
507 507
508 g0.adult[row] <- sum(vector.matrix[,1] == 0) 508 P.adult[row] <- sum(vector.matrix[,1] == 0)
509 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)) 509 F1.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))
510 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)) 510 F2.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))
511 511
512 N.newborn[row] <- num_insects.newborn 512 N.newborn[row] <- num_insects.newborn
513 N.death[row] <- num_insects.death 513 N.death[row] <- num_insects.death
514 N.adult[row] <- num_insects.adult 514 N.adult[row] <- num_insects.adult
515 } # end of days specified in the input temperature data 515 } # end of days specified in the input temperature data
525 S5.replications[,N.replications] <- S5 525 S5.replications[,N.replications] <- S5
526 newborn.replications[,N.replications] <- N.newborn 526 newborn.replications[,N.replications] <- N.newborn
527 death.replications[,N.replications] <- N.death 527 death.replications[,N.replications] <- N.death
528 adult.replications[,N.replications] <- N.adult 528 adult.replications[,N.replications] <- N.adult
529 pop.replications[,N.replications] <- tot.pop 529 pop.replications[,N.replications] <- tot.pop
530 g0.replications[,N.replications] <- gen0.pop 530 P.replications[,N.replications] <- gen0.pop
531 g1.replications[,N.replications] <- gen1.pop 531 F1.replications[,N.replications] <- gen1.pop
532 g2.replications[,N.replications] <- gen2.pop 532 F2.replications[,N.replications] <- gen2.pop
533 g0a.replications[,N.replications] <- g0.adult 533 P_adults.replications[,N.replications] <- P.adult
534 g1a.replications[,N.replications] <- g1.adult 534 F1_adults.replications[,N.replications] <- F1.adult
535 F2_adults.replications[,N.replications] <- g2.adult 535 F2_adults.replications[,N.replications] <- F2.adult
536 } 536 }
537 537
538 # Mean value for adults. 538 # Mean value for adults.
539 adults <- apply((S3.replications + S4.replications + S5.replications), 1, mean) 539 adults <- apply((S3.replications + S4.replications + S5.replications), 1, mean)
540 # Mean value for nymphs. 540 # Mean value for nymphs.
541 nymphs <- apply((S1.replications + S2.replications), 1, mean) 541 nymphs <- apply((S1.replications + S2.replications), 1, mean)
542 # Mean value for eggs. 542 # Mean value for eggs.
543 eggs <- apply(S0.replications, 1, mean) 543 eggs <- apply(S0.replications, 1, mean)
544 # Mean value for P. 544 # Mean value for P.
545 g0 <- apply(g0.replications, 1, mean) 545 P <- apply(P.replications, 1, mean)
546 # Mean value for F1. 546 # Mean value for F1.
547 g1 <- apply(g1.replications, 1, mean) 547 F1 <- apply(F1.replications, 1, mean)
548 # Mean value for F2. 548 # Mean value for F2.
549 g2 <- apply(g2.replications, 1, mean) 549 F2 <- apply(F2.replications, 1, mean)
550 # Mean value for P adults. 550 # Mean value for P adults.
551 g0a <- apply(g0a.replications, 1, mean) 551 P_adults <- apply(P_adults.replications, 1, mean)
552 # Mean value for F1 adults. 552 # Mean value for F1 adults.
553 g1a <- apply(g1a.replications, 1, mean) 553 F1_adults <- apply(F1_adults.replications, 1, mean)
554 # Mean value for F2 adults. 554 # Mean value for F2 adults.
555 F2_adults <- apply(F2_adults.replications, 1, mean) 555 F2_adults <- apply(F2_adults.replications, 1, mean)
556 556
557 # Standard error for adults. 557 # Standard error for adults.
558 adults.std_error <- apply((S3.replications + S4.replications + S5.replications), 1, sd) / sqrt(opt$replications) 558 adults.std_error <- apply((S3.replications + S4.replications + S5.replications), 1, sd) / sqrt(opt$replications)
559 # Standard error for nymphs. 559 # Standard error for nymphs.
560 nymphs.std_error <- apply((S1.replications + S2.replications) / sqrt(opt$replications), 1, sd) 560 nymphs.std_error <- apply((S1.replications + S2.replications) / sqrt(opt$replications), 1, sd)
561 # Standard error for eggs. 561 # Standard error for eggs.
562 eggs.std_error <- apply(S0.replications, 1, sd) / sqrt(opt$replications) 562 eggs.std_error <- apply(S0.replications, 1, sd) / sqrt(opt$replications)
563 # Standard error value for P. 563 # Standard error value for P.
564 g0.std_error <- apply(g0.replications, 1, sd) / sqrt(opt$replications) 564 P.std_error <- apply(P.replications, 1, sd) / sqrt(opt$replications)
565 # Standard error for F1. 565 # Standard error for F1.
566 g1.std_error <- apply(g1.replications, 1, sd) / sqrt(opt$replications) 566 F1.std_error <- apply(F1.replications, 1, sd) / sqrt(opt$replications)
567 # Standard error for F2. 567 # Standard error for F2.
568 g2.std_error <- apply(g2.replications, 1, sd) / sqrt(opt$replications) 568 F2.std_error <- apply(F2.replications, 1, sd) / sqrt(opt$replications)
569 # Standard error for P adult. 569 # Standard error for P adult.
570 g0a.std_error <- apply(g0a.replications, 1, sd) / sqrt(opt$replications) 570 P_adults.std_error <- apply(P_adults.replications, 1, sd) / sqrt(opt$replications)
571 # Standard error for F1 adult. 571 # Standard error for F1 adult.
572 g1a.std_error <- apply(g1a.replications, 1, sd) / sqrt(opt$replications) 572 F1_adults.std_error <- apply(F1_adults.replications, 1, sd) / sqrt(opt$replications)
573 # Standard error for F2 adult. 573 # Standard error for F2 adult.
574 F2_adults.std_error <- apply(F2_adults.replications, 1, sd) / sqrt(opt$replications) 574 F2_adults.std_error <- apply(F2_adults.replications, 1, sd) / sqrt(opt$replications)
575 575
576 dev.new(width=20, height=30) 576 dev.new(width=20, height=30)
577 577
609 lines (day.all, eggs-eggs.std_error, col=4, lty=2) 609 lines (day.all, eggs-eggs.std_error, col=4, lty=2)
610 } 610 }
611 611
612 # Subfigure 2: population size by generation 612 # Subfigure 2: population size by generation
613 title <- paste(opt$insect, ": Total pop. by generation :", opt$location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ") 613 title <- paste(opt$insect, ": Total pop. by generation :", opt$location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ")
614 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) 614 plot(day.all, P, main=title, type="l", ylim=c(0, max(F2)), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3)
615 lines(day.all, g1, lwd = 2, lty = 1, col=2) 615 lines(day.all, F1, lwd = 2, lty = 1, col=2)
616 lines(day.all, g2, lwd = 2, lty = 1, col=4) 616 lines(day.all, F2, lwd = 2, lty = 1, col=4)
617 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")) 617 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"))
618 axis(2, cex.axis=3) 618 axis(2, cex.axis=3)
619 legend("topleft", c("P", "F1", "F2"), lty=c(1, 1, 1), col=c(1, 2, 4), cex=3) 619 legend("topleft", c("P", "F1", "F2"), lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)
620 if (opt$std_error_plot == 1) { 620 if (opt$std_error_plot == 1) {
621 # Add Standard error lines to plot 621 # Add Standard error lines to plot
622 # Standard error for adults 622 # Standard error for adults
623 lines (day.all, g0+g0.std_error, lty=2) 623 lines (day.all, P+P.std_error, lty=2)
624 lines (day.all, g0-g0.std_error, lty=2) 624 lines (day.all, P-P.std_error, lty=2)
625 # Standard error for nymphs 625 # Standard error for nymphs
626 lines (day.all, g1+g1.std_error, col=2, lty=2) 626 lines (day.all, F1+F1.std_error, col=2, lty=2)
627 lines (day.all, g1-g1.std_error, col=2, lty=2) 627 lines (day.all, F1-F1.std_error, col=2, lty=2)
628 # Standard error for eggs 628 # Standard error for eggs
629 lines (day.all, g2+g2.std_error, col=4, lty=2) 629 lines (day.all, F2+F2.std_error, col=4, lty=2)
630 lines (day.all, g2-g2.std_error, col=4, lty=2) 630 lines (day.all, F2-F2.std_error, col=4, lty=2)
631 } 631 }
632 632
633 # Subfigure 3: adult population size by generation 633 # 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=" ") 634 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(F2_adults) + 100), main=title, type="l", axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3) 635 plot(day.all, P_adults, ylim=c(0, max(F2_adults) + 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) 636 lines(day.all, F1_adults, lwd = 2, lty = 1, col=2)
637 lines(day.all, F2_adults, lwd = 2, lty = 1, col=4) 637 lines(day.all, F2_adults, 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")) 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"))
639 axis(2, cex.axis=3) 639 axis(2, cex.axis=3)
640 legend("topleft", c("P", "F1", "F2"), lty=c(1, 1, 1), col=c(1, 2, 4), cex=3) 640 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) { 641 if (opt$std_error_plot == 1) {
642 # Add Standard error lines to plot 642 # Add Standard error lines to plot
643 # Standard error for adults 643 # Standard error for adults
644 lines (day.all, g0a+g0a.std_error, lty=2) 644 lines (day.all, P_adults+P_adults.std_error, lty=2)
645 lines (day.all, g0a-g0a.std_error, lty=2) 645 lines (day.all, P_adults-P_adults.std_error, lty=2)
646 # Standard error for nymphs 646 # Standard error for nymphs
647 lines (day.all, g1a+g1a.std_error, col=2, lty=2) 647 lines (day.all, F1_adults+F1_adults.std_error, col=2, lty=2)
648 lines (day.all, g1a-g1a.std_error, col=2, lty=2) 648 lines (day.all, F1_adults-F1_adults.std_error, col=2, lty=2)
649 # Standard error for eggs 649 # Standard error for eggs
650 lines (day.all, F2_adults+F2_adults.std_error, col=4, lty=2) 650 lines (day.all, F2_adults+F2_adults.std_error, col=4, lty=2)
651 lines (day.all, F2_adults-F2_adults.std_error, col=4, lty=2) 651 lines (day.all, F2_adults-F2_adults.std_error, col=4, lty=2)
652 } 652 }
653 653