Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 96:f0c65c9f1078 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 10:57:59 -0500
parents 89492de287a6
children 58bc1a2ca936
comparison
equal deleted inserted replaced
95:08571981adb3 96:f0c65c9f1078
109 } 109 }
110 else { 110 else {
111 dh[i] <- T[i] - 8.4 111 dh[i] <- T[i] - 8.4
112 } 112 }
113 } 113 }
114 else if (i > settime) { 114 else if (i > settime) {
115 n <- i - settime 115 n <- i - settime
116 T[i] = dnp + (ts - dnp) * exp( - b * n / z) 116 T[i] = dnp + (ts - dnp) * exp( - b * n / z)
117 if (T[i] < 8.4) { 117 if (T[i] < 8.4) {
118 dh[i] <- 0 118 dh[i] <- 0
119 } 119 }
199 # so get the value from the first row. 199 # so get the value from the first row.
200 latitude <- temperature_data_frame$LATITUDE[1] 200 latitude <- temperature_data_frame$LATITUDE[1]
201 201
202 cat("Number of days: ", opt$num_days, "\n") 202 cat("Number of days: ", opt$num_days, "\n")
203 203
204 # Initialize matrix for results from all replications. 204 # Initialize matrices for results from all replications.
205 S0.rep <- S1.rep <- S2.rep <- S3.rep <- S4.rep <- S5.rep <- matrix(rep(0, opt$num_days * opt$replications), ncol = opt$replications) 205 S0.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
206 newborn.rep <- death.rep <- adult.rep <- pop.rep <- g0.rep <- g1.rep <- g2.rep <- g0a.rep <- g1a.rep <- g2a.rep <- matrix(rep(0, opt$num_days * opt$replications), ncol=opt$replications) 206 S1.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
207 S2.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
208 S3.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
209 S4.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
210 S5.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
211 newborn.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
212 death.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
213 adult.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
214 pop.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
215 g0.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
216 g1.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
217 g2.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
218 g0a.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
219 g1a.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
220 g2a.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
207 221
208 # Loop through replications. 222 # Loop through replications.
209 for (N.rep in 1:opt$replications) { 223 for (N.rep in 1:opt$replications) {
210 # During each replication start with 1000 individuals. 224 # During each replication start with 1000 individuals.
211 # TODO: user definable as well? 225 # TODO: user definable as well?
272 } 286 }
273 # (or dependent on temperature and life stage?) 287 # (or dependent on temperature and life stage?)
274 u.d <- runif(1) 288 u.d <- runif(1)
275 if (u.d < death.prob) { 289 if (u.d < death.prob) {
276 death.vec <- c(death.vec, i) 290 death.vec <- c(death.vec, i)
277 } 291 }
278 else { 292 else {
279 # Aggregrate index of dead bug. 293 # Aggregrate index of dead bug.
280 # Event 1 end of diapause. 294 # Event 1 end of diapause.
281 if (vec.ind[1] == 0 && vec.ind[2] == 3) { 295 if (vec.ind[1] == 0 && vec.ind[2] == 3) {
282 # Overwintering adult (previttelogenic). 296 # Overwintering adult (previttelogenic).
410 else { 424 else {
411 # Add 1 day in current stage. 425 # Add 1 day in current stage.
412 vec.ind[4] <- vec.ind[4] + 1 426 vec.ind[4] <- vec.ind[4] + 1
413 } 427 }
414 vec.mat[i,] <- vec.ind 428 vec.mat[i,] <- vec.ind
415 } 429 }
416 # Event 3.3 old nymph to adult: previttelogenic or diapausing? 430 # Event 3.3 old nymph to adult: previttelogenic or diapausing?
417 if (vec.ind[2] == 2) { 431 if (vec.ind[2] == 2) {
418 # Old nymph stage. 432 # Old nymph stage.
419 # Add to degree_days. 433 # Add to degree_days.
420 vec.ind[3] <- vec.ind[3] + degree_days.temp 434 vec.ind[3] <- vec.ind[3] + degree_days.temp
452 # Remove record of dead. 466 # Remove record of dead.
453 # Find how many new born. 467 # Find how many new born.
454 num_insects.newborn <- length(birth.vec[,1]) 468 num_insects.newborn <- length(birth.vec[,1])
455 vec.mat <- rbind(vec.mat, birth.vec) 469 vec.mat <- rbind(vec.mat, birth.vec)
456 # Update population size for the next day. 470 # Update population size for the next day.
457 num_insects <- num_insects - num_insects.death + num_insects.newborn 471 num_insects <- num_insects - num_insects.death + num_insects.newborn
458 472
459 # Aggregate results by day. 473 # Aggregate results by day.
460 tot.pop <- c(tot.pop, num_insects) 474 tot.pop <- c(tot.pop, num_insects)
461 # Egg. 475 # Egg.
462 s0 <- sum(vec.mat[,2] == 0) 476 sum_eggs <- sum(vec.mat[,2] == 0)
463 # Young nymph. 477 # Young nymph.
464 s1 <- sum(vec.mat[,2] == 1) 478 sum_young_nymphs <- sum(vec.mat[,2] == 1)
465 # Old nymph. 479 # Old nymph.
466 s2 <- sum(vec.mat[,2] == 2) 480 sum_old_nymphs <- sum(vec.mat[,2] == 2)
467 # Previtellogenic. 481 # Previtellogenic.
468 s3 <- sum(vec.mat[,2] == 3) 482 sum_previtellogenic <- sum(vec.mat[,2] == 3)
469 # Vitellogenic. 483 # Vitellogenic.
470 s4 <- sum(vec.mat[,2] == 4) 484 sum_vitellogenic <- sum(vec.mat[,2] == 4)
471 # Diapausing. 485 # Diapausing.
472 s5 <- sum(vec.mat[,2] == 5) 486 sum_diapausing <- sum(vec.mat[,2] == 5)
473 # Overwintering adult. 487 # Overwintering adult.
474 gen0 <- sum(vec.mat[,1] == 0) 488 gen0 <- sum(vec.mat[,1] == 0)
475 # First generation. 489 # First generation.
476 gen1 <- sum(vec.mat[,1] == 1) 490 gen1 <- sum(vec.mat[,1] == 1)
477 # Second generation. 491 # Second generation.
482 # Generation 0 pop size. 496 # Generation 0 pop size.
483 gen0.pop[row] <- gen0 497 gen0.pop[row] <- gen0
484 gen1.pop[row] <- gen1 498 gen1.pop[row] <- gen1
485 gen2.pop[row] <- gen2 499 gen2.pop[row] <- gen2
486 500
487 S0[row] <- s0 501 S0[row] <- sum_eggs
488 S1[row] <- s1 502 S1[row] <- sum_young_nymphs
489 S2[row] <- s2 503 S2[row] <- sum_old_nymphs
490 S3[row] <- s3 504 S3[row] <- sum_previtellogenic
491 S4[row] <- s4 505 S4[row] <- sum_vitellogenic
492 S5[row] <- s5 506 S5[row] <- sum_diapausing
493 507
494 g0.adult[row] <- sum(vec.mat[,1] == 0) 508 g0.adult[row] <- sum(vec.mat[,1] == 0)
495 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)) 509 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))
496 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)) 510 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))
497 511
519 g0a.rep[,N.rep] <- g0.adult 533 g0a.rep[,N.rep] <- g0.adult
520 g1a.rep[,N.rep] <- g1.adult 534 g1a.rep[,N.rep] <- g1.adult
521 g2a.rep[,N.rep] <- g2.adult 535 g2a.rep[,N.rep] <- g2.adult
522 } 536 }
523 537
524 # Mean value for adults 538 # Mean value for adults.
525 mean_value_adult <- apply((S3.rep + S4.rep + S5.rep), 1, mean) 539 mean_value_adult <- apply((S3.rep + S4.rep + S5.rep), 1, mean)
526 # Mean value for nymphs 540 # Mean value for nymphs.
527 mean_value_nymph <- apply((S1.rep + S2.rep), 1, mean) 541 mean_value_nymph <- apply((S1.rep + S2.rep), 1, mean)
528 # Mean value for eggs 542 # Mean value for eggs.
529 mean_value_egg <- apply(S0.rep, 1, mean) 543 mean_value_egg <- apply(S0.rep, 1, mean)
530 # Mean value for P 544 # Mean value for P.
531 g0 <- apply(g0.rep, 1, mean) 545 g0 <- apply(g0.rep, 1, mean)
532 # Mean value for F1 546 # Mean value for F1.
533 g1 <- apply(g1.rep, 1, mean) 547 g1 <- apply(g1.rep, 1, mean)
534 # Mean value for F2 548 # Mean value for F2.
535 g2 <- apply(g2.rep, 1, mean) 549 g2 <- apply(g2.rep, 1, mean)
536 # Mean value for P adult 550 # Mean value for P adult.
537 g0a <- apply(g0a.rep, 1, mean) 551 g0a <- apply(g0a.rep, 1, mean)
538 # Mean value for F1 adult 552 # Mean value for F1 adult.
539 g1a <- apply(g1a.rep, 1, mean) 553 g1a <- apply(g1a.rep, 1, mean)
540 # Mean value for F2 adult 554 # Mean value for F2 adult.
541 g2a <- apply(g2a.rep, 1, mean) 555 g2a <- apply(g2a.rep, 1, mean)
542 556
543 # Standard error for adults 557 # Standard error for adults.
544 mean_value_adult.std_error <- apply((S3.rep + S4.rep + S5.rep), 1, sd) / sqrt(opt$replications) 558 mean_value_adult.std_error <- apply((S3.rep + S4.rep + S5.rep), 1, sd) / sqrt(opt$replications)
545 # Standard error for nymphs 559 # Standard error for nymphs.
546 mean_value_nymph.std_error <- apply((S1.rep + S2.rep) / sqrt(opt$replications), 1, sd) 560 mean_value_nymph.std_error <- apply((S1.rep + S2.rep) / sqrt(opt$replications), 1, sd)
547 # Standard error for eggs 561 # Standard error for eggs.
548 mean_value_egg.std_error <- apply(S0.rep, 1, sd) / sqrt(opt$replications) 562 mean_value_egg.std_error <- apply(S0.rep, 1, sd) / sqrt(opt$replications)
549 # Standard error value for P 563 # Standard error value for P.
550 g0.std_error <- apply(g0.rep, 1, sd) / sqrt(opt$replications) 564 g0.std_error <- apply(g0.rep, 1, sd) / sqrt(opt$replications)
551 # Standard error for F1 565 # Standard error for F1.
552 g1.std_error <- apply(g1.rep, 1, sd) / sqrt(opt$replications) 566 g1.std_error <- apply(g1.rep, 1, sd) / sqrt(opt$replications)
553 # Standard error for F2 567 # Standard error for F2.
554 g2.std_error <- apply(g2.rep, 1, sd) / sqrt(opt$replications) 568 g2.std_error <- apply(g2.rep, 1, sd) / sqrt(opt$replications)
555 # Standard error for P adult 569 # Standard error for P adult.
556 g0a.std_error <- apply(g0a.rep, 1, sd) / sqrt(opt$replications) 570 g0a.std_error <- apply(g0a.rep, 1, sd) / sqrt(opt$replications)
557 # Standard error for F1 adult 571 # Standard error for F1 adult.
558 g1a.std_error <- apply(g1a.rep, 1, sd) / sqrt(opt$replications) 572 g1a.std_error <- apply(g1a.rep, 1, sd) / sqrt(opt$replications)
559 # Standard error for F2 adult 573 # Standard error for F2 adult.
560 g2a.std_error <- apply(g2a.rep, 1, sd) / sqrt(opt$replications) 574 g2a.std_error <- apply(g2a.rep, 1, sd) / sqrt(opt$replications)
561 575
562 dev.new(width=20, height=30) 576 dev.new(width=20, height=30)
563 577
564 # Start PDF device driver to save charts to output. 578 # Start PDF device driver to save charts to output.
581 lines(day.all, mean_value_egg, lwd=2, lty=1, col=4) 595 lines(day.all, mean_value_egg, lwd=2, lty=1, col=4)
582 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"))
583 axis(2, cex.axis=3) 597 axis(2, cex.axis=3)
584 leg.text <- c("Egg", "Nymph", "Adult") 598 leg.text <- c("Egg", "Nymph", "Adult")
585 legend("topleft", leg.text, lty=c(1, 1, 1), col=c(4, 2, 1), cex=3) 599 legend("topleft", leg.text, lty=c(1, 1, 1), col=c(4, 2, 1), cex=3)
586 if (opt$se_plot == 1) { 600 if (opt$std_error_plot == 1) {
587 # Add Standard error lines to plot 601 # Add Standard error lines to plot
588 # Standard error for adults 602 # Standard error for adults
589 lines (day.all, mean_value_adult+mean_value_adult.std_error, lty=2) 603 lines (day.all, mean_value_adult+mean_value_adult.std_error, lty=2)
590 lines (day.all, mean_value_adult-mean_value_adult.std_error, lty=2) 604 lines (day.all, mean_value_adult-mean_value_adult.std_error, lty=2)
591 # Standard error for nymphs 605 # Standard error for nymphs
592 lines (day.all, mean_value_nymph+mean_value_nymph.std_error, col=2, lty=2) 606 lines (day.all, mean_value_nymph+mean_value_nymph.std_error, col=2, lty=2)
593 lines (day.all, mean_value_nymph-mean_value_nymph.std_error, col=2, lty=2) 607 lines (day.all, mean_value_nymph-mean_value_nymph.std_error, col=2, lty=2)
594 # Standard error for eggs 608 # Standard error for eggs
595 lines (day.all, mean_value_egg+mean_value_egg.std_error, col=4, lty=2) 609 lines (day.all, mean_value_egg+mean_value_egg.std_error, col=4, lty=2)
603 lines(day.all, g2, lwd = 2, lty = 1, col=4) 617 lines(day.all, g2, lwd = 2, lty = 1, col=4)
604 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(1, at=c(1:12) * 30 - 15, cex.axis=3, labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))
605 axis(2, cex.axis=3) 619 axis(2, cex.axis=3)
606 leg.text <- c("P", "F1", "F2") 620 leg.text <- c("P", "F1", "F2")
607 legend("topleft", leg.text, lty=c(1, 1, 1), col=c(1, 2, 4), cex=3) 621 legend("topleft", leg.text, lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)
608 if (opt$se_plot == 1) { 622 if (opt$std_error_plot == 1) {
609 # Add Standard error lines to plot 623 # Add Standard error lines to plot
610 # Standard error for adults 624 # Standard error for adults
611 lines (day.all, g0+g0.std_error, lty=2) 625 lines (day.all, g0+g0.std_error, lty=2)
612 lines (day.all, g0-g0.std_error, lty=2) 626 lines (day.all, g0-g0.std_error, lty=2)
613 # Standard error for nymphs 627 # Standard error for nymphs
614 lines (day.all, g1+g1.std_error, col=2, lty=2) 628 lines (day.all, g1+g1.std_error, col=2, lty=2)
615 lines (day.all, g1-g1.std_error, col=2, lty=2) 629 lines (day.all, g1-g1.std_error, col=2, lty=2)
616 # Standard error for eggs 630 # Standard error for eggs
617 lines (day.all, g2+g2.std_error, col=4, lty=2) 631 lines (day.all, g2+g2.std_error, col=4, lty=2)
629 legend("topleft", leg.text, lty=c(1, 1, 1), col=c(1, 2, 4), cex=3) 643 legend("topleft", leg.text, lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)
630 if (opt$std_error_plot == 1) { 644 if (opt$std_error_plot == 1) {
631 # Add Standard error lines to plot 645 # Add Standard error lines to plot
632 # Standard error for adults 646 # Standard error for adults
633 lines (day.all, g0a+g0a.std_error, lty=2) 647 lines (day.all, g0a+g0a.std_error, lty=2)
634 lines (day.all, g0a-g0a.std_error, lty=2) 648 lines (day.all, g0a-g0a.std_error, lty=2)
635 # Standard error for nymphs 649 # Standard error for nymphs
636 lines (day.all, g1a+g1a.std_error, col=2, lty=2) 650 lines (day.all, g1a+g1a.std_error, col=2, lty=2)
637 lines (day.all, g1a-g1a.std_error, col=2, lty=2) 651 lines (day.all, g1a-g1a.std_error, col=2, lty=2)
638 # Standard error for eggs 652 # Standard error for eggs
639 lines (day.all, g2a+g2a.std_error, col=4, lty=2) 653 lines (day.all, g2a+g2a.std_error, col=4, lty=2)