Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 99:51045fde125e 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:04:11 -0500
parents d33a401c0153
children 52114847afff
comparison
equal deleted inserted replaced
98:d33a401c0153 99:51045fde125e
215 g0.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)
216 g1.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)
217 g2.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)
218 g0a.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)
219 g1a.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)
220 g2a.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.
225 # TODO: user definable as well? 225 # TODO: user definable as well?
465 vector.ind[4] <- vector.ind[4] + 1 465 vector.ind[4] <- vector.ind[4] + 1
466 vector.matrix[i,] <- vector.ind 466 vector.matrix[i,] <- vector.ind
467 } 467 }
468 } # Else if it is still alive. 468 } # Else if it is still alive.
469 } # End of the individual bug loop. 469 } # End of the individual bug loop.
470 # Find how many died. 470 # Find the number of deaths.
471 num_insects.death <- length(death.vector) 471 num_insects.death <- length(death.vector)
472 if (num_insects.death > 0) { 472 if (num_insects.death > 0) {
473 # Remove record of dead.
473 vector.matrix <- vector.matrix[-death.vector, ] 474 vector.matrix <- vector.matrix[-death.vector, ]
474 } 475 }
475 # Remove record of dead. 476 # Find the number of births.
476 # Find how many new born.
477 num_insects.newborn <- length(birth.vector[,1]) 477 num_insects.newborn <- length(birth.vector[,1])
478 vector.matrix <- rbind(vector.matrix, birth.vector) 478 vector.matrix <- rbind(vector.matrix, birth.vector)
479 # Update population size for the next day. 479 # Update population size for the next day.
480 num_insects <- num_insects - num_insects.death + num_insects.newborn 480 num_insects <- num_insects - num_insects.death + num_insects.newborn
481 481
482 # Aggregate results by day. 482 # Aggregate results by day.
483 tot.pop <- c(tot.pop, num_insects) 483 tot.pop <- c(tot.pop, num_insects)
484 # Egg. 484
485 sum_eggs <- sum(vector.matrix[,2] == 0) 485 # All adults population size.
486 # Young nymph.
487 sum_young_nymphs <- sum(vector.matrix[,2] == 1)
488 # Old nymph.
489 sum_old_nymphs <- sum(vector.matrix[,2] == 2)
490 # Previtellogenic.
491 sum_previtellogenic <- sum(vector.matrix[,2] == 3)
492 # Vitellogenic.
493 sum_vitellogenic <- sum(vector.matrix[,2] == 4)
494 # Diapausing.
495 sum_diapausing <- sum(vector.matrix[,2] == 5)
496 # Overwintering adult.
497 sum_overwintering_adults <- sum(vector.matrix[,1] == 0)
498 # First generation.
499 sum_first_generation <- sum(vector.matrix[,1] == 1)
500 # Second generation.
501 sum_second_generation <- sum(vector.matrix[,1] == 2)
502 # Sum of all adults.
503 num_insects.adult <- sum(vector.matrix[,2] == 3) + sum(vector.matrix[,2] == 4) + sum(vector.matrix[,2] == 5) 486 num_insects.adult <- sum(vector.matrix[,2] == 3) + sum(vector.matrix[,2] == 4) + sum(vector.matrix[,2] == 5)
504 487
505 # Population sizes. 488 # Overwintering adult population size.
506 gen0.pop[row] <- sum_overwintering_adults 489 gen0.pop[row] <- sum(vector.matrix[,1] == 0)
507 gen1.pop[row] <- sum_first_generation 490 # First generation population size.
508 gen2.pop[row] <- sum_second_generation 491 gen1.pop[row] <- sum(vector.matrix[,1] == 1)
509 492 # Second generation population size.
510 S0[row] <- sum_eggs 493 gen2.pop[row] <- sum(vector.matrix[,1] == 2)
511 S1[row] <- sum_young_nymphs 494
512 S2[row] <- sum_old_nymphs 495 # Egg population size.
513 S3[row] <- sum_previtellogenic 496 S0[row] <- sum(vector.matrix[,2] == 0)
514 S4[row] <- sum_vitellogenic 497 # Young nymph population size.
515 S5[row] <- sum_diapausing 498 S1[row] <- sum(vector.matrix[,2] == 1)
499 # Old nymph population size.
500 S2[row] <- sum(vector.matrix[,2] == 2)
501 # Previtellogenic population size.
502 S3[row] <- sum(vector.matrix[,2] == 3)
503 # Vitellogenic population size.
504 S4[row] <- sum(vector.matrix[,2] == 4)
505 # Diapausing population size.
506 S5[row] <- sum(vector.matrix[,2] == 5)
516 507
517 g0.adult[row] <- sum(vector.matrix[,1] == 0) 508 g0.adult[row] <- sum(vector.matrix[,1] == 0)
518 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 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))
519 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 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))
520 511
539 g0.replications[,N.replications] <- gen0.pop 530 g0.replications[,N.replications] <- gen0.pop
540 g1.replications[,N.replications] <- gen1.pop 531 g1.replications[,N.replications] <- gen1.pop
541 g2.replications[,N.replications] <- gen2.pop 532 g2.replications[,N.replications] <- gen2.pop
542 g0a.replications[,N.replications] <- g0.adult 533 g0a.replications[,N.replications] <- g0.adult
543 g1a.replications[,N.replications] <- g1.adult 534 g1a.replications[,N.replications] <- g1.adult
544 g2a.replications[,N.replications] <- g2.adult 535 F2_adults.replications[,N.replications] <- g2.adult
545 } 536 }
546 537
547 # Mean value for adults. 538 # Mean value for adults.
548 adult_mean <- apply((S3.replications + S4.replications + S5.replications), 1, mean) 539 adults <- apply((S3.replications + S4.replications + S5.replications), 1, mean)
549 # Mean value for nymphs. 540 # Mean value for nymphs.
550 nymph_mean <- apply((S1.replications + S2.replications), 1, mean) 541 nymphs <- apply((S1.replications + S2.replications), 1, mean)
551 # Mean value for eggs. 542 # Mean value for eggs.
552 egg_mean <- apply(S0.replications, 1, mean) 543 eggs <- apply(S0.replications, 1, mean)
553 # Mean value for P. 544 # Mean value for P.
554 g0 <- apply(g0.replications, 1, mean) 545 g0 <- apply(g0.replications, 1, mean)
555 # Mean value for F1. 546 # Mean value for F1.
556 g1 <- apply(g1.replications, 1, mean) 547 g1 <- apply(g1.replications, 1, mean)
557 # Mean value for F2. 548 # Mean value for F2.
558 g2 <- apply(g2.replications, 1, mean) 549 g2 <- apply(g2.replications, 1, mean)
559 # Mean value for P adult. 550 # Mean value for P adults.
560 g0a <- apply(g0a.replications, 1, mean) 551 g0a <- apply(g0a.replications, 1, mean)
561 # Mean value for F1 adult. 552 # Mean value for F1 adults.
562 g1a <- apply(g1a.replications, 1, mean) 553 g1a <- apply(g1a.replications, 1, mean)
563 # Mean value for F2 adult. 554 # Mean value for F2 adults.
564 F2_adult_mean <- apply(g2a.replications, 1, mean) 555 F2_adults <- apply(F2_adults.replications, 1, mean)
565 556
566 # Standard error for adults. 557 # Standard error for adults.
567 adult_mean.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)
568 # Standard error for nymphs. 559 # Standard error for nymphs.
569 nymph_mean.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)
570 # Standard error for eggs. 561 # Standard error for eggs.
571 egg_mean.std_error <- apply(S0.replications, 1, sd) / sqrt(opt$replications) 562 eggs.std_error <- apply(S0.replications, 1, sd) / sqrt(opt$replications)
572 # Standard error value for P. 563 # Standard error value for P.
573 g0.std_error <- apply(g0.replications, 1, sd) / sqrt(opt$replications) 564 g0.std_error <- apply(g0.replications, 1, sd) / sqrt(opt$replications)
574 # Standard error for F1. 565 # Standard error for F1.
575 g1.std_error <- apply(g1.replications, 1, sd) / sqrt(opt$replications) 566 g1.std_error <- apply(g1.replications, 1, sd) / sqrt(opt$replications)
576 # Standard error for F2. 567 # Standard error for F2.
578 # Standard error for P adult. 569 # Standard error for P adult.
579 g0a.std_error <- apply(g0a.replications, 1, sd) / sqrt(opt$replications) 570 g0a.std_error <- apply(g0a.replications, 1, sd) / sqrt(opt$replications)
580 # Standard error for F1 adult. 571 # Standard error for F1 adult.
581 g1a.std_error <- apply(g1a.replications, 1, sd) / sqrt(opt$replications) 572 g1a.std_error <- apply(g1a.replications, 1, sd) / sqrt(opt$replications)
582 # Standard error for F2 adult. 573 # Standard error for F2 adult.
583 g2a.std_error <- apply(g2a.replications, 1, sd) / sqrt(opt$replications) 574 F2_adults.std_error <- apply(F2_adults.replications, 1, sd) / sqrt(opt$replications)
584 575
585 dev.new(width=20, height=30) 576 dev.new(width=20, height=30)
586 577
587 # Start PDF device driver to save charts to output. 578 # Start PDF device driver to save charts to output.
588 pdf(file=opt$output, width=20, height=30, bg="white") 579 pdf(file=opt$output, width=20, height=30, bg="white")
595 start_date <- temperature_data_frame$DATE[1] 586 start_date <- temperature_data_frame$DATE[1]
596 end_date <- temperature_data_frame$DATE[opt$num_days] 587 end_date <- temperature_data_frame$DATE[opt$num_days]
597 588
598 # Subfigure 1: population size by life stage 589 # Subfigure 1: population size by life stage
599 title <- paste(opt$insect, ": Total pop. by life stage :", opt$location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ") 590 title <- paste(opt$insect, ": Total pop. by life stage :", opt$location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ")
600 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) 591 plot(day.all, adults, main=title, type="l", ylim=c(0, max(eggs + eggs.std_error, nymphs + nymphs.std_error, adults + adults.std_error)), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3)
601 # Young and old nymphs. 592 # Young and old nymphs.
602 lines(day.all, nymph_mean, lwd=2, lty=1, col=2) 593 lines(day.all, nymphs, lwd=2, lty=1, col=2)
603 # Eggs 594 # Eggs
604 lines(day.all, egg_mean, lwd=2, lty=1, col=4) 595 lines(day.all, eggs, lwd=2, lty=1, col=4)
605 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")) 596 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"))
606 axis(2, cex.axis=3) 597 axis(2, cex.axis=3)
607 legend("topleft", c("Egg", "Nymph", "Adult"), lty=c(1, 1, 1), col=c(4, 2, 1), cex=3) 598 legend("topleft", c("Egg", "Nymph", "Adult"), lty=c(1, 1, 1), col=c(4, 2, 1), cex=3)
608 if (opt$std_error_plot == 1) { 599 if (opt$std_error_plot == 1) {
609 # Add Standard error lines to plot 600 # Add Standard error lines to plot
610 # Standard error for adults 601 # Standard error for adults
611 lines (day.all, adult_mean+adult_mean.std_error, lty=2) 602 lines (day.all, adults+adults.std_error, lty=2)
612 lines (day.all, adult_mean-adult_mean.std_error, lty=2) 603 lines (day.all, adults-adults.std_error, lty=2)
613 # Standard error for nymphs 604 # Standard error for nymphs
614 lines (day.all, nymph_mean+nymph_mean.std_error, col=2, lty=2) 605 lines (day.all, nymphs+nymphs.std_error, col=2, lty=2)
615 lines (day.all, nymph_mean-nymph_mean.std_error, col=2, lty=2) 606 lines (day.all, nymphs-nymphs.std_error, col=2, lty=2)
616 # Standard error for eggs 607 # Standard error for eggs
617 lines (day.all, egg_mean+egg_mean.std_error, col=4, lty=2) 608 lines (day.all, eggs+eggs.std_error, col=4, lty=2)
618 lines (day.all, egg_mean-egg_mean.std_error, col=4, lty=2) 609 lines (day.all, eggs-eggs.std_error, col=4, lty=2)
619 } 610 }
620 611
621 # Subfigure 2: population size by generation 612 # Subfigure 2: population size by generation
622 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=" ")
623 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, 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)
639 lines (day.all, g2-g2.std_error, col=4, lty=2) 630 lines (day.all, g2-g2.std_error, col=4, lty=2)
640 } 631 }
641 632
642 # Subfigure 3: adult population size by generation 633 # Subfigure 3: adult population size by generation
643 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=" ")
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) 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)
645 lines(day.all, g1a, lwd = 2, lty = 1, col=2) 636 lines(day.all, g1a, lwd = 2, lty = 1, col=2)
646 lines(day.all, F2_adult_mean, lwd = 2, lty = 1, col=4) 637 lines(day.all, F2_adults, lwd = 2, lty = 1, col=4)
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")) 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"))
648 axis(2, cex.axis=3) 639 axis(2, cex.axis=3)
649 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)
650 if (opt$std_error_plot == 1) { 641 if (opt$std_error_plot == 1) {
651 # Add Standard error lines to plot 642 # Add Standard error lines to plot
654 lines (day.all, g0a-g0a.std_error, lty=2) 645 lines (day.all, g0a-g0a.std_error, lty=2)
655 # Standard error for nymphs 646 # Standard error for nymphs
656 lines (day.all, g1a+g1a.std_error, col=2, lty=2) 647 lines (day.all, g1a+g1a.std_error, col=2, lty=2)
657 lines (day.all, g1a-g1a.std_error, col=2, lty=2) 648 lines (day.all, g1a-g1a.std_error, col=2, lty=2)
658 # Standard error for eggs 649 # Standard error for eggs
659 lines (day.all, F2_adult_mean+g2a.std_error, col=4, lty=2) 650 lines (day.all, F2_adults+F2_adults.std_error, col=4, lty=2)
660 lines (day.all, F2_adult_mean-g2a.std_error, col=4, lty=2) 651 lines (day.all, F2_adults-F2_adults.std_error, col=4, lty=2)
661 } 652 }
662 653
663 # Turn off device driver to flush output. 654 # Turn off device driver to flush output.
664 dev.off() 655 dev.off()