library(rotations)
  library(ggplot2)
  library(gridExtra)

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

10,000 runs comparing rotations and sequences.

#eval=FALSE just to temporarily hide this text

# Starting question : is there any benefit to rotations when there are no costs and no migration to an untreated area ?
# 
# Test 1. No refugia (i.e. Coverage=100%) or fitness costs. Would expect all policies to be equivalent.
# 
# Test 2. Dominance of resistance and fitness costs both set to 0.5.
# I suspect one of the advantages of rotations (and the presence of refugia) is that they maintain lower resistance allele frequencies which slows the spread of recessive resistance alleles. I suspect/guess/intuit that if we make them semi-dominate these effects should disappear and all policies should be the same. I’m probably wrong but it’s a nice test…..
# 

#\pagebreak
#results='hide',
#comment=NA removes ## before text output

#start small & build it up
#maybe create an equiv. of resistance::sensiAnPaperPart.r later

# example for 1 run
# linputs <- list(max_gen = 50, n_insecticides = 2)
# do.call(run_rot, linputs)
# dfinputs <- as.data.frame(linputs)

# So one way of doing is to create an linputs list for each run & pass to run_rot()
# so run_rot() only ever runs for 1 scenario
# in contrast resistance::runModel2() accepted an input matrix & would run for multiple scenarios
# not sure which is better

# set random seed so that results can be reproduced
#set.seed(1)
set.seed(2)
max_gen         <- 500
plot_scenarios  <- FALSE
n_scenarios     <- 10 #000 
already_ran     <- FALSE
#already_ran     <- TRUE

