Mercurial > repos > greg > insect_phenology_model
changeset 16:309954bbe999 draft
Uploaded
| author | greg |
|---|---|
| date | Tue, 06 Mar 2018 08:26:42 -0500 |
| parents | dd86ee185113 |
| children | 6f31ade9a0f4 |
| files | insect_phenology_model.R |
| diffstat | 1 files changed, 141 insertions(+), 96 deletions(-) [+] |
line wrap: on
line diff
--- a/insect_phenology_model.R Mon Mar 05 08:27:27 2018 -0500 +++ b/insect_phenology_model.R Tue Mar 06 08:26:42 2018 -0500 @@ -11,7 +11,7 @@ make_option(c("--insects_per_replication"), action="store", dest="insects_per_replication", type="integer", help="Number of insects with which to start each replication"), make_option(c("--life_stages"), action="store", dest="life_stages", help="Selected life stages for plotting"), make_option(c("--life_stages_adult"), action="store", dest="life_stages_adult", default=NULL, help="Adult life stages for plotting"), - make_option(c("--life_stage_nymph"), action="store", dest="life_stage_nymph", default=NULL, help="Nymph life stages for plotting"), + make_option(c("--life_stages_nymph"), action="store", dest="life_stages_nymph", default=NULL, help="Nymph life stages for plotting"), make_option(c("--location"), action="store", dest="location", help="Selected location"), make_option(c("--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"), make_option(c("--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"), @@ -219,7 +219,7 @@ render_chart = function(date_labels, chart_type, plot_std_error, insect, location, latitude, start_date, end_date, days, maxval, replications, life_stage, group, group_std_error, group2=NULL, group2_std_error=NULL, group3=NULL, group3_std_error=NULL, - life_stages_adult=NULL, life_stage_nymph=NULL) { + life_stages_adult=NULL, life_stages_nymph=NULL) { if (chart_type=="pop_size_by_life_stage") { if (life_stage=="Total") { title = paste(insect, ": Reps", replications, ":", life_stage, "Pop :", location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" "); @@ -248,9 +248,9 @@ legend_text = c(life_stage); columns = c(4); } else if (life_stage=="Nymph") { - stage = paste(life_stage_nymph, "Nymph Pop :", sep=" "); + stage = paste(life_stages_nymph, "Nymph Pop :", sep=" "); title = paste(insect, ": Reps", replications, ":", stage, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" "); - legend_text = c(paste(life_stage_nymph, life_stage, sep=" ")); + legend_text = c(paste(life_stages_nymph, life_stage, sep=" ")); columns = c(2); } else if (life_stage=="Adult") { stage = paste(life_stages_adult, "Adult Pop", sep=" "); @@ -274,7 +274,7 @@ } else if (life_stage=="Egg") { title_str = ": Egg Pop by Gen :"; } else if (life_stage=="Nymph") { - title_str = paste(":", life_stage_nymph, "Nymph Pop by Gen", ":", sep=" "); + title_str = paste(":", life_stages_nymph, "Nymph Pop by Gen", ":", sep=" "); } else if (life_stage=="Adult") { title_str = paste(":", life_stages_adult, "Adult Pop by Gen", ":", sep=" "); } @@ -327,20 +327,25 @@ if (life_stage=="Total") { process_eggs = TRUE; process_nymphs = TRUE; - life_stage_nymph = "Total"; process_adults = TRUE; - life_stages_adult = "Total"; } else if (life_stage=="Egg") { process_eggs = TRUE; } else if (life_stage=="Nymph") { process_nymphs = TRUE; - life_stage_nymph = opt$life_stage_nymph; } else if (life_stage=="Adult") { process_adults = TRUE; - life_stages_adult = opt$life_stages_adult; } } - +if (process_adults) { + # Split life_stages_adult into a list of strings for plots. + life_stages_adult_str = as.character(opt$life_stages_adult); + life_stages_adult = strsplit(life_stages_adult_str, ",")[[1]]; +} +if (process_nymphs) { +# Split life_stages_nymph into a list of strings for plots. + life_stages_nymph_str = as.character(opt$life_stages_nymph); + life_stages_nymph = strsplit(life_stages_nymph_str, ",")[[1]]; +} # Initialize matrices. if (process_eggs) { Eggs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); @@ -808,45 +813,49 @@ } if (process_nymphs) { # Calculate nymph populations for selected life stage. - if (life_stage_nymph=="Total") { - # Mean value for all nymphs. - nymphs = apply((YoungNymphs.replications+OldNymphs.replications), 1, mean); - # Standard error for all nymphs. - nymphs.std_error = apply((YoungNymphs.replications+OldNymphs.replications) / sqrt(opt$replications), 1, sd); - } else if (life_stage_nymph=="Young") { - # Mean value for young nymphs. - nymphs = apply(YoungNymphs.replications, 1, mean); - # Standard error for young nymphs. - nymphs.std_error = apply(YoungNymphs.replications / sqrt(opt$replications), 1, sd); - } else if (life_stage_nymph=="Old") { - # Mean value for old nymphs. - nymphs = apply(OldNymphs.replications, 1, mean); - # Standard error for old nymphs. - nymphs.std_error = apply(OldNymphs.replications / sqrt(opt$replications), 1, sd); + for (life_stage_nymph in life_stages_nymph) { + if (life_stages_nymph=="Total") { + # Mean value for all nymphs. + total_nymphs = apply((YoungNymphs.replications+OldNymphs.replications), 1, mean); + # Standard error for all nymphs. + total_nymphs.std_error = apply((YoungNymphs.replications+OldNymphs.replications) / sqrt(opt$replications), 1, sd); + } else if (life_stages_nymph=="Young") { + # Mean value for young nymphs. + young_nymphs = apply(YoungNymphs.replications, 1, mean); + # Standard error for young nymphs. + young_nymphs.std_error = apply(YoungNymphs.replications / sqrt(opt$replications), 1, sd); + } else if (life_stages_nymph=="Old") { + # Mean value for old nymphs. + old_nymphs = apply(OldNymphs.replications, 1, mean); + # Standard error for old nymphs. + old_nymphs.std_error = apply(OldNymphs.replications / sqrt(opt$replications), 1, sd); + } } } if (process_adults) { # Calculate adult populations for selected life stage. - if (life_stages_adult=="Total") { - # Mean value for all adults. - adults = apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, mean); - # Standard error for all adults. - adults.std_error = apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, sd) / sqrt(opt$replications); - } else if (life_stages_adult == "Pre-vittelogenic") { - # Mean value for previtellogenic adults. - adults = apply(Previtellogenic.replications, 1, mean); - # Standard error for previtellogenic adults. - adults.std_error = apply(Previtellogenic.replications, 1, sd) / sqrt(opt$replications); - } else if (life_stages_adult == "Vittelogenic") { - # Mean value for vitellogenic adults. - adults = apply(Vitellogenic.replications, 1, mean); - # Standard error for vitellogenic adults. - adults.std_error = apply(Vitellogenic.replications, 1, sd) / sqrt(opt$replications); - } else if (life_stages_adult == "Diapausing") { - # Mean value for vitellogenic adults. - adults = apply(Diapausing.replications, 1, mean); - # Standard error for vitellogenic adults. - adults.std_error = apply(Diapausing.replications, 1, sd) / sqrt(opt$replications); + for (life_stage_adult in life_stages_adult) { + if (life_stages_adult=="Total") { + # Mean value for all adults. + total_adults = apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, mean); + # Standard error for all adults. + total_adults.std_error = apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, sd) / sqrt(opt$replications); + } else if (life_stages_adult == "Pre-vittelogenic") { + # Mean value for previtellogenic adults. + previttelogenic_adults = apply(Previtellogenic.replications, 1, mean); + # Standard error for previtellogenic adults. + previttelogenic_adults.std_error = apply(Previtellogenic.replications, 1, sd) / sqrt(opt$replications); + } else if (life_stages_adult == "Vittelogenic") { + # Mean value for vitellogenic adults. + vittelogenic_adults = apply(Vitellogenic.replications, 1, mean); + # Standard error for vitellogenic adults. + vittelogenic_adults.std_error = apply(Vitellogenic.replications, 1, sd) / sqrt(opt$replications); + } else if (life_stages_adult == "Diapausing") { + # Mean value for vitellogenic adults. + diapausing_adults = apply(Diapausing.replications, 1, mean); + # Standard error for vitellogenic adults. + diapausing_adults.std_error = apply(Diapausing.replications, 1, sd) / sqrt(opt$replications); + } } } @@ -931,33 +940,37 @@ # Turn off device driver to flush output. dev.off(); } else if (life_stage == "Nymph") { - # Start PDF device driver. - dev.new(width=20, height=30); - file_name = paste(tolower(life_stage_nymph), "nymph_pop_by_generation.pdf", sep="_"); - file_path = paste(output_dir, file_name, sep="/"); - pdf(file=file_path, width=20, height=30, bg="white"); - par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); - # Nymph population size by generation. - maxval = max(F2_nymphs) + 100; - render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, - opt$replications, life_stage, group=P_nymphs, group_std_error=P_nymphs.std_error, group2=F1_nymphs, group2_std_error=F1_nymphs.std_error, - group3=F2_nymphs, group3_std_error=F2_nymphs.std_error, life_stage_nymph=life_stage_nymph); - # Turn off device driver to flush output. - dev.off(); + for (life_stage_nymph in life_stages_nymph) { + # Start PDF device driver. + dev.new(width=20, height=30); + file_name = paste(tolower(life_stage_nymph), "nymph_pop_by_generation.pdf", sep="_"); + file_path = paste(output_dir, file_name, sep="/"); + pdf(file=file_path, width=20, height=30, bg="white"); + par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); + # Nymph population size by generation. + maxval = max(F2_nymphs) + 100; + render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, + opt$replications, life_stage, group=P_nymphs, group_std_error=P_nymphs.std_error, group2=F1_nymphs, group2_std_error=F1_nymphs.std_error, + group3=F2_nymphs, group3_std_error=F2_nymphs.std_error, life_stages_nymph=life_stage_nymph); + # Turn off device driver to flush output. + dev.off(); + } } else if (life_stage == "Adult") { - # Start PDF device driver. - dev.new(width=20, height=30); - file_name = paste(tolower(life_stages_adult), "adult_pop_by_generation.pdf", sep="_"); - file_path = paste(output_dir, file_name, sep="/"); - pdf(file=file_path, width=20, height=30, bg="white"); - par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); - # Adult population size by generation. - maxval = max(F2_adults) + 100; - render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, - opt$replications, life_stage, group=P_adults, group_std_error=P_adults.std_error, group2=F1_adults, group2_std_error=F1_adults.std_error, - group3=F2_adults, group3_std_error=F2_adults.std_error, life_stages_adult=life_stages_adult); - # Turn off device driver to flush output. - dev.off(); + for (life_stage_adult in life_stages_adult) { + # Start PDF device driver. + dev.new(width=20, height=30); + file_name = paste(tolower(life_stage_adult), "adult_pop_by_generation.pdf", sep="_"); + file_path = paste(output_dir, file_name, sep="/"); + pdf(file=file_path, width=20, height=30, bg="white"); + par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); + # Adult population size by generation. + maxval = max(F2_adults) + 100; + render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, + opt$replications, life_stage, group=P_adults, group_std_error=P_adults.std_error, group2=F1_adults, group2_std_error=F1_adults.std_error, + group3=F2_adults, group3_std_error=F2_adults.std_error, life_stages_adult=life_stage_adult); + # Turn off device driver to flush output. + dev.off(); + } } else if (life_stage == "Total") { # Start PDF device driver. dev.new(width=20, height=30); @@ -987,31 +1000,63 @@ # Turn off device driver to flush output. dev.off(); } else if (life_stage == "Nymph") { - # Start PDF device driver. - dev.new(width=20, height=30); - file_name = paste(tolower(life_stage_nymph), "nymph_pop.pdf", sep="_"); - file_path = paste(output_dir, file_name, sep="/"); - pdf(file=file_path, width=20, height=30, bg="white"); - par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); - # Nymph population size. - maxval = max(nymphs+nymphs.std_error); - render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, - opt$replications, life_stage, group=nymphs, group_std_error=nymphs.std_error, life_stage_nymph=life_stage_nymph); - # Turn off device driver to flush output. - dev.off(); + for (life_stage_nymph in life_stages_nymph) { + # Start PDF device driver. + dev.new(width=20, height=30); + file_name = paste(tolower(life_stage_nymph), "nymph_pop.pdf", sep="_"); + file_path = paste(output_dir, file_name, sep="/"); + pdf(file=file_path, width=20, height=30, bg="white"); + par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); + if (life_stage_nymph=="Total") { + # Total nymph population size. + group = total_nymphs; + group_std_error = total_nymphs.std_error; + } else if (life_stage_nymph=="Young") { + # Young nymph population size. + group = young_nymphs; + group_std_error = young_nymphs.std_error; + } else if (life_stage_nymph=="Old") { + # Old nymph population size. + group = old_nymphs; + group_std_error = old_nymphs.std_error; + } + maxval = max(group+group.std_error); + render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, + opt$replications, life_stage, group=group, group_std_error=group_std_error, life_stages_nymph=life_stage_nymph); + # Turn off device driver to flush output. + dev.off(); + } } else if (life_stage == "Adult") { - # Start PDF device driver. - dev.new(width=20, height=30); - file_name = paste(tolower(life_stages_adult), "adult_pop.pdf", sep="_"); - file_path = paste(output_dir, file_name, sep="/"); - pdf(file=file_path, width=20, height=30, bg="white"); - par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); - # Adult population size. - maxval = max(adults+adults.std_error); - render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, - opt$replications, life_stage, group=adults, group_std_error=adults.std_error, life_stages_adult=life_stages_adult); - # Turn off device driver to flush output. - dev.off(); + for (life_stage_adult in life_stages_adult) { + # Start PDF device driver. + dev.new(width=20, height=30); + file_name = paste(tolower(life_stages_adult), "adult_pop.pdf", sep="_"); + file_path = paste(output_dir, file_name, sep="/"); + pdf(file=file_path, width=20, height=30, bg="white"); + par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); + if (life_stage_adult=="Total") { + # Total adult population size. + group = total_adults; + group_std_error = total_adults.std_error + } else if (life_stage_adult=="Pre-vittelogenic") { + # Pre-vittelogenic adult population size. + group = previttelogenic_adults; + group_std_error = previttelogenic_adults.std_error + } else if (life_stage_adult=="Vittelogenic") { + # Vittelogenic adult population size. + group = vittelogenic_adults; + group_std_error = vittelogenic_adults.std_error + } else if (life_stage_adult=="Diapausing") { + # Diapausing adult population size. + group = diapausing_adults; + group_std_error = diapausing_adults.std_error + } + maxval = max(group+group_std_error); + render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, + opt$replications, life_stage, group=group, group_std_error=group_std_error, life_stages_adult=life_stage_adult); + # Turn off device driver to flush output. + dev.off(); + } } else if (life_stage == "Total") { # Start PDF device driver. dev.new(width=20, height=30); @@ -1021,7 +1066,7 @@ # Total population size. maxval = max(eggs+eggs.std_error, nymphs+nymphs.std_error, adults+adults.std_error); render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, - opt$replications, life_stage, group=adults, group_std_error=adults.std_error, group2=nymphs, group2_std_error=nymphs.std_error, group3=eggs, + opt$replications, life_stage, group=total_adults, group_std_error=total_adults.std_error, group2=total_nymphs, group2_std_error=total_nymphs.std_error, group3=eggs, group3_std_error=eggs.std_error); # Turn off device driver to flush output. dev.off();