R/testHostSwitch.R

In HostSwitch: Simulate the Extent of Host Switching by Consumers

#' Test for the significance of the difference between two difference of two
#' HostSwitch objects
#'
#' @param simulated_quantities1 An object created by [simHostSwitch()]
#' @param simulated_quantities2 An object created by [simHostSwitch()]
#' @param parameter Quantity of interest, possible values are:
#' \itemize{
#'   \item "j" for total number of dispersing events (or **j**umps)
#'   \item "s" for total number of **s**uccessful host switch events
#'   \item "d" for **d**istance between the pRes_sim and pRes_new_sim in
#'         case of host switch
#' }
#' @param test Statistical test, available tests are:
#' \itemize{
#'   \item ''t'' for **t**-test (parametric)
#'   \item ''w'' for **W**ilcoxon-test (non-parametric)
#' }
#' @param warmup Number of warmup steps to be excluded when comparing models,
#' see details. Possible value are NULL (default) or
#' positive integer (min=1,max=50).
#' @param plot If *TRUE*, a boxplot is drawn. Default = FALSE.
#' @details This function tests the significance of the difference between
#' two objects generated by simHostSwitch function.
#' Warmup represents the initial condition that we want to exclude from the
#' test. The initial condition corresponds to the number of
#' generations (n_generations): warmup = 1 means that the generation at time 0
#' is excluded from comparison; warmup = 2 means generations at times 0 and 1
#' are excluded and so on.
#' If warmup = NULL all generations are considered for comparison, i.e. initial
#' condition is not considered.
#' @return An object of class testHostSwitch
#' @examples
#' m1 = simHostSwitch(n_generations=100,n_sim=100)
#' m2 = simHostSwitch(n_generations=50,n_sim=50)
#' testHostSwitch(simulated_quantities1=m1,simulated_quantities2=m2,
#' parameter="j",test="t",plot=TRUE)
#' @export