# if already done runs, just load data for plots
if (already_ran)
{
  #C:/rsprojects/rotations/inst/documents/

  #load("df_res_all_rot10000.rda")
  load("df_res_all_rot_mig1000.rda")

}else
{
  df_res_all <- data.frame()

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

    # inputs (could also be saved to a dataframe first & extracted one row at a time)
    linputs <- list()

    ## constant inputs ##
    linputs$coverage <- 1
    linputs$migration <- 0
    #linputs$cost <- 0 #variable cost below
    linputs$max_gen <- max_gen

    #beware min_rwr_interval, it's to stop repeated small rotations
    linputs$min_rwr_interval <- 5
    linputs$min_gens_switch_back <- 1
    linputs$no_r_below_start <- FALSE #BEWARE this has been set to TRUE in all earlier runs
    linputs$no_r_below_mut <- FALSE #BEWARE trying out new mutation-selection balance option
    #set plot to F if using grid.arrange as in phaps-rotations-better.Rmd
    linputs$plot <-     FALSE

    ## variable inputs ##
    linputs$cost <-           runif(1, min=0,   max=0.1) 
    #linputs$cost <-           0 
    linputs$expo_hi <-        runif(1, min=0.4, max=0.9)  #min=0.1, max=0.9
    linputs$male_expo_prop <- runif(1, min=0,   max=1)
    linputs$eff <-            runif(1, min=0.5, max=1) # min=0.3, max=1
    linputs$rr  <-            runif(1, min=0.1, max=0.9)
    linputs$dom_sel <-        runif(1, min=0,   max=1)
    linputs$dom_cos <-        runif(1, min=0,   max=1) 
    #TODO make this lognormal to get lower frequencies too
    linputs$start_freqs <-    runif(1, min=0.005,   max=0.1) #note not log yet   
    linputs$n_insecticides <- sample( 2:5, 1 )  #beware set replace=T if more than 1
    # dispersal (initially with cost set to 0)
    #linputs$coverage <-       runif(1, min=0.1, max=0.9)
    #linputs$migration <-      runif(1, min=0.1, max=0.9)


    # create null objects used later for plotting
    ggrot <- ggseq <- NULL

    # to compare sequence to rotation run all the scenarios above for both seq & rot 
    for(rot_or_not in 0:1)
    {  
      #to get rot_int of 0 every other run
      if (rot_or_not == 0) linputs$rot_interval <- 0
      else   #random rot interval integer
          linputs$rot_interval <-   sample( 5:50, 1 )  #beware set replace=T if more than 1

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

      # summarise results per scenario
      # note this is also currently calc at end of run_rot() here is probably better

      res <- gens_under_thresh(dfres)

      # to plot results for checking
      if ( linputs$rot_interval == 0 )
      {
        gens_seq <- as.numeric(res)
      } else
      {
        gens_rot <- as.numeric(res)
      } 


      if (plot_scenarios)
      {
        # make plot
        ggres <- rot_plot_resistance(dfres, plot_refuge=FALSE, logy=FALSE, plot=FALSE)

        # to plot results for checking
        if ( linputs$rot_interval == 0 )
        {
          ggseq <- ggres
          gens_seq <- as.numeric(res)
        } else
        {
          ggrot <- ggres
          gens_rot <- as.numeric(res)
        }
      }

    # TODO I may want to concat all dfres into one object (adda column for scenario num) to save for later analysis
    # e.g. if I wanted to calculate average resistance frequencies over time in addition to gens below threshold

    } #end rot_or_not loop    

    # make an output row for each scenario with inputs, gens_rot, gens_seq
    # need to make sure I retain the rot_interval for rot

    # add scenario id column onto results
    df_in_out <- as.data.frame(c(linputs, gens_rot=gens_rot, gens_seq=gens_seq, id=i))
    df_in_out$rot_minus_seq <- df_in_out$gens_rot - df_in_out$gens_seq

    # previously I used rbind this scenario inputs & results onto previous
    # not a good speed or memory way of doing
    #df_res_all <- rbind(df_res_all, df_in_out)

    # now replaced with this better way
    # set a dataframe for the outputs to the required size in scenario1
    if (i==1)
    {
      df_res_all <- data.frame(matrix(NA, ncol=length(names(df_in_out)), nrow=n_scenarios))
      names(df_res_all) <- names(df_in_out)
    }

    #putting results into one row of overall results dataframe
    df_res_all[i,] <- df_in_out 

    if (plot_scenarios)
    {  
      # plot input vals
      # extract out those between 0 & 1
      dfi <- data.frame(t(as.data.frame(linputs)))
      names(dfi) <- 'value'
      dfi$inputs <- rownames(dfi)

      #take out columns not between 0 & 1
      dfi <- filter(dfi, !inputs %in% c('max_gen','rot_interval','n_insecticides','plot','min_rwr_interval'))
      #also remove coverage & migration both set to 0 for now
      dfi <- filter(dfi, !inputs %in% c('coverage','migration'))

      ggins <- ggplot(dfi, aes(x=inputs,y=value)) + 
                geom_point() +
                theme_minimal() +
                ylim(0,1) +
                theme(axis.text.x = element_text(angle = 20,hjust = 1, vjust = 1))  

      grid.arrange(ggrot, ggseq, ggins, ncol=1)

      cat('scenario',i,'sum generations deployed insecticide under thresh rotation:', gens_rot,' sequence:',gens_seq,"\n\n")

      #can I give the code to rerun this scenario
      cat(paste0("run_rot(n_insecticides=",linputs$n_insecticides,", ",
                  "cost=",signif(linputs$cost,2),", ",
                  "start_freqs=",signif(linputs$start_freqs,2),", ", 
                  "rot_interval=",linputs$rot_interval,", ",
                  "eff=",signif(linputs$eff,2),",\n",
                  "dom_sel=",signif(linputs$dom_sel,2),", ",
                  "dom_cos=",signif(linputs$dom_cos,2),", ",
                  "rr=",signif(linputs$rr,2),", ",
                  "expo_hi=",signif(linputs$expo_hi,2),", ",
                  "coverage=",signif(linputs$coverage,2),",\n",
                  "migration=",signif(linputs$migration,2),", ",
                  "max_gen=",linputs$max_gen,", ", 
                  "min_rwr_interval=",linputs$min_rwr_interval,            
                  ")\n") )

    }

    #I could save partial outputs every 200 scenarios
    #just in case doesn't make to end
    if (i%%1000 == 0)
    {
      #df_res_all$rot_minus_seq <- df_res_all$gens_rot - df_res_all$gens_seq

      # save object containing inputs and results as rda for analysis
      save(df_res_all, file='df_res_all_rot.rda') #paste0(outFolder,'*.rda'))    
    } 

  } #end scenario loop  

}

proportion of scenarios not reaching resistance thresholds

