R/multisim.R

In rucky151/shapp1: Scenario Analysis for Seed Health

#' Create a matrix with output from multiple simulations for a single parameter combination
#'
#' This function creates a matrix for output from multiple simulations for a single parameter combination across nseasons
#'
#' The truncated random normal variables were generated by \code{\link{altruncnorm}}. Refer \code{\link{onesim}} and \code{\link{multipar}}for relevant information.
#'
#' @inheritParams onesim
#' @param nsim the number of simulations, numeric or numeric vector

#' @keywords seed health
#' @importFrom stats median quantile rnorm var
#' @export
#' @examples
#' multisim() # to be added


# to do - GENERAL TESTING
# to do - check whether parameter list is correct

# Columns of output matrix
# col 1 - fHP healthy plant number
# col 2 - fDP diseased plant number (after roguing)
# col 3 - fHS healthy seed number
# col 4 - fDS diseased seed number
# col 5 - fpHS proportion healthy seed
# col 6 - fpDS proportion diseased seed
# col 7 - HPtrans Season in which HP first transitions below HPcut*Kx, if HPtrans is NA i.e., HP never less than HPcut, set to max seasons tested
# col 8 - pHStrans Season in which pHS first transitions below pHScut, if pHStrans is NA i.e., pHS never less than pHScut, set to max seasons tested
# col 9 - HPpseas Healthy plants are calculated from season 1 onwards
# col 10 - pHSpseas Proportion seasons with pHS below pHScut
# col 11 - fYld end of season yield
# col 12 - fYL end of season yield loss
# col 13 - DPbr diseased plants before roguing

# Step 1B. Create matrix for output from one simulation (stochastic model)



# set.seed(1234)


multisim <- function(pHSinit=0.8, Kx = 100, betax=0.02, wxtnormm=0.8, wxtnormsd=0.3, hx=1, mxtnormm=1, 
                     mxtnormsd=0.1, axtnormm=1, axtnormsd=0.1, rx=0.1, zxtnormm=1, zxtnormsd= 0.1, gx=4, 
                     cx=0.9, phix=0, nseasons=10, nsim=10, HPcut=0.5,pHScut=0.5,maY=100, miY=0, thetax=0.2, Ex=0) {

  # nsim - number of simulations

  outmf <- as.data.frame(matrix(data=-999, nrow=nsim, ncol=12, 
                                dimnames = list(1:nsim,c('fHP', 'fDP', 'fHS', 'fDS', 'fpHS', 'fpDS', 'HPtrans', 'pHStrans', 'HPpseas', 'pHSpseas', 'fYld', 'fYL'))))

  for(si in 1:nsim) {
    temp <- onesim(pHSinit=pHSinit, Kx = Kx, betax=betax, wxtnormm=wxtnormm, wxtnormsd=wxtnormsd, 
                   hx=hx, mxtnormm=mxtnormm, mxtnormsd=mxtnormsd, axtnormm=axtnormm, axtnormsd=axtnormsd, 
                   rx=rx,zxtnormm=zxtnormm, zxtnormsd= zxtnormsd, gx=gx, cx=cx, phix=phix, nseasons=nseasons, 
                   HPcut=HPcut, pHScut=pHScut, maY=maY, miY=miY, thetax=thetax, Ex=Ex)$outfin
    
    outmf$fHP[si] <- temp$HP
    outmf$fDP[si] <- temp$DP
    outmf$fHS[si] <- temp$HS
    outmf$fDS[si] <- temp$DS
    outmf$fpHS[si] <- temp$pHS
    outmf$fpDS[si] <- temp$pDS
    outmf$HPtrans[si] <- temp$HPtrans
    outmf$pHStrans[si] <- temp$pHStrans
    outmf$HPpseas[si] <- temp$HPpseas
    outmf$pHSpseas[si] <- temp$pHSpseas
    outmf$fYld[si]<-temp$Yld
    outmf$fYL[si]<-temp$YL
  }

  quantile0.05 <- function(x){quantile(x,probs=0.05,na.rm=T)}
  quantile0.95 <- function(x){quantile(x,probs=0.95,na.rm=T)}

  outfsum <- as.data.frame(matrix(data=-999, nrow=5, ncol=12, 
                                  dimnames = list(c('mean','median','var','q0.05','q0.95'), c('fHP', 'fDP', 'fHS', 'fDS', 'fpHS', 'fpDS',  'HPtrans', 'pHStrans', 'HPpseas', 'pHSpseas', 'fYld', 'fYL'))))

  # the first row gives the mean for each of the responses, the second the median, etc.
  outfsum[1,] <- apply(outmf,MARGIN=2,FUN=mean)
  outfsum[2,] <- apply(outmf,MARGIN=2,FUN=median)
  outfsum[3,] <- apply(outmf,MARGIN=2,FUN=var)
  outfsum[4,] <- apply(outmf,MARGIN=2,FUN=quantile0.05)
  outfsum[5,] <- apply(outmf,MARGIN=2,FUN=quantile0.95)

  #-----------------------------------------------
  # warning message
  if ( pHSinit < 0 | pHSinit > 1){
    warning(paste('pHSinit: your input value is', pHSinit,', it must be between 0 and 1'))
  } else if (betax < 0.001 ) {
    warning(paste('betax: your input value is', betax,', it must be greater than or equal to 0.001'))
  } else if (wxtnormm < 0 | wxtnormm > 1) {
    warning(paste('wxtnormm: your input value is', wxtnormm,', it must be between 0 and 1'))
  } else if (hx < 0 | hx > 1) {
    warning(paste('hx: your input value is', hx,', it must be between 0 and 1'))
  } else if (mxtnormm < 0 | mxtnormm > 1) {
    warning(paste('mxtnormm: your input value is', mxtnormm,', it must be between 0 and 1'))
  } else if (axtnormm < 0 | axtnormm > 1) {
    warning(paste('axtnormm: your input value is', axtnormm,', it must be between 0 and 1'))
  } else if (rx < 0 | rx > 1) {
    warning(paste('rx: your input value is', rx,', it must be between 0 and 1'))
  } else if (zxtnormm < 0 | zxtnormm > 1) {
    warning(paste('zxtnormm: your input value is', zxtnormm,', it must be between 0 and 1'))
  } else if (gx < 0 ) {
    warning(paste('gx: your input value is', gx,', it must be 0 or a positive integer'))
  } else if (cx < 0 | cx > 1) {
    warning(paste('cx: your input value is', cx,', it must be between 0 and 1'))
  } else if (phix < 0 | phix > 1) {
    warning(paste('phix: your input value is', phix,', it must be between 0 and 1'))
  } else if (thetax < -1 | thetax > 1) {
    warning(paste('thetax: your input value is', thetax,', it must be between -1 and 1'))
  } else if (Ex < 0 ) {
    warning(paste('Ex: your input value is', Ex,', it must be greater than or equal to 0'))
  }else {

    list(outmf=outmf,outfsum=outfsum)
  }
  #-----------------------------------------------

}