testHostSwitch = function(simulated_quantities1,simulated_quantities2,
                          parameter,test,warmup=NULL,plot=FALSE){

  # input checks
  checkmate::assert_class(simulated_quantities1,"HostSwitch") #class HostSwitch
  checkmate::assert_class(simulated_quantities2,"HostSwitch") #class HostSwitch
  checkmate::assertChoice(parameter, c("d", "j","s")) # parameter
  checkmate::assertChoice(test, c("w","t")) # test
  checkmate::assertChoice(plot, c(TRUE,FALSE)) # plot
  checkmate::assertCount(warmup,positive=TRUE,null.ok = TRUE)
  checkmate::assertNumeric(warmup,upper=50,null.ok = TRUE) # warmup


  values <- NULL # global variables

  out=list()
  n_sim1 = simulated_quantities1$n_sim; n_sim2 = simulated_quantities2$n_sim

  if (n_sim1 <2){
    warning("The number of simulations for 'simulated_quantities1' is < 2. For comparisons, we recommend at least n_sim = 2", call. = F)
  }

  if (n_sim2 <2){
    warning("The number of simulations for 'simulated_quantities2' is < 2. For comparisons, we recommend at least n_sim = 2", call. = F)
  }

  n_generations1 = simulated_quantities1$n_generations
  n_generations2 = simulated_quantities2$n_generations

  if (length(warmup)>0){
    # model 1
    simulated_quantities1[["pRes_sim"]]      = lapply(
      simulated_quantities1[["pRes_sim"]], function(x) x[-c(1:warmup)])
    simulated_quantities1[["pRes_new_sim"]]  = lapply(
      simulated_quantities1[["pRes_new_sim"]], function(x) x[-c(1:warmup)])
    simulated_quantities1[["pInd_jump_sim"]] = lapply(
      simulated_quantities1[["pInd_jump_sim"]], function(x) x[-c(1:warmup)])
    simulated_quantities1[["pInd_whichsurv_sim"]] = lapply(
      simulated_quantities1[["pInd_whichsurv_sim"]],
      function(x) x[-c(1:warmup)])
    for (i in 1:n_sim1){
      simulated_quantities1[["pInd_sim"]][[i]] =
        simulated_quantities1[["pInd_sim"]][[i]][-c(1:warmup)]
    }
    n_generations1 = n_generations1 - warmup

    # model 2
    simulated_quantities2[["pRes_sim"]]      = lapply(
      simulated_quantities2[["pRes_sim"]], function(x) x[-c(1:warmup)])
    simulated_quantities2[["pRes_new_sim"]]  = lapply(
      simulated_quantities2[["pRes_new_sim"]], function(x) x[-c(1:warmup)])
    simulated_quantities2[["pInd_jump_sim"]] = lapply(
      simulated_quantities2[["pInd_jump_sim"]], function(x) x[-c(1:warmup)])
    simulated_quantities2[["pInd_whichsurv_sim"]] = lapply(
      simulated_quantities2[["pInd_whichsurv_sim"]],
      function(x) x[-c(1:warmup)])
    for (i in 1:n_sim2){
      simulated_quantities2[["pInd_sim"]][[i]] =
        simulated_quantities2[["pInd_sim"]][[i]][-c(1:warmup)]
    }
    n_generations2 = n_generations2 - warmup
  }





# choice of parameter
## JUMPS
  if(parameter == "j"){
    title = "Comparison of averages of dispersion events by the Consumer"
   x = unlist(lapply(simulated_quantities1$pInd_jump_sim,
                   function(x) length(which(x>0))))
   y = unlist(lapply(simulated_quantities2$pInd_jump_sim,
                   function(x) length(which(x>0))))
  }
## SUCCESSFUL HOST SWITCHES
  if(parameter == "s"){
    title = "Comparison of averages of host switch by the Consumer"
    x = rep(0,n_sim1)
    y = rep(0,n_sim2)

    for (i in 1:n_sim1){
      dat = lapply(simulated_quantities1[c(1,2)], `[[`, i)
      x[i] = length(which(dat$pRes_sim[-1]==dat$pRes_new_sim))
    }

    for (i in 1:n_sim2){
      dat = lapply(simulated_quantities2[c(1,2)], `[[`, i)
      y[i] = length(which(dat$pRes_sim[-1]==dat$pRes_new_sim))
    }
  }



## Distance between mean parasite and new host phenotype
  if(parameter == "d"){
    title = "Comparison of averages of phenotypic distance for consumers
             in current and novel Resource"

    # compare pRes and pRes_new in case of successful jump

    ## sim1
    survPosition1 = list()
    for (i in 1:n_sim1){
      dat = simulated_quantities1$pInd_whichsurv_sim[i]
      dat = flatten2(dat)
      survPosition1[[i]] = (which(
        unlist(lapply(dat,function(x) length(which(x>0))))>0))-1
    }

    ## pRes
    if(n_sim1>1){
    pRes_when_Survived =mapply(FUN = function(x,y) {d <- x[y]},
                               x = simulated_quantities1$pRes_sim,
                               # get only pRes when jump occurred
                               y = survPosition1,SIMPLIFY = FALSE)
    pRes=unlist(flatten2(pRes_when_Survived))

    ## pRes_new
    pRes_new_when_Survived =mapply(FUN = function(x,y) {d <- x[y]},
                                   x = simulated_quantities1$pRes_new_sim,
                                   # get only pRes when jump occurred
                                   y = survPosition1,SIMPLIFY = FALSE)
    pRes_new=unlist(flatten2(pRes_new_when_Survived))
    #if(is.matrix(pRes_new_when_Survived)){pRes_new=pRes_new_when_Survived}
    }
    # 7, 228, [249,]
    if(n_sim1==1){
      pRes =  unlist(simulated_quantities1$pRes_sim)[unlist(survPosition1)]
      pRes_new =  unlist(simulated_quantities1$pRes_new_sim)[unlist(survPosition1)]
    }
    x = abs(pRes-pRes_new)




    # sim2
    survPosition2= list()

    for (i in 1:n_sim2){
      # for each simulation, records position (=generation) of survived individuals
      dat = simulated_quantities2$pInd_whichsurv_sim[i]
      dat = flatten2(dat)
      survPosition2[[i]] = (which(unlist(lapply(dat,function(x) length(which(x>0))))>0))-1
    }

    ## pRes
    if(n_sim2>1){
    pRes_when_Survived =mapply(FUN = function(x,y) {d <- x[y]},
                               x = simulated_quantities2$pRes_sim,
                               # get only pRes when jump occurred
                               y = survPosition2,SIMPLIFY = FALSE)
    pRes=unlist(flatten2(pRes_when_Survived))
    #if(is.matrix(pRes_when_Survived)){pRes=pRes_when_Survived}

    ## pRes_new
    pRes_new_when_Survived =mapply(FUN = function(x,y) {d <- x[y]},
                                   x = simulated_quantities2$pRes_new_sim,
                                   # get only pRes when jump occurred
                                   y = survPosition2,SIMPLIFY = FALSE)
    pRes_new=unlist(flatten2(pRes_new_when_Survived))
}
    if(n_sim2==1){
      pRes =  unlist(simulated_quantities2$pRes_sim)[unlist(survPosition2)]
      pRes_new =  unlist(simulated_quantities2$pRes_new_sim)[unlist(survPosition2)]
    }

    y = abs(pRes-pRes_new)
  }

# statistical test
if(test == "t"){
  out$result = stats::t.test(x,y)
}

if(test == "w"){
    out$result = stats::wilcox.test(x,y)
}

df1 = data.frame(x=rep("sim1",length(x)),values=x) # simulated_quantities1
df2 = data.frame(x=rep("sim2",length(y)),values=y) # simulated_quantities2
plotInput = data.frame(rbind(df1,df2))
out$result=append(out$result,plotInput)



# plot
if(plot == TRUE){

  g=ggplot2::ggplot(data=plotInput,aes(x=x,y=values,group=x)) +
    geom_boxplot() +
  labs(title = title) +
  ggplot2::theme_bw()
 print(g)
}

methods::new("testHostSwitch", out$result)
}

#' Set class and method
#'
#' This is a build-time dependency on methods, as opposed to a run-time
#' dependency, thus requiring the importFrom tag to avoid a NOTE when checking
#' the package on CRAN.
#'
#' @keywords internal
#' @export
#'
methods::setClass("testHostSwitch", representation=representation("list"))
methods::setMethod("show",signature = "testHostSwitch",
                   definition = function(object) {
  cat("An object of class ", class(object), "\n", sep = "")
  cat("Test result comparing 2 HostSwitch simulations using",
      object$method,"\n\n")
  cat(names(object$statistic),": ",object$statistic,", df:",object$parameter,
      ", p.value:",object$p.value,"\n\n",sep="")
  invisible(NULL)
})