# check proportion of scenarios not reaching thresholds by max_runs for both seq & rot

    nruns_fail <- df_res_all %>% 
                  filter(gens_seq == max_gen & gens_rot == max_gen) %>%
                  summarise(n=n()) %>% 
                  unlist()

    cat('proportion runs not reaching thresholds : ', nruns_fail/nrow(df_res_all) ) 

    # remove runs where thresholds not reached for seq & rot 
    df_res_all <- df_res_all %>% 
                  filter(gens_seq != max_gen | gens_rot != max_gen) 


    nruns_same <- df_res_all %>% 
                  filter(gens_seq == gens_rot) %>%
                  summarise(n=n()) %>% 
                  unlist()

    cat('proportion runs (that did reach threshold) where rotations and sequence the same  : ', nruns_same/nrow(df_res_all) )  


    nruns_rot_20pc_better <- df_res_all %>% 
                  filter(rot_minus_seq/gens_seq >= 0.2 ) %>%
                  summarise(n=n()) %>% 
                  unlist()

    cat('proportion runs (that did reach threshold) where rotations >20 % better than sequence  : ', nruns_rot_20pc_better/nrow(df_res_all) )   

difference between rotations and sequences across all runs

mostly (~67%) no difference but infrequent where rotations better

~ 7% of runs, rotation at least 20% better than sequence

    print( ggplot(df_res_all, aes(x=rot_minus_seq)) + 

             geom_histogram(binwidth=10) +
             xlab('rot_minus_seq : rotation generations minus sequence') +
             ylab('count (sqrt scale to show low values)') +
             scale_y_sqrt(breaks=c(0,10,20,50,200,500,2000,5000)) + #can't use log scale because anchored at 0
             theme_minimal() +
             theme(panel.grid.minor.y=element_blank()))

\pagebreak

Examples of runs where rotations > 20% better than sequences

#eval false because smaller runs don't have these

scenarios_to_plot <- c(624,2167,602,4028,5771)

for(scenario in scenarios_to_plot)
{
  cat(scenario,"\n")

  linputs <- as.list( df_res_all[scenario, 1:15] )

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

  rot_plot_resistance(dfres, plot_refuge=FALSE, logy=FALSE, plot=TRUE)

  #now run as sequence
  linputs$rot_interval <- 0
  dfres <- do.call(run_rot, linputs)

  rot_plot_resistance(dfres, plot_refuge=FALSE, logy=FALSE, plot=TRUE) 


  dfi <- data.frame(t(as.data.frame(linputs)))
  names(dfi) <- 'value'
  dfi$inputs <- rownames(dfi)

  #take out columns not between 0 & 1
  dfi <- filter(dfi, !inputs %in% c('max_gen','rot_interval','n_insecticides','plot','min_rwr_interval'))
  #also remove coverage & migration both set to 0 for now
  dfi <- filter(dfi, !inputs %in% c('coverage','migration'))

  ggins <- ggplot(dfi, aes(x=inputs,y=value)) + 
            geom_point() +
            theme_minimal() +
            ylim(0,1) +
            theme(axis.text.x = element_text(angle = 20,hjust = 1, vjust = 1))
  plot(ggins)

}

How does the difference between rotations and sequences respond to model inputs ?

# go through selected inputs and plot them on x with rot_minus_seq on y 
# initially with 200 reps doesn't show anything obvious, try it with 2000 overnight runs
in_cols <- c('cost','dom_cos','dom_sel','rr','expo_hi','male_expo_prop','eff','start_freqs')

y <- 'rot_minus_seq'  

for(i in in_cols)
{

    #plotnum <- plotnum + 1

    #print( ggplot(ggInsOuts, aes_string(x=i, y=y, colour="strategy")) + 
    #plotlist[[plotnum]] <- 
    print( ggplot(df_res_all, aes_string(x=i, y=y)) + 

             geom_point(shape=1, size=1, alpha=0.5, show.legend=FALSE) +
             theme_minimal() +
             geom_smooth(linetype='dashed',size=0.5) )

}

PRCC

  library(sensitivity)

# pcc(X, y, rank = FALSE, nboot = 0, conf = 0.95)
# 
# Arguments
# X a data frame containing model input variables
# y a vector containing the responses 
# rank logical. If TRUE, the analysis is done on the ranks.
# nboot the number of bootstrap replicates.
# conf the confidence level of the bootstrap confidence intervals.

