R/testHostSwitch.R

Defines functions testHostSwitch

Documented in testHostSwitch

#' Test for the significance of the difference between two difference of two HostSwitch objects
#'
#' @param simulated_quantities1 An object created by \code{\link{simHostSwitch}}
#' @param simulated_quantities2 An object created by \code{\link{simHostSwitch}}
#' @param parameter Quantity of interest, possible values are:
#' \itemize{
#'   \item "j" for total number of dispersing events (or \strong{j}umps)
#'   \item "s" for total number of \strong{s}uccessful host switch events
#'   \item "d" for \strong{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 \strong{t}-test (parametric)
#'   \item ''w'' for \strong{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 \emph{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)
#' @import checkmate
#' @import stats
#' @import ggplot2
#' @importFrom purrr map
#' @importFrom plyr laply
#' @importFrom utils head
#' @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("t","w")) # 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
  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 = plyr::laply(simulated_quantities1$pInd_jump_sim,function(x) length(which(x>0)))
   y = plyr::laply(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(plyr::laply(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, y = survPosition1,SIMPLIFY = FALSE) # get only pRes when jump occurred
    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_quantities1$pRes_new_sim, y = survPosition1,SIMPLIFY = FALSE) # get only pRes when jump occurred
    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(plyr::laply(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, y = survPosition2,SIMPLIFY = FALSE) # get only pRes when jump occurred
    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, y = survPosition2,SIMPLIFY = FALSE) # get only pRes when jump occurred
    pRes_new=unlist(flatten2(pRes_new_when_Survived))
    #if(is.matrix(pRes_new_when_Survived)){pRes_new=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() +  ##replaced data=out$plotInput with data=plotInput
  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
#' @importFrom methods setClass setMethod
#' @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)
})

Try the HostSwitch package in your browser

Any scripts or data that you put into this service are public.

HostSwitch documentation built on Sept. 27, 2021, 9:06 a.m.