### SURVIVAL FUNCTIONS ### ###################################################################### ###################################################################### ## Logistic mortality function with default ## parameters from Styer et al. The "mort" ## parameter determines whether exponential or ## logistic mortality survival curves are used. mort.day <- function(xx, aa=.0018, bb=.1416, cc=0, ss=1.0730, mort="log", mu=1/32.6694) { ans <- (aa * exp(bb*xx) ) / (1 + (aa*ss/bb)*(exp(bb*xx) - 1) ) + cc if(mort=="exp") ans <- mu ans } ## Survival function This function gives you the ## probability you live > "xx" days past day "start" ## given that you do live up to day "start" s2day <- function(xx, # xx is the age to live to (not number of days to live) start=0, mort="log", adulticide=0, mu = 1/32.6694, aa = .0018, bb=.1416, cc=0, ss=1.0730) { if(start<0) stop("Start time is negative!") if(xx<=0) ans <- 1 if(xx>=1) ans <- prod( (1-adulticide) * (1 - sapply( (start+1):(xx), mort.day, # Note from start+1 to xx mort = mort, mu = mu, aa = aa, bb = bb, cc = cc, ss = ss))) ans } ## Test that the probability of living up to day ## x, living up to day x-1 and dying on day x, and ## dying before day x all add up to 1 test.fxn.1 <- function(xx, start=0) { s2day(xx,start) + s2day(xx-1,start)*mort.day(xx+start) + (1 - s2day(xx - 1,start) ) } ## The following line does this test for x=1:10. ## sapply(1:10, test.fxn.1) ## In original version had a p.2sdays function ## that calculated the probability of living ## exactly s days past start given living up to ## day start. Not in here ## Expected number of days to live (life ## expectancy function). exp.days <- function(start, mort = "log", adulticide = 0, mu = 1/32.6694, aa = .0018, bb=.1416, cc=0, ss=1.0730) { ages.to.live.to <- 0:100 + start survival.curve <- sapply(ages.to.live.to, s2day, start = start, mort = mort, adulticide = adulticide, mu = mu, aa = aa, bb = bb, cc = cc, ss = ss) ans <- sum(survival.curve) ans } ## PLOT SURVIVAL FUNCTIONS ### ###################################################################### ###################################################################### ## Set working directory here. graphics.off() ## Set universal graphics parameters legend <- c("constant mortality", "age-dependent mortality") cex.lab <- 1 ## Plot survival function curves for logistic & ## exponential mortality functions ## This is central to the whole project. Plots will all be functions ## of "mort" which can either be constant (e.g. exponential mortality: ## "exp") or age-dependent (e.g. logistic mortality: "log"); ## "adulticide", "larvicide", and it is also possible to change shape ## parameters of the original survival curves ("aa", "ss", "cc", "bb", ## "mu") to explore the effects of different types of age ## distributions on everything. plot.surv <- function(mort = "log", # Type of mortality assumption adulticide = 0, # Daily mosquito mortality rate due to adult control larvicide = 0, # Proportional reduction in adult mosquito recruitment due to larval control. mu = 1/32.6694, # Next five parameters describe mosquito mortality rate functions, see text for details. aa = .0018, bb = .1416, cc = 0, ss = 1.0730, age.range = 1:80, # Age range over which to plot. normalize = TRUE, # Whether to normalize so that total mosquito population size, m, equals 1 when there is no control. pdf = FALSE, # Whether to make a pdf() of the plot inside the function. quartz = TRUE, # Whether to plot in quartz() inside the function (Mac only) new.plot = TRUE, # Whether to plot on a new plotting device. file.name = "survival curve.pdf", # Name of figure if saving as pdf. directory = "" # Where to save figure if saving. plot.width = 8, plot.height = 6, legend.true = FALSE, # Whether to plot the legend. legend.loc = "topright", # Where to plot the legend if plotting. legend = c("survival curve"), # Legend label if plotting. cex.lab = 1, # Axis label scaling parameter. col = "black", add = FALSE, # Whether to add curve to existing plot polygon = TRUE, # Whether to plot as a shaded polygon or not (just a line). label.yax = FALSE, # Whether to label y-axis label.xax = FALSE, # Whether to label x-axis xlim.auto = TRUE, # Whether to set x-axis limits automatically xlim = c(0,100), # X-axis limits if not setting automatically. ylim = c(0,1), xlab = "age class", ylab = "relative density", main = "", debug = FALSE) # Shows more detail for debugging. { begin.dir <- getwd() # Store Global working directory to reset at the end of function. setwd(directory) # Change working directory for saving. if(debug) {print("Changing or saving directory")} if(pdf) pdf(file = file.name, width = plot.width, height = plot.height) if(!pdf & new.plot & quartz) quartz(width = plot.width, height = plot.height) if(debug) {print("Initializing pdf")} surv.curv <- sapply(age.range, # Adult control incorporated in survival curves s2day, start = 0, mort = mort, adulticide = adulticide, mu = mu, aa = aa, bb = bb, cc = cc, ss = ss) if(debug) {print("Calculating survival curve")} if(normalize) { sum.surv.curv <- sum(surv.curv) if(adulticide == 0 & mort == "exp") save(sum.surv.curv, file="norm.scal.surv.exp") # Save the scalars for adulticide = 0 situations if(adulticide == 0 & mort == "log") save(sum.surv.curv, file="norm.scal.surv.log") # Then can load those for control scenarios. if(adulticide != 0 & mort == "exp") load(file="norm.scal.surv.exp") if(adulticide != 0 & mort == "log") load(file="norm.scal.surv.log") surv.curv <- surv.curv / sum.surv.curv ylim <- c(0,.04) } surv.curv.larv <- surv.curv * (1 - larvicide) # Larval control reduces all age-classes by a proportion. if(debug) {print("Calculating survival curve reduced by larvicide")} if(xlim.auto) xlim <- c(0, max(age.range)) # Sets the xlim to be the age range. if(new.plot) { plot(age.range, surv.curv.larv, type = "n", # Do not plot so that add can be done with points(). col = col, cex.lab = cex.lab, xlab = xlab, ylab = ylab, main = main, xlim = xlim, ylim = ylim, axes = F) ## Axes labels y.ax.ticks <- c(0,.5,1) # Where to label y axis. if(normalize) y.ax.ticks <- c(0,.02,.04) x.ax.ticks <- pretty(age.range, 5) y.ax.labels <- y.ax.ticks x.ax.labels <- x.ax.ticks if(!label.yax) y.ax.labels <- NA # Or don't label it if(!label.xax) x.ax.labels <- NA axis(1, at = x.ax.ticks, labels = x.ax.labels) axis(2, at = y.ax.ticks, labels = y.ax.labels) if(debug) {print("Finished with plot()")} } points(age.range, # This is where the actual points are put on. surv.curv.larv, type = "l", col = col) if(polygon) # This shades the survival curve in. { polygon(c(age.range, rev(age.range)), c(surv.curv.larv, rep(0, length(age.range))), col = col) } if(legend.true) { legend(legend.loc, legend, col = col) if(debug) {print("Finished with legend")} } if(pdf) dev.off() if(debug) {print("Finished turning graphics off")} setwd(begin.dir) if(debug) {print("Finished resetting directory")} } ## Age Class Vectorial Capacity: FUNCTIONS ########################### ###################################################################### ###################################################################### ## Set working directory here graphics.off() ## This script will plot scaled vectorial capacity contributions (C*) ## by age class as a function of "EIP", "mort", "adulticide", ## "larvicide", and all survival curve shape paraemters ## ("mu","aa","bb","cc","ss"). ## Recall we are looking at scaled vectorial capacity (c*=mp^n/-ln(p) ## for constant mortality model). This means that it is written as a ## function of relative mosquito density (scaled to m=1 in no control ## scenarios), p (the daily survival probability), and n (the EIP). ## In the age-dependent model p becomes a function of age and gets ## slightly more complicated but C* eventually gets reduced to a more ## simple equation which is: ## C* = sum(from x = n +1 to infinity of [ exp.days(x) * prod( from i =1 to x of p(i))] ## This says that C* merely equals the proportion of mosquitoes hwo ## live to n+1 days or older, mulitiplied by their life expectancy at ## that age, and summed across age classes. plot.vc <- function(mort = "log", # Type of mortality model, "exp" for constant mortality and "log" for age-dependent. EIP = 3, # Extrinsic incubation period adulticide = 0, # Daily probability of adult mosquito death due to adult control larvicide = 0, # Proportional reduction in adult mosquito recruitment due to larval control mu = 1/32.6694, # Next five parameters describe mortality curves, see text for detals aa = .0018, bb = .1416, cc = 0, ss = 1.0730, age.range = 1:80, # Age range over which to plot pdf = FALSE, # Whether to plot pdf in function. quartz = TRUE, # Whether to plot quartz in function. new.plot = TRUE, # Whether to plot figure on new plotting device file.name = "age-spec vc curve.pdf", # File name if saving figure as pdf. directory = "", # Where to save figure if saving. plot.width = 8, plot.height = 6, class.spec = TRUE, # If true plots C* by whole age class, if false plots per-individual C*. legend.true = FALSE, # Whether to plot legend. legend.loc = "topright", # Where to plot legend if plotting. legend = c("vectorial capacity curve"), # Legend label cex.lab = 1, # Scaling parameter for axes labels col = "black", add = FALSE, # Whether to add curve to existing plot polygon = TRUE, # Whether to plot as a shaded polygon (or if not as a line). lty = 1, xlim = c(0,100), ylim.auto = TRUE, # Whether to set y-axis limits automatically xlim.auto = TRUE, # Whether to set x-axis limits automatically label.yax = TRUE, # Whether to label y-axis. label.xax = TRUE, # Whether to label x-axis. ylim = c(0,3), # Y-axis limits if not setting automatically. xlab = "age class", ylab = "scaled vectorial capacity", main = "", debug = FALSE, # Shows detail if debugging. verbose = FALSE, # Shows even more detail if debugging. browse = FALSE) # Step through function one line at a time. { if(browse) browser(); begin.dir <- getwd() # Store Global working directory to reset at the end of function. setwd(directory) # Change working directory for saving. if(debug) {print("Error changing or saving directory")} if(pdf) pdf(file = file.name, width = plot.width, height = plot.height) if(!pdf & new.plot & quartz) quartz(width = plot.width, height = plot.height) if(debug) {print("Error initializing pdf")} ## Section caclulates Sigma_x, (e.g. scaled survival curve) surv.curv.0 <- sapply(age.range, s2day, start = 0, mort = mort, adulticide = 0, # Need to calculate survival curve with no control so can scale that to 1 and all other control scenarios relative to that. mu = mu, aa = aa, bb = bb, cc = cc, ss = ss) norm.scalar <- sum(surv.curv.0) surv.curv <- sapply(age.range, # Adult control incorporated in survival curves s2day, start = 0, mort = mort, adulticide = adulticide, mu = mu, aa = aa, bb = bb, cc = cc, ss = ss) surv.curv.larv <- surv.curv * (1 - larvicide) # Larval control reduces all age-classes by a proportion. surv.curv.norm <- surv.curv.larv / norm.scalar #Scales survival curve to be demographic distribution. if(debug) {print("Error calculating survival curve")} if(debug) {print("Just finished calculating survival curve reduced by larvicide")} ## Section calculates probability of living through the EIP given current age age.after.eip <- age.range + EIP # Calculate age they must survive to finish the EIP. p.surv.eip <- mapply(s2day, xx = age.after.eip, start = age.range, # Even though we want to start at age.range + 1, this is done in s2day() already. MoreArgs = list(mort = mort, adulticide = adulticide, mu = mu, aa = aa, bb = bb, cc = cc, ss = ss) ) if(debug) {print("Just finished calculating probability of living through EIP")} ## Section calculates e(x+n), life expectancy given surviving EIP life.exp.vec <- sapply(age.after.eip, exp.days, mort = mort, adulticide = adulticide, mu = mu, aa = aa, bb = bb, cc = cc, ss = ss) if(debug) {print("Just finished calculating life expectancies.")} ## Check that these three vectors are the same length error.check <- length(surv.curv.norm) != length(p.surv.eip) | length(surv.curv.norm) != length(life.exp.vec) if(error.check) { print(paste("survival curve vector is length", length(surv.curv.norm))) print(paste("probability of surviving EIP vector is length", length(p.surv.eip))) print(paste("life expectancy vector is length", length(life.exp.vec))) stop("These vectors cannot be different lengths") } ## Multiply probability of surviving EIP & life expectancy vector ## to get per individual contribution to vectorial capacity by age age.vc <- p.surv.eip * life.exp.vec ## Multiply these three vectors to get the age-class specific vectorial capacities class.vc <- surv.curv.norm * p.surv.eip * life.exp.vec if(class.spec) vec.to.plot <- class.vc if(!class.spec) vec.to.plot <- age.vc if(ylim.auto) ylim <- c(0,1.2 * max(vec.to.plot)) if(xlim.auto) xlim <- range(age.range) names.arg <- age.range # Label all the ages if(names.arg.set) # Or label only 5 of them to make it pretty { nums <- pretty(age.range,5) names.arg[-nums] <- NA names.arg[1] <- 0 } rep(NA, length(age.range))[pretty(age.range, 5)+1] if(new.plot) { plot(age.range, # Don't plot here so can use points(). vec.to.plot, type = "n", col = col, cex.lab = cex.lab, axes = FALSE, xlab = xlab, ylab = ylab, main = main, xlim = xlim, ylim = ylim) if(class.spec) y.ax.ticks <- seq(0,1, by = .5) if(!class.spec) y.ax.ticks <- seq(0,30, by = 10) x.ax.ticks <- pretty(age.range,5) if(label.yax) y.ax.labels <- y.ax.ticks if(!label.yax) y.ax.labels <- NA if(label.xax) x.ax.labels <- x.ax.ticks if(!label.xax) x.ax.labels <- NA if(verbose) print(y.ax.ticks) if(verbose) print(x.ax.ticks) if(verbose) print(y.ax.labels) if(verbose) print(x.ax.labels) axis(1, at = x.ax.ticks, labels = x.ax.labels) axis(2, at = y.ax.ticks, labels = y.ax.labels) } points(age.range, vec.to.plot, col = col, lty = lty, type = "l") if(polygon) { polygon(c(age.range, rev(age.range)), c(vec.to.plot, rep(0, length(age.range))), col = col, lty = lty) } if(debug) {print("Just finished with plot()")} if(legend.true) { legend(legend.loc, legend, col = col) if(debug) {print("Just finished with legend")} } if(pdf) dev.off() if(debug) {print("Just finished turning graphics off")} setwd(begin.dir) if(debug) {print("Just finished resetting directory")} } ## Univariate Sensitivity Analysis of VC to Control FUNCTIONS ######## ###################################################################### ###################################################################### ## Set working directory here graphics.off() ## The goal of this script is to plot sensitivity analyses of scaled ## vectorial capacity (C*, expected number of infectious biting-days) to ## increased mortality of adults (adult control) or decreased input of ## adults (larval control). ## Plot will be a curve showing how VC changes for different ## mosquito control parameters. The function can flexibly handle being ## an adult control sensitivity or larval control sensitivity with the ## "control" parameter. ## Should be able to plot absolute VC, and % of max VC (e.g. from no ## control scenario). Each of these has different implications, ## e.g. whether you are trying to reduce transmission below a certain ## threshold (look at absolute), or looking at effectiveness of control. univ.control.sens <- function(mort = "log", EIP = 2, control = "adult", adulticide.vec = seq(0, .5, length.out = 11), # Adult control vector for which to conduct sensitivity analysis. larvicide.vec = seq(0, 1, length.out = 11), # Larval control vector for which to conduct sensitivity analysis. mu = 1/32.6694, # Next five parameters describe mosquito mortality rates, see text for details. aa = .0018, bb = .1416, cc = 0, ss = 1.0730, age.range = 1:80, # Age range over which to calculate C* (past 80 contribution is negligible. VC.perc = TRUE, # Determines whether VC is plotted as a percentage or absolute value. pdf = FALSE, # Whether to produce a pdf figure inside the function. new.plot = TRUE, # Whether to plot figure on a new plotting device. quartz = TRUE, # Whether to plot figure through quartz (Mac only) plot.file.name = "VC cont plot.pdf", # Name to save pdf as if saving pdf. save.array = TRUE, # Whether to save the output array for future plotting/manipulation. array.file.name = "VC array", # Name of array file if saved. directory = "", # Where to save or load array from if saving or loading array. plot.width = 8, # Width of plot if using pdf() or quartz() in the function. plot.height = 6, # Height of plot if using pdf() or quartz() in the function. legend.true = FALSE, # Put legend in the plot? cex.lab = 1, # Legend scaling parameter. plot.only = FALSE, # Skips calculations and just plots a loaded array. plot.true = TRUE, # If FALSE, does not produce any graphics. array.to.load = "VC array", # Array to load if plot.only=TRUE. xlab.auto = TRUE, # Automatically set x-axis label. xlab = "reduced adult survival", ylab = "scaled vectorial capacity", ylim.auto = TRUE, # Automatically set y-axis limits. ylim = c(0, 10), xlim = c(0, .1), label.xax = TRUE, # Label x-axis? label.yax = TRUE, # Label y-axis? main = "", # Plot title type = "l", lwd = 2, lty = 1, col = "black", legend.lab = "", # Legend label debug = FALSE, # If debugging code, shows more. verbose = FALSE, # If debugging code, shows even more. browse = FALSE) # Step through function inside function. { if(browse) browser; begin.dir <- getwd() # Store Global working directory to reset at the end of function. setwd(directory) # Change working directory for saving. if(debug) print("Just finished changing or saving directory") ## Deal with flexibility of which type of control to do sensitivity to. if(control != "adult" & control != "larval") stop("Wrong specification of control, not 'adult' or 'larval'.") if(control == "adult") control.vec <- adulticide.vec if(control == "larval") control.vec <- larvicide.vec if(verbose) print(paste("control.vec is", round(control.vec,3))) if(!plot.only) { sys.time.begin <- proc.time() # Record beginning time to calculate simulation time at end. ## Initialize vector to fill with C* values and mosquito densities. vc.vec <- rep(NA, length(control.vec)) moz.dens <- rep(NA, length(control.vec)) ## Loop through control vector for(ii in 1:length(control.vec)) { if(control == "adult") # If looping through adulticide, set larvicide to 0. { larvicide <- 0 adulticide <- control.vec[ii] } if(control == "larval") # If looping through larvicide, set adulticide to 0. { adulticide <- 0 larvicide <- control.vec[ii] } if(debug) print("Just finished initializing loop over moz control parameters") if(verbose) print(paste("Looping at ", ii,"-th value of ", length(control.vec), control, "values.")) ## Section caclulates Sigma_x, (e.g. survival curve) surv.curv <- sapply(age.range, # Adult control incorporated in survival curves s2day, start = 0, mort = mort, adulticide = adulticide, mu = mu, aa = aa, bb = bb, cc = cc, ss = ss) if(debug) {print("Just finished calculating survival curve")} surv.curv.larv <- surv.curv * (1 - larvicide) # Larval control reduces all age-classes by a proportion. if(debug) print("Just finished calculating survival curve reduced by larvicide") ## Need to normalize mosquito population to size 1 in ## no control scenario and use that to compare other ## sizes to. This means 0 must be part of ## adulticide.vec & larvicide.vec if(control.vec[1] != 0) stop("First value of control parameter vector is not 0. Need that to be so to allow normalization of mosquito population size.") if(adulticide == 0 & larvicide == 0) norm.scalar <- sum(surv.curv.larv) surv.curv.norm <- surv.curv.larv / norm.scalar # Divide by scalar so that in control scenario total moz density = 1. if(debug) print("Just finished scaling mosquito densities") ## Section calculates probability of living through the EIP given current age age.after.eip <- age.range + EIP # Calculate age they must survive to finish the EIP. p.surv.eip <- mapply(s2day, xx = age.after.eip, start = age.range, # Even though we want to start at age.range + 1, this is done in s2day() already. MoreArgs = list(mort = mort, adulticide = adulticide, mu = mu, aa = aa, bb = bb, cc = cc, ss = ss) ) if(debug) print("Just finished calculating probability of living through EIP") ## Section calculates e(x+n), life expectancy given surviving EIP life.exp.vec <- sapply(age.after.eip, exp.days, mort = mort, adulticide = adulticide, mu = mu, aa = aa, bb = bb, cc = cc, ss = ss) if(debug) print("Just finished calculating life expectancies.") ## Check that these three vectors are the same length error.check <- length(surv.curv.norm) != length(p.surv.eip) | length(surv.curv.norm) != length(life.exp.vec) if(error.check) { print(paste("survival curve vector is length", length(surv.curv.norm))) print(paste("probability of surviving EIP vector is length", length(p.surv.eip))) print(paste("life expectancy vector is length", length(life.exp.vec))) stop("These vectors cannot be different lengths") } ## Multiply these three vectors to get the age-class specific vectorial capacities class.vc <- surv.curv.norm * p.surv.eip * life.exp.vec tot.vc <- sum(class.vc) vc.vec[ii] <- tot.vc moz.dens[ii] <- sum(surv.curv.norm) if(debug) print("Just finished putting tot.vc into vc.vec.") if(verbose) print(paste("C* on this iteration is", round(tot.vc,3))) } sys.time.end <- proc.time() processing.time <- sys.time.end - sys.time.begin # Calculate processing time. print(processing.time[3]) ## Save the array. control.vec.sav <- control.vec # Need to resave these as different names so can load later and not have if(save.array) save(vc.vec, moz.dens, control.vec.sav, processing.time, file = array.file.name) } if(plot.true) { ## Load an array if just plotting if(plot.only) load(array.to.load) ## Initialize graphics if(pdf) pdf(file = plot.file.name, width = plot.width, height = plot.height) if(!pdf & new.plot & quartz) quartz(width = plot.width, height = plot.height) if(debug) print("Just finished initializing pdf") if(verbose) print(vc.vec) ## Plot C* vs the control parameter. if(xlab.auto) # Set the xlab to match "control". { if(control == "adult") xlab <- "daily probability of death due to control" if(control == "larval") xlab <- "reduced recruitment of adults" } if(VC.perc) vc.vec <- vc.vec / vc.vec[1] * 100 if(ylim.auto) ylim <- c(0, 1.1*max(vc.vec)) if(VC.perc) ylim <- c(0, 100) if(new.plot) { plot(control.vec, vc.vec, xlim = xlim, ylim = ylim, xlab = xlab, ylab = ylab, main = main, type = "n", axes = FALSE, cex.lab = cex.lab) ## Axes labels y.ax.ticks <- pretty(ylim, 5) if(VC.perc) y.ax.ticks <- c(0,25,50,75,100) x.ax.ticks <- pretty(xlim, 5) y.ax.labels <- y.ax.ticks x.ax.labels <- x.ax.ticks if(!label.yax) y.ax.labels <- NA # Or don't label it if(!label.xax) x.ax.labels <- NA axis(1, at = x.ax.ticks, labels = x.ax.labels) axis(2, at = y.ax.ticks, labels = y.ax.labels) } points(control.vec, vc.vec, col = col, type = type, lwd = lwd, lty = lty) if(pdf) dev.off() } } ###################################################################### ###################################################################### ###################################################################### ###################################################################### #################### FIGURE MAKING CODE ############################## ###################################################################### ###################################################################### ###################################################################### ## Figure 1 ## Figure 1 ## 1A Mortality Curves ## 1B Life expectancy curves ## 1C Exponential age distribution ## 1D Logistic age distribution ## Must set the directory here ## Here we will set the survival curve parameters: ## This value of mu is used because it yields same life expectancy as ## the age-dependent model when mosquitoes age through discretized age ## classes. mu <- 1/32.6694 aa <- .0018 bb <- .1416 cc <- 0 ss <- 1.0730 age.range <- 1:80 adulticide <- 0 larvicide <- 0 ## Plot two panel plot with survival curves for log/exp mortality. setwd(directory) graphics.off() pdf("1 Figure.pdf", width = 6.89, height = 5.5) par(mfrow=c(2,2), mar = c(2, 5, 2, .5), oma = c(2, 0, 2, 0)) legend <- c("constant mortality", "age-dependent mortality") ## 1A - Mortality curves curve(mort.day(x, mort="log"), 0, 80, ylim = c(0, .2), cex.lab = cex.lab, ylab="daily probability of death", xlab = NA, main = "A", lwd = 2, lty = 1, axes = F) points(sapply(0:80, mort.day, mort="exp"), lty=2, lwd = 2, type="l") axis(side = 1, at = c(0, 20, 40, 60, 80), labels = FALSE, tick = TRUE) axis(side = 2, at = c(0, .1, .2), labels = TRUE, tick = TRUE) legend("top", legend=legend, lwd=2, lty=c(2,1), bty = "o", cex= .8) ## 1B - Life expectancy curves days <- 0:80 life.xp.log <- sapply(days, exp.days, mort="log") life.xp.exp <- sapply(days, exp.days, mort="exp") plot(days, life.xp.log, xlab = NA, ylab="expected days to live", main = "B", type="l", ylim = c(0, 40), lty=1, cex.lab = cex.lab, lwd=2, axes = F) points(days, life.xp.exp, type="l", lty=2, col="black", lwd=2) axis(side = 1, at = c(0, 20, 40, 60, 80), labels = FALSE, tick = TRUE) axis(side = 2, at = seq(0, 40, by = 10), labels = TRUE, tick = TRUE) ## 1C & 1D - Age distributions for(ii in c("exp","log")) { ylab <- "normalized density" main <- ifelse(ii=="exp", "C", "D") plot.surv(mort = ii, #Looping through mort type adulticide = adulticide, larvicide = larvicide, mu = mu, aa = aa, bb = bb, cc = cc, ss = ss, age.range = age.range, normalize = TRUE, pdf = FALSE, quartz = FALSE, new.plot = TRUE, directory = directory, plot.width = 8, plot.height = 6, legend.true = FALSE, legend.loc = "topright", legend = c("survival curve"), cex.lab = 1, col = "black", add = FALSE, polygon = TRUE, xlim.auto = TRUE, xlim = c(0,80), ylim = c(0,1), label.yax = TRUE, label.xax = TRUE, xlab = NA, ylab = ylab, main = main, debug = FALSE) } ## Outside axis labels. mtext("age in days", side = 1, line = .5, outer = TRUE, adj = .28, cex = .8) mtext("age in days", side = 1, line = .5, outer = TRUE, adj = .84, cex = .8) ## Finish PDF dev.off() ## Figure 2 ## Shows the individual level and age-class level ## contributions to C* (scaled vectorial capacity). ## 2A - Individual-level C* for EIP = 2, 11, EXP. ## 2B - Individual-level C* for EIP = 2, 11, LOG. ## 2C - Class-level C* for EIP = 2, 11, EXP. ## 2D - Class-level C* for EIP = 2, 11, LOG. ## Must set the directory here ## Here we will set the survival curve parameters: mu <- 1/32.6694 aa <- .0018 bb <- .1416 cc <- 0 ss <- 1.0730 age.range <- 1:80 adulticide <- 0 larvicide <- 0 directory <- "/figures" setwd(directory) ## Set graphics parameters graphics.off() pdf("2 Figure.pdf", width = 6.89, height = 5.5) par(mfrow=c(2,2), mar = c(2, 5, 2, .5), oma = c(2, 0, 2, 0)) ## Set looping parameters class.spec.vec <- c(FALSE, TRUE) ylab.vec <- c("per-individual vectorial capacity", "age-class vectorial capacity") xlab.vec <- c(NA, "age in days") mort.vec <- c("exp","log") EIP.vec <- c(3, 11) ymax.vec <- c(30, 1) col.vec <- c("blue","green") label.yax.vec <- c(TRUE, FALSE) label.xax.vec <- c(FALSE, TRUE) main.vec <- t(matrix(c("A", "B", "C","D"), 2, 2)) ## Plotting loop for(ii in 1:length(class.spec.vec)) { for(jj in 1:length(mort.vec)) { for(kk in 1:length(EIP.vec)) { if(jj==2) ylab <- NA if(jj==1) ylab <- ylab.vec[ii] new.plot <- ifelse(kk==1, TRUE, FALSE) # Makes second EIP plot on same figure as first plot.vc(mort = mort.vec[jj], # LOOPING through this par EIP = EIP.vec[kk], adulticide = 0, larvicide = 0, mu = mu, aa = aa, bb = bb, cc = cc, ss = ss, age.range = age.range, pdf = FALSE, quartz = FALSE, new.plot = new.plot, file.name = "age-spec vc curve.pdf", directory = directory, plot.width = 8, plot.height = 6, class.spec = class.spec.vec[ii], # If true plots VC by whole age class, if false plots per-indiv VC legend.true = FALSE, legend.loc = "topright", legend = c("vectorial capacity curve"), cex.lab = 1, col = col.vec[kk], add = FALSE, polygon = TRUE, lty = 1, xlim = c(0,80), ylim.auto = FALSE, xlim.auto = TRUE, label.yax = label.yax.vec[jj], label.xax = label.xax.vec[ii], ylim = c(0, ymax.vec[ii]), xlab = xlab.vec[ii], ylab = ylab, main = main.vec[ii, jj], debug = FALSE, verbose = FALSE, browse = FALSE) if(ii == 1 & jj ==2) legend("top", paste("EIP = ", rev(EIP.vec)), col = rev(col.vec), lty = 1, cex = .8) } } } ## Add labels. mtext("age in days", side = 1, line = .5, outer = TRUE, adj = .28, cex = .8) mtext("age in days", side = 1, line = .5, outer = TRUE, adj = .82, cex = .8) ## Finish PDF. dev.off() ## Figure 3 ## Shows affect of adult control on several plots. ## MORT = EXP ## Figure 3A - Age distribution ## Figure 3B - Life expectancy ## Figure 3C - DENG C* ## Figure 3D - CHIKV C* ## MORT = LOG ## Figure 3E - Age distribution ## Figure 3F - Life expectancy ## Figure 3G - DENG C* ## Figure 3H - CHIKV C* ## Here we will set the survival curve parameters: mu <- 1/32.6694 aa <- .0018 bb <- .1416 cc <- 0 ss <- 1.0730 age.range <- 1:80 adulticide <- 0 larvicide <- 0 ## Start up graphics graphics.off() pdf("3 Figure.pdf", width = 4.92, height = 7) par(mfrow=c(4,2), mar = c(2, 5, 2, .5), oma = c(2, 0, 2, 0)) ## Set up looping parameters mort.vec <- c("exp","log") adulticide <- .05 EIP.vec <- c(11, 3) label.yax.vec <- c(TRUE, FALSE) ## Figure 3A-B - Age distributions for(ii in 1:length(mort.vec)) { ylab <- "normalized density" if(ii>1) ylab <- "" main <- ifelse(ii == 1, "A", "B") plot.surv(mort = mort.vec[ii], # Plot age distribution (survival curve) for no control adulticide = 0, larvicide = 0, mu = mu, aa = aa, bb = bb, cc = cc, ss = ss, age.range = age.range, normalize = TRUE, pdf = FALSE, quartz = FALSE, new.plot = TRUE, file.name = "survival curve.pdf", directory = directory, plot.width = 8, plot.height = 6, legend.true = FALSE, legend.loc = "topright", legend = c("survival curve"), cex.lab = 1, col = "black", add = FALSE, xlim = c(0,80), xlim.auto = TRUE, label.yax = label.yax.vec[ii], label.xax = FALSE, ylim = c(0,1), names.arg.set = TRUE, xlab = "", ylab = ylab, main = main, debug = FALSE) plot.surv(mort = mort.vec[ii], adulticide = adulticide, larvicide = 0, mu = mu, aa = aa, bb = bb, cc = cc, ss = ss, age.range = age.range, normalize = TRUE, pdf = FALSE, quartz = FALSE, new.plot = FALSE, # Add this to the last plot. file.name = "survival curve.pdf", directory = "/Users/sbellan/Documents/R files/MPH thesis/PLOSOne/figures/", plot.width = 8, plot.height = 6, legend.true = FALSE, legend.loc = "topright", legend = c("survival curve"), cex.lab = 1, col = "red", add = TRUE, xlim = c(0,80), xlim.auto = TRUE, names.arg.set = TRUE, ylim = c(0,1), xlab = "", ylab = "", label.yax = FALSE, label.xax = FALSE, main = main, debug = FALSE) if(mort.vec[ii] == "exp") load(file="norm.scal.surv.exp") if(mort.vec[ii] == "log") load(file="norm.scal.surv.log") y.chik <- s2day(5, mort = mort.vec[ii])/sum.surv.curv y.deng <- s2day(12, mort = mort.vec[ii])/sum.surv.curv lines(c(5, 5), c(0, y.chik), lty = 2, col = "blue") lines(c(12, 12), c(0, y.deng), lty = 2, col = "green") } ## End Figure 3A-B ################################################### ###################################################################### ## Figure 3C-D - Life expectancy plots for(ii in 1:length(mort.vec)) { main <- ifelse(ii == 1, "C", "D") ylab <- "expected days to live" if(ii>1) ylab <- "" life.xp.no.cont <- sapply(age.range, exp.days, mort = mort.vec[ii], adulticide = 0, # Without control mu = mu, aa = aa, bb = bb, cc = cc, ss = ss) life.xp <- sapply(age.range, exp.days, mort = mort.vec[ii], adulticide = adulticide, # With control mu = mu, aa = aa, bb = bb, cc = cc, ss = ss) ylim <- c(0,1.2*max(life.xp.no.cont)) plot(age.range, life.xp.no.cont, xlab = "", ylab = ylab, main = main, type = "l", axes = FALSE, col = "black", lty = 1, ylim = ylim, cex.lab = cex.lab, lwd = 2) if(ii == 1) y.ticks <- c(0,10,20,30) if(ii > 1) y.ticks <- NA axis(1, at = pretty(age.range, 5), labels = NA) axis(2, at = c(0,10,20, 30), labels = y.ticks) points(age.range, life.xp, lty = 1, col = "red", type = "l", lwd = 2) y.chik <- exp.days(5, mort = mort.vec[ii]) y.deng <- exp.days(12, mort = mort.vec[ii]) lines(c(5, 5), c(0, y.chik), lty = 2, col = "blue") lines(c(12, 12), c(0, y.deng), lty = 2, col = "green") } ## End of Figure 3C-D ################################################ ###################################################################### ## Figure 3E-H - C* curves label.xax.vec <- c(FALSE, TRUE) label.yax.vec <- c(TRUE, FALSE, FALSE, FALSE) main.vec <- t(matrix(c("E","F","G","H"), 2, 2)) for(jj in 1:length(EIP.vec)) { for(ii in 1:length(mort.vec)) { main <- main.vec[jj, ii] ylab <- "C*" if(ii>1) ylab <- "" xlab <- "" if(jj>1) xlab <- "age in days" plot.vc(mort = mort.vec[ii], # LOOP EIP = EIP.vec[jj], # LOOP adulticide = 0, # NO LOOP because no control larvicide = 0, # NO LOOP because no control mu = mu, aa = aa, bb = bb, cc = cc, ss = ss, age.range = age.range, pdf = FALSE, quartz = FALSE, new.plot = TRUE, # Initialize plot here file.name = "age-spec vc curve.pdf", directory = "/Users/sbellan/Documents/R files/MPH thesis/PLOSOne/figures/", plot.width = 8, plot.height = 6, class.spec = TRUE, # If true plots VC by whole age class, if false plots per-indiv VC legend.true = FALSE, legend.loc = "topright", legend = c("vectorial capacity curve"), cex.lab = 1, col = "black", # Plot no control in black. add = FALSE, polygon = TRUE, xlim = c(0,100), ylim.auto = FALSE, xlim.auto = TRUE, label.xax = label.xax.vec[jj], label.yax = label.yax.vec[ii], names.arg.set = TRUE, ylim = c(0,1), xlab = xlab, ylab = ylab, main = main, debug = FALSE, browse = FALSE) plot.vc(mort = mort.vec[ii], # LOOP EIP = EIP.vec[jj], # LOOP adulticide = adulticide, larvicide = 0, mu = mu, aa = aa, bb = bb, cc = cc, ss = ss, age.range = age.range, pdf = FALSE, quartz = TRUE, new.plot = FALSE, # Add to no control plot file.name = "age-spec vc curve.pdf", directory = "/Users/sbellan/Documents/R files/MPH thesis/PLOSOne/figures/", plot.width = 8, plot.height = 6, class.spec = TRUE, # If true plots VC by whole age class, if false plots per-indiv VC legend.true = FALSE, legend.loc = "topright", legend = c("vectorial capacity curve"), cex.lab = 1, col = "red", # Plot controlled population in red. add = FALSE, polygon = TRUE, xlim = c(0,100), ylim.auto = TRUE, xlim.auto = TRUE, label.xax = FALSE, label.yax = FALSE, names.arg.set = TRUE, ylim = c(0,3), xlab = "", ylab = "scaled vectorial capacity", main = "", debug = FALSE, browse = FALSE) } } ## Text to add at end mtext("age in days", side = 1, line = .5, outer = TRUE, adj = .28, cex = .8) mtext("age in days", side = 1, line = .5, outer = TRUE, adj = .84, cex = .8) ## Finish pdf dev.off() ## Figure 3 ## Shows affect of larval control on several plots. ## MORT = EXP ## Figure 3A - Age distribution ## Figure 3B - Life expectancy ## Figure 3C - DENG C* ## Figure 3D - CHIKV C* ## MORT = LOG ## Figure 3E - Age distribution ## Figure 3F - Life expectancy ## Figure 3G - DENG C* ## Figure 3H - CHIKV C* ## Here we will set the survival curve parameters: mu <- 1/32.6694 aa <- .0018 bb <- .1416 cc <- 0 ss <- 1.0730 age.range <- 1:80 adulticide <- 0 larvicide <- 0 ## Start up graphics graphics.off() pdf("4 Figure.pdf", width = 4.92, height = 7) par(mfrow=c(4,2), mar = c(2, 5, 2, .5), oma = c(2, 0, 2, 0)) ## Set up looping parameters mort.vec <- c("exp","log") larvicide <- .5 EIP.vec <- c(11, 3) label.yax.vec <- c(TRUE, FALSE) ## Figure 3A-B - Age distributions for(ii in 1:length(mort.vec)) { ylab <- "normalized density" if(ii>1) ylab <- "" main <- ifelse(ii == 1, "A", "B") plot.surv(mort = mort.vec[ii], # Plot age distribution (survival curve) for no control adulticide = 0, larvicide = 0, mu = mu, aa = aa, bb = bb, cc = cc, ss = ss, age.range = age.range, normalize = TRUE, pdf = FALSE, quartz = FALSE, new.plot = TRUE, file.name = "survival curve.pdf", directory = directory, plot.width = 8, plot.height = 6, legend.true = FALSE, legend.loc = "topright", legend = c("survival curve"), cex.lab = 1, col = "black", add = FALSE, xlim = c(0,80), xlim.auto = TRUE, label.yax = label.yax.vec[ii], label.xax = FALSE, ylim = c(0,1), names.arg.set = TRUE, xlab = "", ylab = ylab, main = main, debug = FALSE) plot.surv(mort = mort.vec[ii], adulticide = 0, larvicide = larvicide, mu = mu, aa = aa, bb = bb, cc = cc, ss = ss, age.range = age.range, normalize = TRUE, pdf = FALSE, quartz = FALSE, new.plot = FALSE, # Add this to the last plot. file.name = "survival curve.pdf", directory = "/Users/sbellan/Documents/R files/MPH thesis/PLOSOne/figures/", plot.width = 8, plot.height = 6, legend.true = FALSE, legend.loc = "topright", legend = c("survival curve"), cex.lab = 1, col = "red", add = TRUE, xlim = c(0,80), xlim.auto = TRUE, names.arg.set = TRUE, ylim = c(0,1), xlab = "", ylab = "", label.yax = FALSE, label.xax = FALSE, main = main, debug = FALSE) if(mort.vec[ii] == "exp") load(file="norm.scal.surv.exp") if(mort.vec[ii] == "log") load(file="norm.scal.surv.log") y.chik <- s2day(5, mort = mort.vec[ii])/sum.surv.curv y.deng <- s2day(12, mort = mort.vec[ii])/sum.surv.curv lines(c(5, 5), c(0, y.chik), lty = 2, col = "blue") lines(c(12, 12), c(0, y.deng), lty = 2, col = "green") } ## End Figure 3A-B ################################################### ###################################################################### ## Figure 3C-D - Life expectancy plots for(ii in 1:length(mort.vec)) { main <- ifelse(ii == 1, "C", "D") ylab <- "expected days to live" if(ii>1) ylab <- "" life.xp.no.cont <- sapply(age.range, exp.days, mort = mort.vec[ii], adulticide = 0, # Without control mu = mu, aa = aa, bb = bb, cc = cc, ss = ss) life.xp <- sapply(age.range, exp.days, mort = mort.vec[ii], adulticide = 0, # With control mu = mu, aa = aa, bb = bb, cc = cc, ss = ss) ylim <- c(0,1.2*max(life.xp.no.cont)) plot(age.range, life.xp.no.cont, xlab = "", ylab = ylab, main = main, type = "l", axes = FALSE, col = "black", lty = 1, ylim = ylim, cex.lab = cex.lab, lwd = 2) if(ii == 1) y.ticks <- c(0,10,20,30) if(ii > 1) y.ticks <- NA axis(1, at = pretty(age.range, 5), labels = NA) axis(2, at = c(0,10,20, 30), labels = y.ticks) points(age.range, life.xp, lty = 2, col = "red", type = "l", lwd = 2) y.chik <- exp.days(5, mort = mort.vec[ii]) y.deng <- exp.days(12, mort = mort.vec[ii]) lines(c(5, 5), c(0, y.chik), lty = 2, col = "blue") lines(c(12, 12), c(0, y.deng), lty = 2, col = "green") } ## End of Figure 3C-D ################################################ ###################################################################### ## Figure 3E-H - C* curves label.xax.vec <- c(FALSE, TRUE) label.yax.vec <- c(TRUE, FALSE, FALSE, FALSE) main.vec <- t(matrix(c("E","F","G","H"), 2, 2)) for(jj in 1:length(EIP.vec)) { for(ii in 1:length(mort.vec)) { main <- main.vec[jj, ii] ylab <- "C*" if(ii>1) ylab <- "" xlab <- "" if(jj>1) xlab <- "age in days" plot.vc(mort = mort.vec[ii], # LOOP EIP = EIP.vec[jj], # LOOP adulticide = 0, # NO LOOP because no control larvicide = 0, # NO LOOP because no control mu = mu, aa = aa, bb = bb, cc = cc, ss = ss, age.range = age.range, pdf = FALSE, quartz = FALSE, new.plot = TRUE, # Initialize plot here file.name = "age-spec vc curve.pdf", directory = "/Users/sbellan/Documents/R files/MPH thesis/PLOSOne/figures/", plot.width = 8, plot.height = 6, class.spec = TRUE, # If true plots VC by whole age class, if false plots per-indiv VC legend.true = FALSE, legend.loc = "topright", legend = c("vectorial capacity curve"), cex.lab = 1, col = "black", # Plot no control in black. add = FALSE, polygon = TRUE, xlim = c(0,100), ylim.auto = FALSE, xlim.auto = TRUE, label.xax = label.xax.vec[jj], label.yax = label.yax.vec[ii], names.arg.set = TRUE, ylim = c(0,1), xlab = xlab, ylab = ylab, main = main, debug = FALSE, browse = FALSE) plot.vc(mort = mort.vec[ii], # LOOP EIP = EIP.vec[jj], # LOOP adulticide = 0, larvicide = larvicide, mu = mu, aa = aa, bb = bb, cc = cc, ss = ss, age.range = age.range, pdf = FALSE, quartz = TRUE, new.plot = FALSE, # Add to no control plot file.name = "age-spec vc curve.pdf", directory = "/Users/sbellan/Documents/R files/MPH thesis/PLOSOne/figures/", plot.width = 8, plot.height = 6, class.spec = TRUE, # If true plots VC by whole age class, if false plots per-indiv VC legend.true = FALSE, legend.loc = "topright", legend = c("vectorial capacity curve"), cex.lab = 1, col = "red", # Plot controlled population in red. add = FALSE, polygon = TRUE, xlim = c(0,100), ylim.auto = TRUE, xlim.auto = TRUE, label.xax = FALSE, label.yax = FALSE, names.arg.set = TRUE, ylim = c(0,3), xlab = "", ylab = "scaled vectorial capacity", main = "", debug = FALSE, browse = FALSE) } } ## Text to add at end mtext("age in days", side = 1, line = .5, outer = TRUE, adj = .28, cex = .8) mtext("age in days", side = 1, line = .5, outer = TRUE, adj = .84, cex = .8) ## Finish pdf dev.off() ## Figure 5 - 1 panel ## Shows how C* % reductions depend on mortality assumption. ## Set working directory here. ## Here we will set the survival curve parameters: mu <- 1/32.6694 aa <- .0018 bb <- .1416 cc <- 0 ss <- 1.0730 age.range <- 1:80 ## Start up graphics graphics.off() pdf("5 Figure.pdf", width = 4.92, height = 4) legend <- c("EIP = 11", "EIP = 3","constant mortality", "age-dependent mortality") ## Set up sensitivity range vectors length.out <- 51 ## The exp term below is to make Figure 6 have more equally spaced points because it is plotted on a log scale kinda. adulticide.vec <- c(0, exp(seq(-18,-7,by=.25)), seq(.1/length.out, .1, length.out = length.out),.4,.6,.8,1) save(adulticide.vec, file = "adulticide.vec") ## Set up looping vectors. mort.vec <- c("exp", "log") control.type <- c("adult", "larval") EIP.vec <- c(11, 3) col.vec <- c("green", "blue") ## Make names to save arrays array.name.vec <- c("C-11", "C-3", "AD-11", "AD-3") save(array.name.vec, file = "Fig5 array.name.vec") array.ii <- 1 for(ii in 1:length(mort.vec)) { for(kk in 1:length(EIP.vec)) { new.plot <- FALSE if(array.ii == 1) new.plot <- TRUE array.name <- array.name.vec[array.ii] univ.control.sens(mort = mort.vec[ii], # LOOP EIP = EIP.vec[kk], # LOOP control = "adult", adulticide.vec = adulticide.vec, larvicide.vec = larvicide.vec, mu = mu, aa = aa, bb = bb, cc = cc, ss = ss, age.range = age.range, VC.perc = TRUE, # LOOP pdf = FALSE, new.plot = new.plot, # LOOP quartz = FALSE, plot.file.name = "VC cont plot.pdf", save.array = TRUE, array.file.name = array.name, directory = "/Users/sbellan/Documents/R files/MPH thesis/PLOSOne/figures/", plot.width = 8, plot.height = 6, legend.true = FALSE, cex.lab = 1, plot.only = TRUE, # Skips calculating and just loads one to plot. plot.true = TRUE, # Plots C* vs control array.to.load = array.name, xlab.auto = FALSE, xlab = "daily probability of adult death due to control", ylab = "% of total C*", ylim.auto = FALSE, ylim = c(0, 30), xlim = c(0, .1), label.xax = TRUE, label.yax = TRUE, main = "", type = "l", lwd = 2, lty = ii, # Code lines by EIP. col = col.vec[kk], legend.lab = "", debug = FALSE, verbose = FALSE, browse = FALSE) array.ii <- array.ii + 1 } } legend("topright", legend, lty = c(1,1,1,2), lwd = 2, col = c("green", "blue", "black", "black"), bty = "o", cex = .8) dev.off() ## Figure 6 ## This figure will display the % extra adult control predicted to ## need to control a disease in the age-dependent model vs the ## constant model vs the % reduction in VC desired to be achieved. ## It calls the arrays that are calculated in Figure 5. ## Set working directory here. pdf("6 Figure.pdf", width = 4.92, height = 4) ## Load arrays load("Fig5 array.name.vec") # Loads array.name.vec load("adulticide.vec") span <- .25 # LOESS parameter newdata <- 0:100 ## Below code produces vectors that are % of total C* for each ## simulation and then uses LOESS regressions to predict the value of ## adulticide necessary to reduce C* by a % given in the newdata ## vector. That way comparisons can be made between simulations for ## the value of control needed to reduce C* to a certain level. load(array.name.vec[1]) # C-11 vc.vec.perc <- vc.vec/vc.vec[1]*100 loess.c11 <- loess(adulticide.vec ~ vc.vec.perc, span = span) pred.c11 <- predict(loess.c11, newdata = newdata) load(array.name.vec[2]) # C-3 vc.vec.perc <- vc.vec/vc.vec[1]*100 loess.c3 <- loess(adulticide.vec ~ vc.vec.perc, span = span) pred.c3 <- predict(loess.c3, newdata = newdata) load(array.name.vec[3]) # AD-11 vc.vec.perc <- vc.vec/vc.vec[1]*100 loess.ad11 <- loess(adulticide.vec ~ vc.vec.perc, span = span) pred.ad11 <- predict(loess.ad11, newdata = newdata) load(array.name.vec[4]) # AD-3 vc.vec.perc <- vc.vec/vc.vec[1]*100 loess.ad3 <- loess(adulticide.vec ~ vc.vec.perc, span = span) pred.ad3 <- predict(loess.ad3, newdata = newdata) ## Calculate % extra adulticide needed to get C* to value in newdata diff11 <- (pred.ad11 - pred.c11)/pred.c11 * 100 diff3 <- (pred.ad3 - pred.c3)/pred.c3 * 100 plot(newdata, diff11, xlim = c(100, 20), ylim = c(0, 120), axes = F, col = "green", lwd = 2, lty = 1, xlab = "% reduction in C*", ylab = "% extra control effort in age-dependent model", type = "l") points(newdata, diff3, col = "blue", lwd = 2, lty = 1, type = "l") ats.x <- seq(100,20,by=-20) ats.y <- seq(0,120,by=20) axis(1, at = ats.x, label = ats.x) axis(2, at = ats.y, label = ats.y) legend("topright", c("EIP = 11", "EIP = 3"), col = c("green","blue"), lty = 1, lwd = 2, bty = "o", cex = .8) dev.off()