#same as above
in_cols <- c('cost','dom_cos','dom_sel','rr','expo_hi','male_expo_prop','eff','start_freqs')


    x <- df_res_all[,in_cols]
    y <- df_res_all['rot_minus_seq']

    pcc_res <- pcc(x, y, rank=TRUE)

    #if you add bootstrap, error bars are added to the default plot & extra columns to the PRCC object
    #pcc_res <- pcc(x, y, rank=TRUE, nboot=100)
    #plot(pcc_res)

    #results are here I can probably rbind them together into a df that I can ggplot
    #pcc_res$PRCC

    to_plot <- pcc_res$PRCC
    #rename column 1 from 'original to PRCC
    names(to_plot)[1] <- 'PRCC'
    to_plot$inputs <- rownames(to_plot)  

    axis_text_x <- element_text(angle = 45,hjust = 1, vjust = 1)
    height <- 0.9

    gg <- ggplot( to_plot, aes_string(x='inputs',y='PRCC') ) +     
             geom_point(shape=1, colour='red') +
             theme(axis.text.x = axis_text_x) +
             geom_hline(yintercept = 0, linetype=3) +
             #theme(panel.grid.major.y = element_line(color = "grey", size = 1)) +
             ylim(-1,1) +
             #ggtitle(paste(strat_name[strategy_num])) +
             xlab(NULL) 

    plot(gg)

  #to put back in some ref gridlines that are removed by cowplot
  #theme_set(theme_gray())
  #theme_set(theme_bw())

  #plot_grid( plotlist[[1]],plotlist[[2]],plotlist[[3]],plotlist[[4]],ncol=1, rel_heights=c(1,1,1,1.8))  
  #plot_grid( plotlist[[1]],plotlist[[2]],plotlist[[3]],ncol=1, rel_heights=c(1,1,1.6), labels='AUTO')    

Runs where rotations are >20% better (in red) are associated with higher dominance of cost and lower dominance of selection

but that does not guarantee, must be interaction with other inputs too.

#'cost','dom_cos','dom_sel','rr','expo_hi','male_expo_prop','eff','start_freqs'

#'cost','dom_cos','dom_sel'

y <- 'rot_minus_seq'  

ggplot(df_res_all, aes(x=dom_sel, y=dom_cos)) + 
     geom_point(shape=3, size=1, alpha=0.2, show.legend=FALSE) +
     theme_minimal() +
     geom_point(data=filter(df_res_all, rot_minus_seq/gens_seq >= 0.3 ), shape=1, size=1, col='red', alpha=0.8)


# looking at modelling it
# library(mgcv)
# mod_lm2 <- gam(rot_minus_seq ~ cost + dom_cos + dom_sel, data=df_res_all)
# summary(mod_lm2)  
# 
# mod_gam2 <- gam(rot_minus_seq ~ s(cost) + s(dom_cos) + s(dom_sel) + s(rr) + s(expo_hi) + s(male_expo_prop) + s(eff) + s(start_freqs), data=df_res_all)
# 
# summary(mod_gam2)  
# #Deviance explained = 19.5%
# 
# mod_gam3 <- gam(rot_minus_seq ~ te(cost,dom_cos,dom_sel), data=df_res_all)
# summary(mod_gam3)

#actually decision trees may be useful
#eval=FALSE tmp because changes above will break this
# summary plot of all runs, would show any differences between sequence & rotations as non horizontal lines

# want to plot line between sequence and rotation value
# any deviation from horizontal will be obvious
# could use rotation interval as the x axis, so that longer intervals will be clearer

# maybe not for publication
# but useful for quick view

#it was useful for me but now that I've moved to one row per scenario it won't work !!

    gg <- ggplot( df_res_all, aes_string(x='rot_interval',
                                         y='tot_gens_dep_under50',
                                         group='scenario',
                                         colour='scenario') ) +
          geom_line( alpha=0.5, lwd=1 )

    plot(gg)


# plot gens againts exposure and effectiveness
# as a quick test to see if I can set the starting values of these higher

    gg <- ggplot( df_res_all, aes_string(x='eff',
                                         y='tot_gens_dep_under50') ) +
          geom_point( alpha=0.5 )
    plot(gg)    

    gg <- ggplot( df_res_all, aes_string(x='expo_hi',
                                         y='tot_gens_dep_under50') ) +
          geom_point( alpha=0.5 )
    plot(gg)    

    gg <- ggplot( df_res_all, aes_string(x='expo_hi*eff*male_expo_prop',
                                         y='tot_gens_dep_under50') ) +
          geom_point( alpha=0.5 )
    plot(gg)    

