# R code to accompany Fletcher, Webber, and Willis (2020) "Modelling the potential efficacy of treatments for white-nose syndrome in bats" # Author: Quinn Fletcher # Contact: q.fletcher@gmail.com # Date: 2020-03-10 #### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### # READ ME: This appendix contains the standard R script that recreates the figures presented in the three Shiny apps. These Shiny apps display the results from Models 1, 2a, and 2b described in the manuscript. The script for each model is structured in the same manner. The script indicates the beginning of each model (e.g. "MODEL 1 – GROWTH RATE OF A TREATED POPULATION"), followed by the libraries that need to be loaded for each model, then the text associated with each model that is presented in the Shiny apps. The next section for each model is where users can enter the scenario values that they would like to model (default values have been provided). Running the code for each model will produce .png figures that recreate those presented in the R Shiny Apps. The script for each model is self-contained, and thus, can be run independently from the script for the other models. We have included comments for the major sections of each model. We assume a general proficiency in R to run or modify the script. #### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### #### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### # MODEL 1 – GROWTH RATE OF A TREATED POPULATION # Load required libraries library(popbio) library(viridis) # Shiny app text: # "Outcome for a treated WNS-affected population (Model 1 in Fletcher et al. (LINK-20XX)): The expected outcome of applying a treatment to a white-nose syndrome (WNS) affected population. The left figure presents a heat map showing what our model predicts the geometric population growth rate (i.e. λ) of a treated WNS-affected population would be (λtreat) resulting from combinations of WNS severity (i.e. λWNS - λ of the WNS-affected population without treatment) and treatment improvement in winter survival (TIWS). Both λWNS and the TIWS are specified using the sliders. The black point in the left figure reflects the combination of the two slider-specified values. The right figure presents a population projection (max λtreat = 1.29; see Appendix S2) of a hypothetical population of 1000 bats, based on the values of WNS severity (λWNS) and the TIWS specified using the sliders." #### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ######## TWO USER-SPECIFIED SCENARIO VALUES #### # User-specified scenario of treatment improvement in winter survival (tiws) and white-nose syndrome severity (population growth rate of WNS-affected populations; wnslambda) values mod1_scenario_tiws <- 0.4 # value between 0 and 1 (increment=0.02); e.g. 0.38, 0.4 or 1, but not 0.41 mod1_scenario_wnslambda <- 0.5 # value between 0.02 and 1 (increment=0.02); e.g. 0.48, 0.5, or 1, but not 0.51 #### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### # Creation of the Lefkovitch matrix look-up table (data frame = matlookup) #### # Quote from Materials and methods: "To run our models, we needed a look-up table specifying Lefkovitch matrix parameters (Rj, Ra, Sj, Sa; Table 2) associated with specific levels of λ. Assumptions of our look-up table (Table 2) were based on a 16-year mark-recapture study of little brown bats pre-WNS (Frick et al. 2010b) and we generally followed Frick et al.’s (2010a) methodology to create the table (Appendix S1)." matlookup <- data.frame(cbind(bj = 0.38, sj = seq(0.01, 0.997, 0.007), sa = seq(0.01, 0.997, 0.007), f = 0.95, ba = 1, lambda = seq(0.01, 0.997, 0.007))) matlookup$sj <- matlookup$sa * 0.47 matlookup$rj <- matlookup$sj * matlookup$bj * matlookup$f matlookup$ra <- matlookup$sa * matlookup$ba * matlookup$f for (i in 1:nrow(matlookup)) { tempmat <- matrix(nrow = 2, ncol = 2, byrow = TRUE, data = c(matlookup$rj[i], matlookup$ra[i], matlookup$sj[i], matlookup$sa[i])) matlookup$lambda[i] <- lambda(tempmat) } # Rounding the look-up lambda values to two decimal places matlookup$lambda_round <- round(matlookup$lambda, 2) # Omitting lines with duplicate $lambda_round values # Quote from Appendix S1: "For the look-up table, we only wanted one record for each λ value; therefore, we omitted the duplicate λ record with the greater Sa and Sj values." matlookup$lambda_round_shift <- c(999, matlookup$lambda_round[1:(length(matlookup$lambda_round))-1]) # "999" is a place filler matlookup$lambda_round_compare <- ifelse(matlookup$lambda_round == matlookup$lambda_round_shift, 1, 0) matlookup <- subset(matlookup, matlookup$lambda_round_compare == 0) matlookup$lambda_round_factor <- as.factor(matlookup$lambda_round) matlookup <- subset(matlookup, select = c(bj, sj, sa, f, ba, rj, ra, lambda, lambda_round, lambda_round_factor)) # A value of 1.29 was set as the maximum growth rate of a treated WNS-affected population (λtreat) associated with combinations of treatment improvement in winter survival (tiws) and white-nose syndrome severity (wnslambda); see Appendix S2 #### growth_limit <- 1.29 # Data frame = all (derived from allsource). All scenarios of white-nose syndrome severity (wnslambda) and treatment improvement in winter survival (tiws) tested in this model. The scenario where wnslambda equals zero has been removed from this model. This data frame also includes the Lefkovitch matrix parameters (rj, ra, sj, and sa) extracted from the look-up table associated with the scenario-specific white-nose syndrome severity (wnslambda) as well as the Lefkovitch matrix parameters that would result from the treatment improvement in winter survival - upd_rj, upd_ra, upd_sj, and upd_sa #### sequence <- seq(0, 1, 0.02) wnslambda <- rep(sequence, 51) tiws <- rep(sequence, each = 51) allsource <- data.frame(wnslambda, tiws) allsource <- subset(allsource, wnslambda != 0) allsource$wnslambda_factor <- as.factor(allsource$wnslambda) all <- merge(allsource, matlookup, by.x = "wnslambda_factor", by.y = "lambda_round_factor", all.x = TRUE) all <- subset(all, select = c(wnslambda, tiws, sj, sa, bj, ba, f)) all <- all[order(all$wnslambda, all$tiws), ] all$upd_sj <- all$sj + all$tiws all$upd_sa <- all$sa + all$tiws all$upd_rj <- all$upd_sj * all$bj * all$f all$upd_ra <- all$upd_sa * all$ba * all$f # The growth rate of a treated WNS-affected population (trtlambda) associated with each model scenario #### all$trtlambda <- as.numeric(NA) for (i in 1:nrow(all)) { tempmat <- matrix(nrow = 2, ncol = 2, byrow = TRUE, data = c(all$upd_rj[i], all$upd_ra[i], all$upd_sj[i], all$upd_sa[i])) all$trtlambda[i] <- lambda(tempmat) } # If trtlambda was greater than one, it was rounded to one. #### # Quote from Figure 1 caption: "All λtreat values greater than 1.00 were rounded to 1.00." all$trtlambda_round1 <- ifelse(all$trtlambda > 1, 1, all$trtlambda) # Steps required to plot the trtlambda_round1 values in the heatmap as an image #### hm <- all$trtlambda_round1 hm <- matrix(hm, nrow = 50, ncol = 51, byrow = T) hm <- t(hm) # Heatmap images range from 0 to 1 on the x-axis and y-axis. This is not a problem for the x-axis because the variable it reflects (tiws) ranges from 0 to 1. However, this is a problem for the y-axis because the variable it reflects (wnslambda) ranges from 0.02 to 1. As a result, the y-axis needs to be rescaled so that the heatmap image range (i.e. 0 to 1) reflects the range of wnslambda (0.02 to 1). We determined the slope and intercept of the relationship (lm1) between the heatmap image endpoints (0 and 1; image_val_end) and the wnslambda end points (0.02 and 1; wnslambda_end) and then used this relationship to determine the heatmap image values (image_val) where we needed to put our rescaled axis tick marks (data frame = ylabels). #### image_val_end <- c(0, 1) rescaled_end <- c(0.02, 1) lm1 <- lm(image_val_end ~ rescaled_end) ylabels <- data.frame(wnslambda_val = c(0.02, 0.2, 0.4, 0.6, 0.8, 1.0)) ylabels$image_val <- coef(lm1)[1] + ylabels$wnslambda_val * coef(lm1)[2] # Merging the user-specified model 1 scenario (data frame = mod1_scenario) with the data frame containing all model 1 scenarios (data frame = all) #### mod1_scenario = data.frame(cbind(mod1_scenario_tiws, mod1_scenario_wnslambda)) mod1_scenario <- merge(mod1_scenario, all, by.x = c("mod1_scenario_wnslambda", "mod1_scenario_tiws"), by.y = c("wnslambda", "tiws"), all.x = T) # If the growth rate of the treated WNS-affected population (trtlambda) was greater than 1.29, then trtlambda value reflected in the projection figure was set to 1.29 (growth_limit). #### mod1_scenario$trtlambda <- ifelse(mod1_scenario$trtlambda > growth_limit, growth_limit, mod1_scenario$trtlambda) # Model 1 heatmap figure #### png("Model1_Heatmap.png", width = 1200, height = 1000, res = 200, units = "px") layout(matrix(c(1,2), 1, 2, byrow = TRUE), c(6,2), c(1,1)) par(mar = c(5.1, 5.1, 4.1, 1.1)) image(hm, col = viridis(100), las = 1, ylim = c(-0.03, 1.03), yaxt = "n", xlim = c(-0.04, 1.04), xaxp = c(0, 1, 5), cex.axis = 1.1) # Rescaled axis axis(side = 2, at = c(ylabels$image_val[1], ylabels$image_val[2], ylabels$image_val[3], ylabels$image_val[4], ylabels$image_val[5], ylabels$image_val[6]), labels = c("0.02", "0.2", "0.4", "0.6", "0.8", "1.0"), las = 1, cex.axis = 1.1) # Axis labels mtext(expression(paste("WNS severity (",lambda[WNS], ")")), side = 2, line = 3.1, cex = 1.1) mtext("Treatment improvement in winter survival", side = 1, line = 2.5, cex = 1.1, at = 0.5) # Rescaled wnslambda point rescaled_wnslambda <- coef(lm1)[1] + mod1_scenario$mod1_scenario_wnslambda * coef(lm1)[2] points(rescaled_wnslambda ~ mod1_scenario$mod1_scenario_tiws, pch = 15, cex = 1.5) # Scale bar par(mar = c(5.1, 1.1, 4.1, 5.1)) zz <- matrix(1:100,nrow = 1) xx <- 1 yy <- seq(from = 0, to = 1, length.out = 100) image(xx, yy, zz, col = viridis(100), axes = FALSE, xlab = "", ylab = "", cex.axis = 1.1) box() axis(side = 4, at = c(0, 0.25, 0.5, 0.75, 1), labels = c("0", "0.25", "0.5", "0.75", "≥1"), las = 2, cex.axis = 1.1) mtext(side = 4, expression(paste(lambda, " of treated population (", lambda[treat], ")")), line = 3.9, cex = 1.1) dev.off() # Population projection (max trtlambda = 1.29) of a hypothetical initial population of 1000 bats (ipop), based on user-specified values of white-nose syndrome severity (wnslambda) and treatment improvement in winter survival (tiws) #### ipop <- 1000 numyears <- 6 projection <- matrix(data = ipop, nrow = numyears, ncol = 1) for (i in 2:(numyears)) { projection[i,1] <- projection[i-1, 1] * mod1_scenario$trtlambda } # Model 1 population projection figure #### png("Model1_Projection.png", width = 1000, height = 1000, res = 200, units = "px") par(mar = c(5.1, 5.1, 4.1, 1.1)) plot(projection[c(1:numyears), 1] ~ c(1:numyears), type = "o", las = 1, cex.axis = 1.1, ylab = "", xlab = "", xaxt = "n", pch = 15, cex = 1.5, ylim = c(0, max(projection))) mtext("Number of bats", side = 2, line = 3.8, cex = 1.1) mtext("Years", side = 1, line = 2.5, cex = 1.1) axis(1, at = c(1, 2, 3, 4, 5, 6), labels = c("0", "1", "2", "3", "4", "5"), cex.axis = 1.1) dev.off() #### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### #### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### # MODEL 2A - TREATING A PORTION OF A POPULATION (CONSTANT WNS SEVERITY) # Load required libraries library(popbio) # Shiny app text: Outcome for a treated WNS-affected population if only a proportion of the population is treated (Model 2a in Fletcher et al. (LINK-20XX)): The predicted outcome of applying a treatment to a white-nose syndrome (WNS) affected population when only a proportion of the population is treated and WNS severity (i.e. λWNS - λ of the WNS-affected population without treatment) remains constant over time. The following 4 values can be specified using the sliders: 1) initial total population size, 2) λWNS, 3) proportion of the population treated, and 4) the treatment improvement in winter survival (TIWS). Population projections are given for both the treated and untreated portions of the population, as well as for the total population. #### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ######## FOUR USER-SPECIFIED SCENARIO VALUES #### # User-specified scenario of initial population size (ipop), white-nose syndrome severity (population growth rate of WNS-affected populations; wnslambda), proportion of the population treated (prop_treated), and treatment improvement in winter survival (tiws) mod2a_scenario_ipop <- 1000 # value between 500 and 10000 (increment=250); e.g. 500, 750, 8000, but not 1100 mod2a_scenario_wnslambda <- 0.5 # value between 0.02 and 1.0 (increment=0.02); e.g. 0.48, 0.5, or 1, but not 0.51 mod2a_scenario_prop_treated <- 0.5 # value between 0.02 and 0.98 (increment=0.02); e.g. 0.50 or 0.6, but not 0.61 mod2a_scenario_tiws <- 0.4 # value between 0 and 1 (increment=0.02); e.g. 0.38, 0.4 or 1, but not 0.41 #### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### # Creation of the Lefkovitch matrix look-up table (data frame = matlookup) #### # Quote from Materials and methods: "To run our models, we needed a look-up table specifying Lefkovitch matrix parameters (Rj, Ra, Sj, Sa; Table 2) associated with specific levels of λ. Assumptions of our look-up table (Table 2) were based on a 16-year mark-recapture study of little brown bats pre-WNS (Frick et al. 2010b) and we generally followed Frick et al.’s (2010a) methodology to create the table (Appendix S1)." matlookup <- data.frame(cbind(bj = 0.38, sj = seq(0.01, 0.997, 0.007), sa = seq(0.01, 0.997, 0.007), f = 0.95, ba = 1, lambda = seq(0.01, 0.997, 0.007))) matlookup$sj <- matlookup$sa * 0.47 matlookup$rj <- matlookup$sj * matlookup$bj * matlookup$f matlookup$ra <- matlookup$sa * matlookup$ba * matlookup$f for (i in 1:nrow(matlookup)) { tempmat <- matrix(nrow = 2, ncol = 2, byrow = TRUE, data = c(matlookup$rj[i], matlookup$ra[i], matlookup$sj[i], matlookup$sa[i])) matlookup$lambda[i] <- lambda(tempmat) } # Rounding the look-up lambda values to two decimal places matlookup$lambda_round <- round(matlookup$lambda, 2) # Omitting lines with duplicate $lambda_round values # Quote from Appendix S1: "For the look-up table, we only wanted one record for each λ value; therefore, we omitted the duplicate λ record with the greater Sa and Sj values." matlookup$lambda_round_shift <- c(999, matlookup$lambda_round[1:(length(matlookup$lambda_round))-1]) # "999" is a place filler matlookup$lambda_round_compare <- ifelse(matlookup$lambda_round == matlookup$lambda_round_shift, 1, 0) matlookup <- subset(matlookup, matlookup$lambda_round_compare == 0) matlookup$lambda_round_factor <- as.factor(matlookup$lambda_round) matlookup <- subset(matlookup, select = c(bj, sj, sa, f, ba, rj, ra, lambda, lambda_round, lambda_round_factor)) # A value of 1.29 was set as the maximum growth rate of a treated WNS-affected population (λtreat) associated with combinations of treatment improvement in winter survival (tiws) and white-nose syndrome severity (wnslambda); see Appendix S2 #### growth_limit <- 1.29 # The Lefkovitch matrix parameters extracted from the look-up table associated with the maximum value for the growth rate of a treated WNS-affected population (1.29; λtreat). #### growth_limit_vals <- subset(matlookup, lambda_round_factor == "1.29", select = c(rj, ra, sj, sa)) matrix_vals_for_growth_limit <- c(growth_limit_vals$rj, growth_limit_vals$ra, growth_limit_vals$sj, growth_limit_vals$sa) # Extracting the Lefkovitch matrix parameters from the look-up table based on the user-specified value of white-nose syndrome severity (wnslambda) #### mod2a_scenario <- data.frame(cbind(mod2a_scenario_ipop, mod2a_scenario_wnslambda, mod2a_scenario_prop_treated, mod2a_scenario_tiws)) mod2a_scenario$mod2a_scenario_wnslambda_factor <- as.factor(mod2a_scenario$mod2a_scenario_wnslambda) mod2a_scenario <- merge(mod2a_scenario, matlookup, by.x = c("mod2a_scenario_wnslambda"), by.y = c("lambda_round_factor"), all.x = T) # Run for 5 years, i.e. 6 time-steps #### numyears <- 6 # Seeding and starting population projection based on user-specified values #### # Seed abundance matrix (i.e. matrix of NAs) pooling the untreated and treated portions at the beginning of each projection year (juvenile: row 1; adult: row 2). beg_pool_abund <- matrix(NA, nrow = 2, ncol = numyears) # Pooled (treated and untreated) abundance matrix for juvenile (row 1) and adult (row 2) bats at the beginning of year 1 beg_pool_abund[1,1] <- mod2a_scenario_ipop * 0.44 # juv beg_pool_abund[2,1] <- mod2a_scenario_ipop * 0.56 # adult # Seed abundance matrices (i.e. matrices of NAs) for the untreated (untrt) and treated (trt) portions at the beginning of each projection year (juvenile: row 1; adult: row 2). beg_untrt_abund <- matrix(NA, nrow = 2, ncol = numyears) beg_trt_abund <- matrix(NA, nrow = 2, ncol = numyears) # Seed abundance matrices (i.e. matrix of NAs) for the untreated (untrt) and treated (trt) portions at the end of each projection year (juvenile: row 1; adult: row 2). end_untrt_abund <- matrix(NA, nrow = 2, ncol = numyears-1) end_trt_abund <- matrix(NA, nrow = 2, ncol = numyears-1) # User-defined scenario Lefkovitch matrix for the untreated portion tempmat_untrt <- matrix(nrow = 2, ncol = 2, byrow = TRUE, data = c(mod2a_scenario$rj, mod2a_scenario$ra, mod2a_scenario$sj, mod2a_scenario$sa)) # Year 1 abundance of juvenile (row 1) and adult (row 2) bats at the beginning of the year for the untreated (untrt) and treated (trt) portions. beg_untrt_abund[1, 1] <- beg_pool_abund[1, 1] * (1 - mod2a_scenario$mod2a_scenario_prop_treated) # juv beg_untrt_abund[2, 1] <- beg_pool_abund[2, 1] * (1 - mod2a_scenario$mod2a_scenario_prop_treated) # adult beg_trt_abund[1, 1] <- beg_pool_abund[1, 1] * (mod2a_scenario$mod2a_scenario_prop_treated) # juv beg_trt_abund[2,1] <- beg_pool_abund[2, 1] * (mod2a_scenario$mod2a_scenario_prop_treated) # adult # Finishing population projection based on user-specified values #### for (j in 2:(numyears)) { # Abundance of juvenile (row 1) and adult (row 2) bats at the end of the year j-1 for the untreated (untrt) portion end_untrt_abund [, j-1] <- tempmat_untrt %*% beg_untrt_abund[, j-1] # User-defined scenario Lefkovitch matrix for the treated portion. If the lambda of this matrix is greater than 1.29 (i.e. growth_limit; see Appendix S2), then assign the scenario Lefkovitch matrix to be the matrix with lambda equal to 1.29 ("matrix_vals_for_growth_limit") test_tempmat_trt <- matrix(nrow = 2, ncol = 2, byrow = TRUE, data = c((mod2a_scenario$sj + mod2a_scenario_tiws) * mod2a_scenario$bj * mod2a_scenario$f, (mod2a_scenario$sa + mod2a_scenario_tiws) * mod2a_scenario$ba * mod2a_scenario$f, (mod2a_scenario$sj + mod2a_scenario_tiws), (mod2a_scenario$sa + mod2a_scenario_tiws))) lambdatest <- lambda(test_tempmat_trt) if (lambdatest > growth_limit) { tempmat_trt <- matrix(nrow = 2, ncol = 2, byrow = TRUE, data = matrix_vals_for_growth_limit) } else { tempmat_trt <- matrix(nrow = 2, ncol = 2, byrow = TRUE, data = c((mod2a_scenario$sj + mod2a_scenario_tiws) * mod2a_scenario$bj * mod2a_scenario$f, (mod2a_scenario$sa + mod2a_scenario_tiws) * mod2a_scenario$ba * mod2a_scenario$f, (mod2a_scenario$sj + mod2a_scenario_tiws), (mod2a_scenario$sa + mod2a_scenario_tiws))) } # Abundance of juvenile (row 1) and adult (row 2) bats at the end of the year j-1 for the treated (trt) portion end_trt_abund[, j-1] <- tempmat_trt %*% beg_trt_abund[, j-1] # Pooled (treated and untreated) abundance of juvenile (row 1) and adult (row 2) bats at the beginning of year j beg_pool_abund[1, j] <- end_trt_abund[1, j-1] + end_untrt_abund[1, j-1] # juv beg_pool_abund[2, j] <- end_trt_abund[2, j-1] + end_untrt_abund[2, j-1] # adult # Abundance of juvenile (row 1) and adult (row 2) bats at the beginning of year j for untreated (untrt) and treated (trt) portions beg_untrt_abund[1, j] <- beg_pool_abund[1, j] * (1-mod2a_scenario$mod2a_scenario_prop_treated) # juv beg_untrt_abund[2, j] <- beg_pool_abund[2, j] * (1-mod2a_scenario$mod2a_scenario_prop_treated) # adult beg_trt_abund[1, j] <- beg_pool_abund[1, j] * (mod2a_scenario$mod2a_scenario_prop_treated) # juv beg_trt_abund[2, j] <- beg_pool_abund[2, j] * (mod2a_scenario$mod2a_scenario_prop_treated) # adult } # end OF j (i.e. Finishing population projection based on user-specified values) # Sum of juvenile plus adult abundances at the beginning (i.e. "beg_...") and end (i.e. "end_...") of each year for the untreated (untrt) and treated (trt) portions #### beg_untrt_abund_s <- colSums(beg_untrt_abund) end_untrt_abund_s <- colSums(end_untrt_abund) beg_trt_abund_s <- colSums(beg_trt_abund) end_trt_abund_s <- colSums(end_trt_abund) # Data formating for projection plot: treated portion (proj_trt) #### proj_trt <- c(beg_trt_abund_s[1], end_trt_abund_s[1]) for (i in 2:5) { temp1 <- beg_trt_abund_s[i] temp2 <- end_trt_abund_s[i] proj_trt <- c(proj_trt , temp1, temp2) } constant <- 0.2 yr <- c(0+constant, 1, 1+constant, 2, 2+constant, 3, 3+constant, 4, 4+constant, 5) proj_trt <- data.frame(yr = yr, proj_trt) # Data formating for projection plot: untreated portion (proj_untrt) #### proj_untrt <- c(beg_untrt_abund_s[1], end_untrt_abund_s[1]) for (i in 2:5) { temp1 <- beg_untrt_abund_s[i] temp2 <- end_untrt_abund_s[i] proj_untrt <- c(proj_untrt , temp1, temp2) } yr <- c(0+constant, 1, 1+constant, 2, 2+constant, 3, 3+constant, 4, 4+constant, 5) proj_untrt <- data.frame(yr = yr, proj_untrt) # Data formating for projection plot: making one data frame for untreated and treated portions (mod2a_scenario_plot) #### mod2a_scenario_plot <- cbind(proj_trt, proj_untrt[, 2]) names(mod2a_scenario_plot) <- c("yr", "proj_trt", "proj_untrt") mod2a_scenario_plot$tot <- mod2a_scenario_plot$proj_trt + mod2a_scenario_plot$proj_untrt mod2a_scenario_plot$tot_plot <- c(1, 1, 0, 1, 0, 1, 0, 1, 0, 1) mod2a_scenario_plot$lambda_pair <- c(1, 1, 2, 2, 3, 3, 4, 4, 5, 5) mod2a_scenario_plot$prop_treated_pair <- c(0, 1, 1, 2, 2, 3, 3, 4, 4, 0) maxv <- max(mod2a_scenario_plot$tot) # Model 2a population projection figure #### png("Model2a_Projection.png", width = 1200, height = 900, res = 200, units = "px") par(mar=c(4.1,5.6,0.5,6.0), xpd=TRUE) plot(tot ~ yr, mod2a_scenario_plot, type = "n", ylim=c(0, maxv), xlim = c(0, 5), las = 1, ylab = "", xlab = "", cex.axis = 1.1) # Treated portion abundance changes due to treatment for (i in 1:5) { temp <- subset(mod2a_scenario_plot, lambda_pair == i) points(proj_trt ~ yr, temp, pch = 15, lty = 1, type = "o") } # Treated portion abundance changes due to partioning population into untreated and treated portions for (i in 1:4) { temp <- subset(mod2a_scenario_plot, prop_treated_pair == i) points(proj_trt ~ yr, temp, pch = 15, lty = 3, type = "o") } # Untreated portion abundance changes due to white-nose syndrome severity for (i in 1:5) { temp <- subset(mod2a_scenario_plot, lambda_pair == i) points(proj_untrt ~ yr, temp, lty = 1, type = "o", pch = 1) } # Untreated portion abundance changes due to partioning population into untreated and treated portions for (i in 1:4) { temp <- subset(mod2a_scenario_plot, prop_treated_pair == i) points(proj_untrt ~ yr, temp, lty = 3, type = "o", pch = 1) } # Total population abundance (i.e. sum of untreated and treated portions) changes points(tot ~ yr, subset(mod2a_scenario_plot, tot_plot == 1), type = "o", pch = 19, lty = 1) legend(x = numyears-1+0.25, y = maxv, c("Total", "Treated", "Untreated"), pch=c(19, 15, 1)) mtext("Number of bats", side = 2, line = 4.5, cex = 1.2) mtext("Years", side = 1, line = 2.5, cex = 1.2) dev.off() #### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### #### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### # MODEL 2B - TREATING A PORTION OF A POPULATION (ATTENUATING WNS SEVERITY) # Load required libraries library(popbio) # Shiny app text: Outcome for a treated WNS-affected population if only a proportion of the population is treated and WNS severity attenutates (Model 2b in Fletcher et al. (LINK-20XX)): The predicted outcome of applying a treatment to a white-nose syndrome (WNS) affected population when only a proportion of the population is treated and WNS severity (i.e. λWNS - λ of the WNS-affected population without treatment) attenuates over time. We assumed that λWNS for both the treated and untreated portions of the population increased by a constant amount each year so that λWNS equaled 1 by year 5 post-invasion. The following 4 values can be specified using the sliders: 1) initial total population size, 2) λWNS, 3) proportion of the population treated, and 4) the treatment improvement in winter survival (TIWS). Population projections are given for both the treated and untreated portions of the population, as well as for the total population. #### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ######## FOUR USER-SPECIFIED SCENARIO VALUES #### # User-specified scenario of initial population size (ipop), white-nose syndrome severity (population growth rate of WNS-affected populations; wnslambda), proportion of the population treated (prop_treated), and treatment improvement in winter survival (tiws) mod2b_scenario_ipop <- 1000 # value between 500 and 10000 (increment=250); e.g. 500, 750, 8000, but not 1100 mod2b_scenario_wnslambda <- 0.5 # value between 0.02 and 1.0 (increment=0.02); e.g. 0.48, 0.5, or 1, but not 0.51 mod2b_scenario_prop_treated <- 0.3 # value between 0.02 and 0.98 (increment=0.02); e.g. 0.50 or 0.6, but not 0.61 mod2b_scenario_tiws <- 0.4 # value between 0 and 1 (increment=0.02); e.g. 0.38, 0.4 or 1, but not 0.41 #### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### # Creation of the Lefkovitch matrix look-up table (data frame = matlookup) #### # Quote from Materials and methods: "To run our models, we needed a look-up table specifying Lefkovitch matrix parameters (Rj, Ra, Sj, Sa; Table 2) associated with specific levels of λ. Assumptions of our look-up table (Table 2) were based on a 16-year mark-recapture study of little brown bats pre-WNS (Frick et al. 2010b) and we generally followed Frick et al.’s (2010a) methodology to create the table (Appendix S1)." matlookup <- data.frame(cbind(bj = 0.38, sj = seq(0.01,0.997,0.007), sa = seq(0.01,0.997,0.007), f = 0.95, ba = 1, lambda = seq(0.01,0.997,0.007))) matlookup$sj <- matlookup$sa * 0.47 matlookup$rj <- matlookup$sj * matlookup$bj * matlookup$f matlookup$ra <- matlookup$sa * matlookup$ba * matlookup$f for (i in 1:nrow(matlookup)) { tempmat <- matrix(nrow = 2, ncol = 2, byrow = TRUE, data = c(matlookup$rj[i], matlookup$ra[i], matlookup$sj[i], matlookup$sa[i])) matlookup$lambda[i] <- lambda(tempmat) } # Rounding the look-up lambda values to two decimal places matlookup$lambda_round <- round(matlookup$lambda, 2) # Omitting lines with duplicate $lambda_round values # Quote from Appendix S1: "For the look-up table, we only wanted one record for each λ value; therefore, we omitted the duplicate λ record with the greater Sa and Sj values." matlookup$lambda_round_shift <- c(999, matlookup$lambda_round[1:(length(matlookup$lambda_round))-1]) # "999" is a place filler matlookup$lambda_round_compare <- ifelse(matlookup$lambda_round == matlookup$lambda_round_shift, 1, 0) matlookup <- subset(matlookup, matlookup$lambda_round_compare == 0) matlookup$lambda_round_factor <- as.factor(matlookup$lambda_round) matlookup <- subset(matlookup, select = c(bj, sj, sa, f, ba, rj, ra, lambda, lambda_round, lambda_round_factor)) # A value of 1.29 was set as the maximum growth rate of a treated WNS-affected population (λtreat) associated with combinations of treatment improvement in winter survival (tiws) and white-nose syndrome severity (wnslambda); see Appendix S2 #### growth_limit <- 1.29 # The Lefkovitch matrix parameters extracted from the look-up table associated with the maximum value for the growth rate of a treated WNS-affected population (1.29; λtreat). #### growth_limit_vals <- subset(matlookup, lambda_round_factor == "1.29", select = c(rj, ra, sj, sa)) matrix_vals_for_growth_limit <- c(growth_limit_vals$rj, growth_limit_vals$ra, growth_limit_vals$sj, growth_limit_vals$sa) # Run for 5 years, i.e. 6 time-steps #### numyears <- 6 # Seeding and starting population projection based on user-specified values #### # Seed abundance matrix (i.e. matrix of NAs) pooling the untreated and treated portions at the beginning of each projection year (juvenile: row 1; adult: row 2). beg_pool_abund <- matrix(NA, nrow = 2, ncol = numyears) # Pooled (treated and untreated) abundance matrix for juvenile (row 1) and adult (row 2) bats at the beginning of year 1 beg_pool_abund[1,1] <- mod2b_scenario_ipop * 0.44 # juv beg_pool_abund[2,1] <- mod2b_scenario_ipop * 0.56 # adult # Seed abundance matrices (i.e. matrices of NAs) for the untreated (untrt) and treated (trt) portions at the beginning of each projection year (juvenile: row 1; adult: row 2). beg_untrt_abund <- matrix(NA, nrow = 2, ncol = numyears) beg_trt_abund <- matrix(NA, nrow = 2, ncol = numyears) # Seed abundance matrices (i.e. matrices of NAs) for the untreated (untrt) and treated (trt) portions at the end of each projection year (juvenile: row 1; adult: row 2). end_untrt_abund <- matrix(NA, nrow = 2, ncol = numyears-1) end_trt_abund <- matrix(NA, nrow = 2, ncol = numyears-1) # Data frame = red_severity. For all scenarios, this data frame allows the white-nose syndrome severity (wnslambda) for both the treated and untreated portions to increase by a constant amount each year such that the wnslambda equals one by year 5 years <- c(1:(numyears)) wnslambdavector <- as.vector(matrix(NA, nrow = 1, ncol = numyears-1)) wnslambdavector[1] <- mod2b_scenario_wnslambda for (k in 2:length(years)) { wnslambdavector[k] <- (1.0 - wnslambdavector[1]) / (numyears-2) + wnslambdavector[k-1] } red_severity <- data.frame(cbind(years, wnslambdavector)) red_severity$wnslambdafactor <- round(red_severity$wnslambdavector, 2) red_severity$wnslambdafactor <- as.factor(red_severity$wnslambdafactor) red_severity <- merge(red_severity, matlookup, by.x = "wnslambdafactor", by.y = "lambda_round_factor", all.x = TRUE) # Year 1 abundance of juvenile (row 1) and adult (row 2) bats at the beginning of the year for the untreated (untrt) and treated (trt) portions. beg_untrt_abund[1,1] <- beg_pool_abund[1,1] * (1 - mod2b_scenario_prop_treated) # juv beg_untrt_abund[2,1] <- beg_pool_abund[2,1] * (1 - mod2b_scenario_prop_treated) # adult beg_trt_abund[1,1] <- beg_pool_abund[1,1] * (mod2b_scenario_prop_treated) # juv beg_trt_abund[2,1] <- beg_pool_abund[2,1] * (mod2b_scenario_prop_treated) # adult # Finishing population projection based on user-specified values #### for (j in 2:numyears) { # User-defined scenario Lefkovitch matrix for the untreated portion incorporating the reducing wnslambda severity values specified in the data frame "red_severity" tempmat_untrt <- matrix(nrow = 2, ncol = 2, byrow = TRUE, data = c(red_severity$rj[j-1], red_severity$ra[j-1], red_severity$sj[j-1], red_severity$sa[j-1])) # Abundance of juvenile (row 1) and adult (row 2) bats at the end of the year j-1 for the untreated (untrt) portion end_untrt_abund [, j-1] <- tempmat_untrt %*% beg_untrt_abund[, j-1] # User-defined scenario Lefkovitch matrix for the treated portion incorporating the reducing severity wnslambda values specified in the data frame "red_severity". If the lambda of this matrix is greater than 1.29 (i.e. growth_limit; see Appendix S2), then assign the scenario Lefkovitch matrix to be the matrix with lambda equal to 1.29 ("matrix_vals_for_growth_limit") test_tempmat_trt <- matrix(nrow = 2, ncol = 2, byrow = TRUE, data = c((red_severity$sj[j-1] + mod2b_scenario_tiws) * red_severity$bj[j-1] * red_severity$f[j-1], (red_severity$sa[j-1] + mod2b_scenario_tiws) * red_severity$ba[j-1] * red_severity$f[j-1], (red_severity$sj[j-1] + mod2b_scenario_tiws), (red_severity$sa[j-1] + mod2b_scenario_tiws))) lambdatest <- lambda(test_tempmat_trt) if (lambdatest > growth_limit) { tempmat_trt <- matrix(nrow = 2, ncol = 2, byrow = TRUE, data = matrix_vals_for_growth_limit) } else { tempmat_trt <- matrix(nrow = 2, ncol = 2, byrow = TRUE, data = c((red_severity$sj[j-1] + mod2b_scenario_tiws) * red_severity$bj[j-1] * red_severity$f[j-1], (red_severity$sa[j-1] + mod2b_scenario_tiws) * red_severity$ba[j-1] * red_severity$f[j-1], (red_severity$sj[j-1] + mod2b_scenario_tiws), (red_severity$sa[j-1] + mod2b_scenario_tiws))) } # Abundance of juvenile (row 1) and adult (row 2) bats at the end of the year j-1 for the treated (trt) portion end_trt_abund [, j-1] <- tempmat_trt %*% beg_trt_abund[, j-1] # Pooled (treated and untreated) abundance of juvenile (row 1) and adult (row 2) bats at the beginning of year j beg_pool_abund[1, j] <- end_trt_abund[1, j-1] + end_untrt_abund[1, j-1] # juv beg_pool_abund[2, j] <- end_trt_abund[2, j-1] + end_untrt_abund[2, j-1] # adult # Abundance of juvenile (row 1) and adult (row 2) bats at the beginning of year j for untreated (untrt) and treated (trt) portions beg_untrt_abund[1, j] <- beg_pool_abund[1, j] * (1 - mod2b_scenario_prop_treated) # juv beg_untrt_abund[2, j] <- beg_pool_abund[2, j] * (1 - mod2b_scenario_prop_treated) # adult beg_trt_abund[1, j] <- beg_pool_abund[1, j] * (mod2b_scenario_prop_treated) # juv beg_trt_abund[2, j] <- beg_pool_abund[2, j] * (mod2b_scenario_prop_treated) # adult } # end OF j (i.e. Finishing population projection based on user-specified values) # Sum of juvenile plus adult abundances at the beginning (i.e. "beg_...") and end (i.e. "end_...") of each year for the untreated (untrt) and treated (trt) portions #### beg_trt_abund_s <- colSums(beg_trt_abund) end_trt_abund_s <- colSums(end_trt_abund) beg_untrt_abund_s <- colSums(beg_untrt_abund) end_untrt_abund_s <- colSums(end_untrt_abund) # Data formating for projection plot: treated portion (proj_trt) #### proj_trt <- c(beg_trt_abund_s[1], end_trt_abund_s[1]) for (i in 2:5) { temp1 <- beg_trt_abund_s[i] temp2 <- end_trt_abund_s[i] proj_trt <- c(proj_trt , temp1, temp2) } constant <- 0.2 yr <- c(0+constant, 1 , 1+constant, 2, 2+constant, 3, 3+constant, 4, 4+constant, 5) proj_trt <- data.frame(yr = yr, proj_trt) # Data formating for projection plot: untreated portion (proj_untrt) #### proj_untrt <- c(beg_untrt_abund_s[1], end_untrt_abund_s[1]) for (i in 2:5) { temp1 <- beg_untrt_abund_s[i] temp2 <- end_untrt_abund_s[i] proj_untrt <- c(proj_untrt , temp1, temp2) } yr <- c(0+constant, 1, 1+constant, 2, 2+constant, 3, 3+constant, 4, 4+constant, 5) proj_untrt <- data.frame(yr = yr, proj_untrt) # Data formating for projection plot: making one data frame for untreated and treated portions (mod2a_scenario_plot) #### mod2b_scenario_plot <- cbind(proj_trt, proj_untrt[,2]) names(mod2b_scenario_plot)=c("yr", "proj_trt", "proj_untrt") mod2b_scenario_plot$tot <- mod2b_scenario_plot$proj_trt + mod2b_scenario_plot$proj_untrt mod2b_scenario_plot$tot_plot <- c(1, 1, 0, 1, 0, 1, 0, 1, 0, 1) mod2b_scenario_plot$lambda_pair <- c(1, 1, 2, 2, 3, 3, 4, 4, 5, 5) mod2b_scenario_plot$prop_treated_pair <- c(0, 1, 1, 2, 2, 3, 3, 4, 4, 0) maxv <- max(mod2b_scenario_plot$tot) # Model 2b population projection figure #### png("Model2b_Projection.png", width = 1200, height = 900, res = 200, units = "px") par(mar=c(4.1,5.6,0.5,6.0), xpd=TRUE) plot(tot ~ yr, mod2b_scenario_plot, type = "n", ylim=c(0, maxv), xlim = c(0, 5), las = 1, ylab = "", xlab = "", cex.axis = 1.1) # Treated portion abundance changes due to treatment for (i in 1:5) { temp <- subset(mod2b_scenario_plot, lambda_pair == i) points(proj_trt ~ yr, temp, pch = 15, lty = 1, type = "o") } # Treated portion abundance changes due to partioning population into untreated and treated portions for (i in 1:4) { temp <- subset(mod2b_scenario_plot, prop_treated_pair == i) points(proj_trt ~ yr, temp, pch = 15, lty = 3, type = "o") } # Untreated portion abundance changes due to white-nose syndrome severity for (i in 1:5) { temp <- subset(mod2b_scenario_plot, lambda_pair == i) points(proj_untrt ~ yr, temp, lty = 1, type = "o", pch = 1) } # Untreated portion abundance changes due to partioning population into untreated and treated portions for (i in 1:4) { temp <- subset(mod2b_scenario_plot, prop_treated_pair == i) points(proj_untrt ~ yr, temp, lty = 3, type = "o", pch = 1) } # Total population abundance (i.e. sum of untreated and treated portions) changes points(tot ~ yr, subset(mod2b_scenario_plot, tot_plot == 1), type = "o", pch = 19, lty = 1) legend(x = numyears-1+0.25, y = maxv, c("Total", "Treated", "Untreated"), pch=c(19, 15, 1)) mtext("Number of bats", side = 2, line = 4.5, cex = 1.2) mtext("Years", side = 1, line = 2.5, cex = 1.2) dev.off()