In ian-hastings/rotations: Simulating the Evolution of Insecticide Resistance under Rotations

An rmarkdown document that can be modified and knit to run scenarios within a specified experiment. Output files are saved to the specified outFolder.

library(rotations)
library(ggplot2)
library(gridExtra)
library(readxl) # for reading the excel file that allows me to keep formatting to see inputs clearly

outFolder <- "C:\\Dropbox\\resistanceResults\\rotations\\"

nscenarios <- 10000 #200
plotyn <- 0 #allow turning plots on and off from here, may want to turn plots off if doing large number of runs

# now specify an experiment from the multi-experiment file
#id_expt <- 'A1'
#id_expt <- 'A2'
#id_expt <- 'A3'
#id_expt <- 'A4'
id_expt <- 'B1'
#id_expt <- 'B2'
#id_expt <- 'B3'
#id_expt <- 'B4'
#id_expt <- 'C1'
#id_expt <- 'C4'
#id_expt <- 'D4'

# to move the pdf of a preceeding run of this file to a saved location
# RUN the line above and below after knitting the file
# ASSUMING id_expt still set to previous run and still same day
# BEWARE that id_expt will not be set correctly in home env if file was knitted
#file.rename("rotation_runs2019.pdf",paste0(outFolder,format(Sys.time(), "%Y-%m-%d-"),id_expt,"-results-figs.pdf"))


# old ways of doing it
# infile_multi <- 'rotations-experiments-creating-201909.xlsx'
# infile_multi <- system.file("extdata", infile_multi, package = "rotations", mustWork = TRUE)
# inex_all <- read_excel(infile_multi)

# older way with individual experiment files
# to add new files currently need to reinstall package because infiles are in extdata
# infile <- "expt-A2-freq-insame-costs-nodisp.csv"
# infile <- "expt-A4-freq-insame-costs-disp.csv"
# # get the id for the expot to include in outfilenames, specifies experiment attributes
# id_expt <- substr(infile,6,7)
# # read an experiment file specifying ranges
# inex <- read_in_expt(infile, nscenarios=nscenarios) 


inex <- read_in_expt(nscenarios = nscenarios, multi=TRUE, id_expt=id_expt)

# text to go at start of document
cat("\nThere are 16 potential experiments :\n\n    ")

cat("Experiment :",id_expt," scenario number is pagenum -1 \n\n    ")

cat("Experiment key :\n    ")
cat("A: resistance frequency threshold, insecticides the same\n    ")
cat("B: resistance frequency threshold, different insecticides\n    ")
cat("C: mortality threshold, insecticides the same\n    ")
cat("D: mortality threshold, different insecticides\n    ")

cat("1: no costs, no dispersal\n    ")
cat("2: costs, no dispersal\n    ")
cat("3: no costs, dispersal\n    ")
cat("4: costs, dispersal\n    ")

cat("\nNumbers on the right are the number of generations that this insecticide is in use and the corresponding resistance is below threshold.\n   ")

cat("\nOn each page a Sequence is at the top and an annual Rotation at the bottom.")

# create an object containing inputs for all scenarios
linmulti <- set_run_inputs(inex=inex)