Have a go at running for insecticides with different parameters

starting to work ...

# set random seed so that results can be reproduced
set.seed(2)
max_gen         <- 500
plot_scenarios  <- FALSE
n_scenarios     <- 100 #1000 #000 

  #df_res_all <- data.frame()

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

    # inputs (could also be saved to a dataframe first & extracted one row at a time)
    linputs <- list()

    ## constant inputs ##
    linputs$coverage <- 1
    linputs$migration <- 0
    #linputs$cost <- 0 #variable cost below
    linputs$max_gen <- max_gen

    #beware min_rwr_interval, it's to stop repeated small rotations
    linputs$min_rwr_interval <- 5
    linputs$min_gens_switch_back <- 1
    linputs$no_r_below_start <- FALSE #BEWARE this has been set to TRUE in all earlier runs
    linputs$no_r_below_mut <- TRUE #BEWARE trying out new mutation-selection balance option
    #set plot to F if using grid.arrange as in phaps-rotations-better.Rmd
    linputs$plot <-     FALSE

    ## variable inputs ##
    linputs$n_insecticides <- sample( 2:5, 1 )  #beware set replace=T if more than 1
    n_ins <- linputs$n_insecticides

    ## new those below can be different for each insecticide
    linputs$cost <-           runif(n_ins, min=0,   max=0.1) 
    #linputs$cost <-           0 
    linputs$expo_hi <-        runif(n_ins, min=0.4, max=0.9)  #min=0.1, max=0.9
    linputs$male_expo_prop <- runif(1, min=0,   max=1)
    linputs$eff <-            runif(n_ins, min=0.5, max=1) # min=0.3, max=1
    linputs$rr  <-            runif(n_ins, min=0.1, max=0.9)
    linputs$dom_sel <-        runif(n_ins, min=0,   max=1)
    linputs$dom_cos <-        runif(n_ins, min=0,   max=1) 
    #DIFFERENT
    linputs$start_freqs <-    runif(n_ins, min=0.0001,   max=0.01) #note not log yet   

    # dispersal (initially with cost set to 0)
    #linputs$coverage <-       runif(1, min=0.1, max=0.9)
    #linputs$migration <-      runif(1, min=0.1, max=0.9)


    # create null objects used later for plotting
    ggrot <- ggseq <- NULL

    # to compare sequence to rotation run all the scenarios above for both seq & rot 
    for(rot_or_not in 0:1)
    {  
      #to get rot_int of 0 every other run
      if (rot_or_not == 0) linputs$rot_interval <- 0
      else   #random rot interval integer
          linputs$rot_interval <-   sample( 5:50, 1 )  #beware set replace=T if more than 1

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

      # summarise results per scenario
      # note this is also currently calc at end of run_rot() here is probably better

      res <- gens_under_thresh(dfres)

      # to plot results for checking
      if ( linputs$rot_interval == 0 )
      {
        gens_seq <- as.numeric(res)
      } else
      {
        gens_rot <- as.numeric(res)
      } 


      if (plot_scenarios)
      {
        # make plot
        ggres <- rot_plot_resistance(dfres, plot_refuge=FALSE, logy=FALSE, plot=FALSE)

        # to plot results for checking
        if ( linputs$rot_interval == 0 )
        {
          ggseq <- ggres
          gens_seq <- as.numeric(res)
        } else
        {
          ggrot <- ggres
          gens_rot <- as.numeric(res)
        }
      }
    } #end rot_or_not loop    

    # make an output row for each scenario with inputs, gens_rot, gens_seq
    # need to make sure I retain the rot_interval for rot

    # to cope with multiple values per scenario convert them to strings
    linputs2 <- lapply(linputs, function(x) ifelse(length(x)>1,toString(x),x))

    #added scenario id column
    df_in_out <- as.data.frame(c(linputs2, gens_rot=gens_rot, gens_seq=gens_seq, id=i), stringsAsFactors=FALSE)

    df_in_out$rot_minus_seq <- df_in_out$gens_rot - df_in_out$gens_seq


    #instead set up the dataframe to the required size in scenario1
    if (i==1)
    {
      df_res_all <- data.frame(matrix(NA, ncol=length(names(df_in_out)), nrow=n_scenarios))
      names(df_res_all) <- names(df_in_out)
    }

    #putting results into one row of overall results dataframe
    df_res_all[i,] <- df_in_out   

    if (plot_scenarios)
    {  
      # plot input vals
      # extract out those between 0 & 1
      dfi <- data.frame(t(as.data.frame(linputs)))
      names(dfi) <- 'value'
      dfi$inputs <- rownames(dfi)

      #take out columns not between 0 & 1
      dfi <- filter(dfi, !inputs %in% c('max_gen','rot_interval','n_insecticides','plot','min_rwr_interval'))
      #also remove coverage & migration both set to 0 for now
      dfi <- filter(dfi, !inputs %in% c('coverage','migration'))

      ggins <- ggplot(dfi, aes(x=inputs,y=value)) + 
                geom_point() +
                theme_minimal() +
                ylim(0,1) +
                theme(axis.text.x = element_text(angle = 20,hjust = 1, vjust = 1))  

      grid.arrange(ggrot, ggseq, ggins, ncol=1)

      cat('scenario',i,'sum generations deployed insecticide under thresh rotation:', gens_rot,' sequence:',gens_seq,"\n\n")

      #can I give the code to rerun this scenario
      cat(paste0("run_rot(n_insecticides=",linputs$n_insecticides,", ",
                  "cost=",signif(linputs$cost,2),", ",
                  "start_freqs=",signif(linputs$start_freqs,2),", ", 
                  "rot_interval=",linputs$rot_interval,", ",
                  "eff=",signif(linputs$eff,2),",\n",
                  "dom_sel=",signif(linputs$dom_sel,2),", ",
                  "dom_cos=",signif(linputs$dom_cos,2),", ",
                  "rr=",signif(linputs$rr,2),", ",
                  "expo_hi=",signif(linputs$expo_hi,2),", ",
                  "coverage=",signif(linputs$coverage,2),",\n",
                  "migration=",signif(linputs$migration,2),", ",
                  "max_gen=",linputs$max_gen,", ", 
                  "min_rwr_interval=",linputs$min_rwr_interval,            
                  ")\n") )

    }

    #save partial outputs just in case doesn't make to end
    if (i%%1000 == 0)
    {
      #df_res_all$rot_minus_seq <- df_res_all$gens_rot - df_res_all$gens_seq

      # save object containing inputs and results as rda for analysis
      save(df_res_all, file='df_res_all_rot.rda') #paste0(outFolder,'*.rda'))    
    } 
  } #end scenario loop  


