Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 27:452e0e189e84 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Fri, 09 Mar 2018 13:41:38 -0500
parents 390d89243074
children afe6d8bac0c0
comparison
equal deleted inserted replaced
26:93209c167a61 27:452e0e189e84
13 make_option(c("--life_stages_adult"), action="store", dest="life_stages_adult", default=NULL, help="Adult life stages for plotting"), 13 make_option(c("--life_stages_adult"), action="store", dest="life_stages_adult", default=NULL, help="Adult life stages for plotting"),
14 make_option(c("--life_stages_nymph"), action="store", dest="life_stages_nymph", default=NULL, help="Nymph life stages for plotting"), 14 make_option(c("--life_stages_nymph"), action="store", dest="life_stages_nymph", default=NULL, help="Nymph life stages for plotting"),
15 make_option(c("--location"), action="store", dest="location", help="Selected location"), 15 make_option(c("--location"), action="store", dest="location", help="Selected location"),
16 make_option(c("--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"), 16 make_option(c("--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"),
17 make_option(c("--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"), 17 make_option(c("--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"),
18 make_option(c("--num_days"), action="store", dest="num_days", type="integer", help="Total number of days in the temperature dataset"),
18 make_option(c("--nymph_mortality"), action="store", dest="nymph_mortality", type="integer", help="Adjustment rate for nymph mortality"), 19 make_option(c("--nymph_mortality"), action="store", dest="nymph_mortality", type="integer", help="Adjustment rate for nymph mortality"),
19 make_option(c("--old_nymph_accumulation"), action="store", dest="old_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (young nymph->old nymph)"), 20 make_option(c("--old_nymph_accumulation"), action="store", dest="old_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (young nymph->old nymph)"),
20 make_option(c("--num_days"), action="store", dest="num_days", type="integer", help="Total number of days in the temperature dataset"), 21 make_option(c("--output"), action="store", dest="output", help="Dataset containing analyzed data"),
21 make_option(c("--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"), 22 make_option(c("--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"),
22 make_option(c("--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"), 23 make_option(c("--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"),
23 make_option(c("--replications"), action="store", dest="replications", type="integer", help="Number of replications"),
24 make_option(c("--plot_generations_separately"), action="store", dest="plot_generations_separately", help="Plot Plot P, F1 and F2 as separate lines or pool across them"), 24 make_option(c("--plot_generations_separately"), action="store", dest="plot_generations_separately", help="Plot Plot P, F1 and F2 as separate lines or pool across them"),
25 make_option(c("--plot_std_error"), action="store", dest="plot_std_error", help="Plot Standard error"), 25 make_option(c("--plot_std_error"), action="store", dest="plot_std_error", help="Plot Standard error"),
26 make_option(c("--replications"), action="store", dest="replications", type="integer", help="Number of replications"),
26 make_option(c("--young_nymph_accumulation"), action="store", dest="young_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (egg->young nymph)") 27 make_option(c("--young_nymph_accumulation"), action="store", dest="young_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (egg->young nymph)")
27 ) 28 )
28 29
29 parser <- OptionParser(usage="%prog [options] file", option_list=option_list); 30 parser <- OptionParser(usage="%prog [options] file", option_list=option_list);
30 args <- parse_args(parser, positional_arguments=TRUE); 31 args <- parse_args(parser, positional_arguments=TRUE);
31 opt <- args$options; 32 opt <- args$options;
32 33
33 add_daylight_length = function(temperature_data_frame, num_columns, num_rows) { 34 add_daylight_length = function(temperature_data_frame, num_rows) {
34 # Return a vector of daylight length (photoperido profile) for 35 # Return a vector of daylight length (photoperido profile) for
35 # the number of days specified in the input temperature data 36 # the number of days specified in the input temperature data
36 # (from Forsythe 1995). 37 # (from Forsythe 1995).
37 p = 0.8333; 38 p = 0.8333;
38 latitude = temperature_data_frame$LATITUDE[1]; 39 latitude = temperature_data_frame$LATITUDE[1];
46 # Compute the length of daylight for the day of the year. 47 # Compute the length of daylight for the day of the year.
47 darkness_length = 24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi))); 48 darkness_length = 24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi)));
48 daylight_length_vector[i] = 24 - darkness_length; 49 daylight_length_vector[i] = 24 - darkness_length;
49 } 50 }
50 # Append daylight_length_vector as a new column to temperature_data_frame. 51 # Append daylight_length_vector as a new column to temperature_data_frame.
51 temperature_data_frame[, num_columns+1] = daylight_length_vector; 52 temperature_data_frame = append_vector(temperature_data_frame, daylight_length_vector, "DAYLEN");
52 return(temperature_data_frame); 53 return(temperature_data_frame);
54 }
55
56 append_vector = function(data_frame, vec, new_column_name) {
57 num_columns = dim(data_frame)[2];
58 current_column_names = colnames(data_frame);
59 # Append vector vec as a new column to data_frame.
60 data_frame[,num_columns+1] = vec;
61 # Reset the column names with the additional column for later access.
62 colnames(data_frame) = append(current_column_names, new_column_name);
63 return(data_frame);
53 } 64 }
54 65
55 get_date_labels = function(temperature_data_frame, num_rows) { 66 get_date_labels = function(temperature_data_frame, num_rows) {
56 # Keep track of the years to see if spanning years. 67 # Keep track of the years to see if spanning years.
57 month_labels = list(); 68 month_labels = list();
242 if (num_columns == 6) { 253 if (num_columns == 6) {
243 # The input data has the following 6 columns: 254 # The input data has the following 6 columns:
244 # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX 255 # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX
245 # Set the column names for access when adding daylight length.. 256 # Set the column names for access when adding daylight length..
246 colnames(temperature_data_frame) = c("LATITUDE","LONGITUDE", "DATE", "DOY", "TMIN", "TMAX"); 257 colnames(temperature_data_frame) = c("LATITUDE","LONGITUDE", "DATE", "DOY", "TMIN", "TMAX");
258 current_column_names = colnames(temperature_data_frame);
247 # Add a column containing the daylight length for each day. 259 # Add a column containing the daylight length for each day.
248 temperature_data_frame = add_daylight_length(temperature_data_frame, num_columns, num_rows); 260 temperature_data_frame = add_daylight_length(temperature_data_frame, num_rows);
249 # Reset the column names with the additional column for later access.
250 colnames(temperature_data_frame) = c("LATITUDE","LONGITUDE", "DATE", "DOY", "TMIN", "TMAX", "DAYLEN");
251 } 261 }
252 return(temperature_data_frame); 262 return(temperature_data_frame);
253 } 263 }
254 264
255 render_chart = function(date_labels, chart_type, plot_std_error, insect, location, latitude, start_date, end_date, days, maxval, 265 render_chart = function(date_labels, chart_type, plot_std_error, insect, location, latitude, start_date, end_date, days, maxval,
346 temperature_data_frame = parse_input_data(opt$input, opt$num_days); 356 temperature_data_frame = parse_input_data(opt$input, opt$num_days);
347 # Get the date labels for plots. 357 # Get the date labels for plots.
348 date_labels = get_date_labels(temperature_data_frame, opt$num_days); 358 date_labels = get_date_labels(temperature_data_frame, opt$num_days);
349 # All latitude values are the same, so get the value for plots from the first row. 359 # All latitude values are the same, so get the value for plots from the first row.
350 latitude = temperature_data_frame$LATITUDE[1]; 360 latitude = temperature_data_frame$LATITUDE[1];
351 # Get the number of days for plots.
352 num_columns = dim(temperature_data_frame)[2];
353 # Determine the specified life stages for processing. 361 # Determine the specified life stages for processing.
354 # Split life_stages into a list of strings for plots. 362 # Split life_stages into a list of strings for plots.
355 life_stages_str = as.character(opt$life_stages); 363 life_stages_str = as.character(opt$life_stages);
356 life_stages = strsplit(life_stages_str, ",")[[1]]; 364 life_stages = strsplit(life_stages_str, ",")[[1]];
357 # Determine the data we need to generate for plotting. 365 # Determine the data we need to generate for plotting.
617 death.vector = c(death.vector, i); 625 death.vector = c(death.vector, i);
618 } 626 }
619 else { 627 else {
620 # End of diapause. 628 # End of diapause.
621 if (vector.individual[1] == 0 && vector.individual[2] == 3) { 629 if (vector.individual[1] == 0 && vector.individual[2] == 3) {
622 # Overwintering adult (previttelogenic). 630 # Overwintering adult (pre-vittelogenic).
623 if (photoperiod > opt$photoperiod && vector.individual[3] > 68 && doy < 180) { 631 if (photoperiod > opt$photoperiod && vector.individual[3] > 68 && doy < 180) {
624 # Add 68C to become fully reproductively matured. 632 # Add 68C to become fully reproductively matured.
625 # Transfer to vittelogenic. 633 # Transfer to vittelogenic.
626 vector.individual = c(0, 4, 0, 0, 0); 634 vector.individual = c(0, 4, 0, 0, 0);
627 vector.matrix[i,] = vector.individual; 635 vector.matrix[i,] = vector.individual;
628 } 636 }
629 else { 637 else {
630 # Add to # Add average temperature for current day. 638 # Add average temperature for current day.
631 vector.individual[3] = vector.individual[3] + averages.temp; 639 vector.individual[3] = vector.individual[3] + averages.temp;
632 # Add 1 day in current stage. 640 # Add 1 day in current stage.
633 vector.individual[4] = vector.individual[4] + 1; 641 vector.individual[4] = vector.individual[4] + 1;
634 vector.matrix[i,] = vector.individual; 642 vector.matrix[i,] = vector.individual;
635 } 643 }
636 } 644 }
637 if (vector.individual[1] != 0 && vector.individual[2] == 3) { 645 if (vector.individual[1] != 0 && vector.individual[2] == 3) {
638 # Not overwintering adult (previttelogenic). 646 # Not overwintering adult (pre-vittelogenic).
639 current.gen = vector.individual[1]; 647 current.gen = vector.individual[1];
640 if (vector.individual[3] > 68) { 648 if (vector.individual[3] > 68) {
641 # Add 68C to become fully reproductively matured. 649 # Add 68C to become fully reproductively matured.
642 # Transfer to vittelogenic. 650 # Transfer to vittelogenic.
643 vector.individual = c(current.gen, 4, 0, 0, 0); 651 vector.individual = c(current.gen, 4, 0, 0, 0);
748 # Add 1 day in current stage. 756 # Add 1 day in current stage.
749 vector.individual[4] = vector.individual[4] + 1; 757 vector.individual[4] = vector.individual[4] + 1;
750 } 758 }
751 vector.matrix[i,] = vector.individual; 759 vector.matrix[i,] = vector.individual;
752 } 760 }
753 # Old nymph to adult: previttelogenic or diapausing? 761 # Old nymph to adult: pre-vittelogenic or diapausing?
754 if (vector.individual[2] == 2) { 762 if (vector.individual[2] == 2) {
755 # Add average temperature for current day. 763 # Add average temperature for current day.
756 vector.individual[3] = vector.individual[3] + averages.temp; 764 vector.individual[3] = vector.individual[3] + averages.temp;
757 if (vector.individual[3] >= (200+opt$adult_accumulation)) { 765 if (vector.individual[3] >= (200+opt$adult_accumulation)) {
758 # From old to adult, degree_days requirement met. 766 # From old to adult, degree_days requirement met.
1039 } 1047 }
1040 1048
1041 if (process_eggs) { 1049 if (process_eggs) {
1042 # Mean value for eggs. 1050 # Mean value for eggs.
1043 eggs = apply(Eggs.replications, 1, mean); 1051 eggs = apply(Eggs.replications, 1, mean);
1052 temperature_data_frame = append_vector(temperature_data_frame, eggs, "EGG");
1044 # Standard error for eggs. 1053 # Standard error for eggs.
1045 eggs.std_error = apply(Eggs.replications, 1, sd) / sqrt(opt$replications); 1054 eggs.std_error = apply(Eggs.replications, 1, sd) / sqrt(opt$replications);
1055 temperature_data_frame = append_vector(temperature_data_frame, eggs.std_error, "EGGSE");
1046 } 1056 }
1047 if (process_nymphs) { 1057 if (process_nymphs) {
1048 # Calculate nymph populations for selected life stage. 1058 # Calculate nymph populations for selected life stage.
1049 for (life_stage_nymph in life_stages_nymph) { 1059 for (life_stage_nymph in life_stages_nymph) {
1050 if (life_stage_nymph=="Total") { 1060 if (life_stage_nymph=="Total") {
1051 # Mean value for all nymphs. 1061 # Mean value for all nymphs.
1052 total_nymphs = apply((YoungNymphs.replications+OldNymphs.replications), 1, mean); 1062 total_nymphs = apply((YoungNymphs.replications+OldNymphs.replications), 1, mean);
1063 temperature_data_frame = append_vector(temperature_data_frame, total_nymphs, "TOTALNYMPH");
1053 # Standard error for all nymphs. 1064 # Standard error for all nymphs.
1054 total_nymphs.std_error = apply((YoungNymphs.replications+OldNymphs.replications) / sqrt(opt$replications), 1, sd); 1065 total_nymphs.std_error = apply((YoungNymphs.replications+OldNymphs.replications) / sqrt(opt$replications), 1, sd);
1066 temperature_data_frame = append_vector(temperature_data_frame, total_nymphs.std_error, "TOTALNYMPHSE");
1055 } else if (life_stage_nymph=="Young") { 1067 } else if (life_stage_nymph=="Young") {
1056 # Mean value for young nymphs. 1068 # Mean value for young nymphs.
1057 young_nymphs = apply(YoungNymphs.replications, 1, mean); 1069 young_nymphs = apply(YoungNymphs.replications, 1, mean);
1070 temperature_data_frame = append_vector(temperature_data_frame, young_nymphs, "YOUNGNYMPH");
1058 # Standard error for young nymphs. 1071 # Standard error for young nymphs.
1059 young_nymphs.std_error = apply(YoungNymphs.replications / sqrt(opt$replications), 1, sd); 1072 young_nymphs.std_error = apply(YoungNymphs.replications / sqrt(opt$replications), 1, sd);
1073 temperature_data_frame = append_vector(temperature_data_frame, young_nymphs.std_error, "YOUNGNYMPHSE");
1060 } else if (life_stage_nymph=="Old") { 1074 } else if (life_stage_nymph=="Old") {
1061 # Mean value for old nymphs. 1075 # Mean value for old nymphs.
1062 old_nymphs = apply(OldNymphs.replications, 1, mean); 1076 old_nymphs = apply(OldNymphs.replications, 1, mean);
1077 temperature_data_frame = append_vector(temperature_data_frame, old_nymphs, "OLDNYMPH");
1063 # Standard error for old nymphs. 1078 # Standard error for old nymphs.
1064 old_nymphs.std_error = apply(OldNymphs.replications / sqrt(opt$replications), 1, sd); 1079 old_nymphs.std_error = apply(OldNymphs.replications / sqrt(opt$replications), 1, sd);
1080 temperature_data_frame = append_vector(temperature_data_frame, old_nymphs.std_error, "OLDNYMPHSE");
1065 } 1081 }
1066 } 1082 }
1067 } 1083 }
1068 if (process_adults) { 1084 if (process_adults) {
1069 # Calculate adult populations for selected life stage. 1085 # Calculate adult populations for selected life stage.
1070 for (life_stage_adult in life_stages_adult) { 1086 for (life_stage_adult in life_stages_adult) {
1071 if (life_stage_adult=="Total") { 1087 if (life_stage_adult=="Total") {
1072 # Mean value for all adults. 1088 # Mean value for all adults.
1073 total_adults = apply((Previttelogenic.replications+Vittelogenic.replications+Diapausing.replications), 1, mean); 1089 total_adults = apply((Previttelogenic.replications+Vittelogenic.replications+Diapausing.replications), 1, mean);
1090 temperature_data_frame = append_vector(temperature_data_frame, total_adults, "TOTALADULT");
1074 # Standard error for all adults. 1091 # Standard error for all adults.
1075 total_adults.std_error = apply((Previttelogenic.replications+Vittelogenic.replications+Diapausing.replications), 1, sd) / sqrt(opt$replications); 1092 total_adults.std_error = apply((Previttelogenic.replications+Vittelogenic.replications+Diapausing.replications), 1, sd) / sqrt(opt$replications);
1093 temperature_data_frame = append_vector(temperature_data_frame, total_adults.std_error, "TOTALADULTSE");
1076 } else if (life_stage_adult == "Pre-vittelogenic") { 1094 } else if (life_stage_adult == "Pre-vittelogenic") {
1077 # Mean value for previttelogenic adults. 1095 # Mean value for previttelogenic adults.
1078 previttelogenic_adults = apply(Previttelogenic.replications, 1, mean); 1096 previttelogenic_adults = apply(Previttelogenic.replications, 1, mean);
1097 temperature_data_frame = append_vector(temperature_data_frame, previttelogenic_adults, "PRE-VITADULT");
1079 # Standard error for previttelogenic adults. 1098 # Standard error for previttelogenic adults.
1080 previttelogenic_adults.std_error = apply(Previttelogenic.replications, 1, sd) / sqrt(opt$replications); 1099 previttelogenic_adults.std_error = apply(Previttelogenic.replications, 1, sd) / sqrt(opt$replications);
1100 temperature_data_frame = append_vector(temperature_data_frame, previttelogenic_adults.std_error, "PRE-VITADULTSE");
1081 } else if (life_stage_adult == "Vittelogenic") { 1101 } else if (life_stage_adult == "Vittelogenic") {
1082 # Mean value for vittelogenic adults. 1102 # Mean value for vittelogenic adults.
1083 vittelogenic_adults = apply(Vittelogenic.replications, 1, mean); 1103 vittelogenic_adults = apply(Vittelogenic.replications, 1, mean);
1104 temperature_data_frame = append_vector(temperature_data_frame, vittelogenic_adults, "VITADULT");
1084 # Standard error for vittelogenic adults. 1105 # Standard error for vittelogenic adults.
1085 vittelogenic_adults.std_error = apply(Vittelogenic.replications, 1, sd) / sqrt(opt$replications); 1106 vittelogenic_adults.std_error = apply(Vittelogenic.replications, 1, sd) / sqrt(opt$replications);
1107 temperature_data_frame = append_vector(temperature_data_frame, vittelogenic_adults.std_error, "VITADULTSE");
1086 } else if (life_stage_adult == "Diapausing") { 1108 } else if (life_stage_adult == "Diapausing") {
1087 # Mean value for vittelogenic adults. 1109 # Mean value for vittelogenic adults.
1088 diapausing_adults = apply(Diapausing.replications, 1, mean); 1110 diapausing_adults = apply(Diapausing.replications, 1, mean);
1111 temperature_data_frame = append_vector(temperature_data_frame, diapausing_adults, "DIAPAUSINGADULT");
1089 # Standard error for vittelogenic adults. 1112 # Standard error for vittelogenic adults.
1090 diapausing_adults.std_error = apply(Diapausing.replications, 1, sd) / sqrt(opt$replications); 1113 diapausing_adults.std_error = apply(Diapausing.replications, 1, sd) / sqrt(opt$replications);
1114 temperature_data_frame = append_vector(temperature_data_frame, diapausing_adults.std_error, "DIAPAUSINGADULTSE");
1091 } 1115 }
1092 } 1116 }
1093 } 1117 }
1094 1118
1095 if (plot_generations_separately) { 1119 if (plot_generations_separately) {
1172 F2_total_adults = m_se[[5]]; 1196 F2_total_adults = m_se[[5]];
1173 F2_total_adults.std_error = m_se[[6]]; 1197 F2_total_adults.std_error = m_se[[6]];
1174 } 1198 }
1175 } 1199 }
1176 1200
1201 # Save the analyzed data.
1202 write.csv(temperature_data_frame, file=opt$output, row.names=F);
1177 # Display the total number of days in the Galaxy history item blurb. 1203 # Display the total number of days in the Galaxy history item blurb.
1178 cat("Number of days: ", opt$num_days, "\n"); 1204 cat("Number of days: ", opt$num_days, "\n");
1179
1180 # Information needed for plots plots. 1205 # Information needed for plots plots.
1181 days = c(1:opt$num_days); 1206 days = c(1:opt$num_days);
1182 start_date = temperature_data_frame$DATE[1]; 1207 start_date = temperature_data_frame$DATE[1];
1183 end_date = temperature_data_frame$DATE[opt$num_days]; 1208 end_date = temperature_data_frame$DATE[opt$num_days];
1184 1209