# run all scenarios
for(i in 1:nscenarios)
{
  #message should output to the R markdown console
  if (i%%200 == 1) message("scenario ",i," of ",nscenarios," ",Sys.time())

  linputs <- get_one_in(linmulti, scen_num=i)
  # creating 2nd version of inputs that retains multi-insectcide strings so that can be saved in an output file as a single row
  linputs_tosave <- get_one_in(linmulti, scen_num=i, convert_diff_insecticide_strings = FALSE)

  if (linputs$rot_interval == 0) stop("in rotation_runs2019 rot_interval in input file needs to be > 0 to allow it to be compared to a sequence")

  linputs$plot <- plotyn #allows setting whether to plot from script
  linputs$add_gens_under50 <- TRUE #to add text to plots of gens under threshold

  #record rotation rot interval locally otherwise it gets overwritten by 0
  rot_rot_interval <- linputs$rot_interval

  #to add label between plots in output doc
  #this adds label underneath, not my intention but does the job
  if (plotyn) cat("\n\n\\pagebreak\n    ") #needs results='asis' in the chunk options
  #cat("scenario",i,"\n     ") stopped plots from aligning
  #instead just remember that scenario number is pagenum -1 


  # to compare sequence to rotation run all the scenarios above for both seq & rot 
  # need to do the rotation 2nd so that rot_interval from the input file can be saved with the outputs
  for(rot_or_not in 0:1)
  {  
    #to get rot_int of 0 in alternate runs - the other run uses the value from the input file
    if (rot_or_not == 0) linputs$rot_interval <- 0
    else linputs$rot_interval <- rot_rot_interval


    ##################
    # run one scenario
    dfres <- do.call(run_rot, linputs)

    # summarise results per scenario
    # this calcs for whole run, at end of run_rot() it calcs per insecticide
    if ( linputs$mort_or_freq == 'freq' ) threshold <- linputs$freq_thresh
    else if ( linputs$mort_or_freq == 'mort' ) threshold <- linputs$mort_thresh 

    gens_u <- gens_under_thresh(dfres=dfres, threshold=threshold, mort_or_freq=linputs$mort_or_freq)

    # 10/12/19 change to recording end generations and mean mortality
    gens <- calc_max_gens(dfres=dfres)
    cum_mort <- cumulat_mortality(dfres=dfres)

    # to be able to save as for rot or seq
    if ( linputs$rot_interval == 0 ) {gens_seq <- as.numeric(gens)
                                      gens_u_seq <- as.numeric(gens_u)
                                      mort_seq <- cum_mort/gens_seq}
    else                             {gens_rot <- as.numeric(gens)
                                      gens_u_rot <- as.numeric(gens_u)
                                      mort_rot <- cum_mort/gens_rot}

    #code to rerun approx this scenario, can be copied from pdf and pasted into R
    #has to cope with converting vectors to strings
    #spaces between args stop the pdf from messing with words
    if (plotyn) cat(paste0( "run_rot(rot_interval=",linputs$rot_interval,", ", 
                "n_insecticides=",linputs$n_insecticides,", ",
                #"mort_or_freq='",linputs$mort_or_freq,"', ", #quotes dont work, freq is now default
                #'mort_or_freq="',linputs$mort_or_freq,'", ',
                "start_freqs=",signif(linputs$start_freqs,2),", ", 
                "cost=c(",paste(signif(linputs$cost,2),collapse=','),"), ", 
                "eff=c(",paste(signif(linputs$eff,2),collapse=','),"), ",               
                "dom_sel=c(",paste(signif(linputs$dom_sel,2),collapse=','),"), ",
                "dom_cos=c(",paste(signif(linputs$dom_cos,2),collapse=','),"), ",
                "rr=c(",paste(signif(linputs$rr,2),collapse=','),"), ",
                "expo_hi=",signif(linputs$expo_hi,2),", ",
                "male_expo_prop=",signif(linputs$male_expo_prop,2),", ",
                "coverage=",signif(linputs$coverage,2),", ",
                "migration=",signif(linputs$migration,2),", ",
                "max_gen=",linputs$max_gen,", ", 
                "change_interval=",linputs$change_interval,
                #"min_rwr_interval=",linputs$min_rwr_interval,  
                ")\n\n") )

  } #end rot_or_not loop

  # add scenario id column and difference onto results
  # df_in_out is one row

  #2020-02-25 beware previously when insecticides are different df_in_out ended up with multiple rows with results repeated
  #which got unintentionally semi-corrected below
  #df_in_out <- as.data.frame(c(linputs, gens_rot=gens_rot, gens_seq=gens_seq, 
  #now using linputs_tosave which contains unchanged strings, and had to add stringsAsFactors = FALSE
  df_in_out <- as.data.frame(c(linputs_tosave, gens_rot=gens_rot, gens_seq=gens_seq,                                
                                        gens_u_rot=gens_u_rot, gens_u_seq=gens_u_seq,
                                        mort_rot=mort_rot, mort_seq=mort_seq, id=i), stringsAsFactors = FALSE)


  df_in_out$gens_rot_minus_seq <- df_in_out$gens_rot - df_in_out$gens_seq
  df_in_out$gens_u_rot_minus_seq <- df_in_out$gens_u_rot - df_in_out$gens_u_seq  
  df_in_out$mort_rot_minus_seq <- df_in_out$mort_rot - df_in_out$mort_seq

  # if the first scenario create df to store results of all scenarios
  if (i==1)
  {
    df_res_all <- data.frame(matrix(NA, ncol=length(names(df_in_out)), nrow=nscenarios))
    names(df_res_all) <- names(df_in_out)
  }

  #2020-02-25 beware (see above) previoulsy when insecticides are different this saved just first row of df_in_out
  #result is that correct results are saved but with just the inputs for the first insecticide
  #this shouldn't effect the output results
  #now corrected above by keeping the values for multiple different insecticides as a string
  #putting results into one row of overall results dataframe
  df_res_all[i,] <- df_in_out


  #can save partial outputs in case doesn't make to end
  if (i==nscenarios | i%%1000 == 0)
  {

    # save object containing inputs and results as rda for analysis
    # by default saved to inst/documents folder
    save(df_res_all, file=paste0(outFolder,format(Sys.time(), "%Y-%m-%d-%H-%M-"),'df_res_all_rot',id_expt,nscenarios,'.rda')) #paste0(outFolder,'*.rda')) 
    # also save the experiment input file
    save(inex, file=paste0(outFolder,format(Sys.time(), "%Y-%m-%d-%H-%M-"),'inex-',id_expt,'.rda'))   
  } 

} #end scenario loop  

ian-hastings/rotations documentation built on Dec. 14, 2020, 11:42 p.m.