#to allow plotting of chosen scenarios by re-running them 
#initially didn't work for diff insecticides because the inputs are saved as strings 
#now the strings get converted back to vectors below

scenarios_to_plot <- c(21,23,32) #big difference ones all have insecticide1 as the good one
scenarios_to_plot <- c(29,38,46)
scenarios_to_plot <- c(13)

for(scenario in scenarios_to_plot)
{
  #cat(scenario,"\n")
  linputs <- as.list( df_res_all[scenario, 1:15] )

  #converting multiple insecticide inputs back from strings
  #have to ignore the ones that aren't strings

  #note the space after the comma in split=", "
  linputs <- lapply(linputs, function(x) if (!grepl(",",x)) x else
                                          as.numeric(unlist(strsplit(as.character(x),split=", "))))

  #TESTING changing min_rwr_interval from 5
  #if min_rwr_interval is set to 0 you get frequent switch backs to other insecticides for 1 generation
  #but when it's set to 5 it pushes frequences above the resistance freq threshold
  linputs$min_rwr_interval <- 5
  #min_gens_switch_back has to be 0 to allow one insecticide to continue when all others dead
  linputs$min_gens_switch_back = 0

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

  ggrot <- rot_plot_resistance(dfres, plot_refuge=FALSE, logy=TRUE, plot=FALSE)

  # now run as sequence
  linputs$rot_interval <- 0
  dfres <- do.call(run_rot, linputs)

  ggseq <- rot_plot_resistance(dfres, plot_refuge=FALSE, logy=TRUE, plot=FALSE) 

  grid.arrange(ggrot, ggseq, ncol=1)

}  


ian-hastings/rotations documentation built on Dec. 5, 2019, 3:14 a.m.