shapp1: Scenario Analysis for Seed Health

Documented in multipar

#' Create a matrix with output for multiple parameter combinations
#'
#' This function creates a matrix with output for multiple parameter combinations across nseasons
#'
#' The truncated random normal variables were generated by \code{\link{altruncnorm}}. Refer \code{\link{onesim}}, and \code{\link{multisim}} for relevant information.
#'dde
#' @inheritParams onesim
#' @param nsim vector of number of simulations
#' @import tibble
#' @importFrom stats median quantile rnorm var
#' @keywords seed health
#' @examples
#' multipar() # to be added
#' @export
#'





# 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

# mean (mean), median (median),  variance (var), quantile 0.05 (0.05), quantile 0.95 (0.95)

# set.seed(1234)

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

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


ncomb <- length(pHSinit) * length(Kx) * length(betax) * length(wxtnormm) * length(wxtnormsd)*length(hx) * length(mxtnormm)*length(mxtnormsd) * length(axtnormm) * length(axtnormsd) * length(rx) * length(zxtnormm) * length(zxtnormsd) * length(gx) * length(cx) * length(phix) * length(nseasons) * length(nsim) * length(HPcut) * length(pHScut) * length(maY)*length(miY)*length(thetax)*length(Ex)

basn <- c('fHP', 'fDP', 'fHS', 'fDS', 'fpHS', 'fpDS', 'HPtrans', 'pHStrans', 'HPpseas', 'pHSpseas', 'fYld', 'fYL')

outmpc <- as.data.frame(matrix(data=-999, nrow=ncomb, ncol=84, dimnames = list(1:ncomb,c('pHSinit','Kx','betax','wxtnormm','wxtnormsd','hx','mxtnormm','mxtnormsd','axtnormm','axtnormsd','rx','zxtnormm','zxtnormsd','gx','cx','phix','nseasons','nsim', 'HPcut', 'pHScut', 'maY','miY','thetax','Ex', paste(basn,'mean',sep=''),paste(basn,'median',sep=''),paste(basn,'var',sep=''),paste(basn,'q0.05',sep=''),paste(basn,'q0.95',sep='')))))

icomb <- 1 # indicates number of current parameter combination, and corresponding row of outmpc

for(i1 in 1:length(pHSinit)){ tpHSinit <- pHSinit[i1]
for(i2 in 1:length(Kx)){tKx <- Kx[i2]
for(i3 in 1:length(betax)){tbetax <- betax[i3]
for(i4 in 1:length(wxtnormm)){twxtnormm <- wxtnormm[i4]
for(i5 in 1:length(wxtnormsd)){twxtnormsd<- wxtnormsd[i5]
for(i6 in 1:length(hx)){thx <- hx[i6]
for(i7 in 1:length(mxtnormm)){tmxtnormm <- mxtnormm[i7]
for(i8 in 1:length(mxtnormsd)){tmxtnormsd <- mxtnormsd[i8]
for(i9 in 1:length(axtnormm)){taxtnormm <-axtnormm[i9]
for(i10 in 1:length(axtnormsd)){taxtnormsd <- axtnormsd[i10]
for(i11 in 1:length(rx)){trx <- rx[i11]
for(i12 in 1:length(zxtnormm)){tzxtnormm <-zxtnormm[i12]
for(i13 in 1:length(zxtnormsd)){tzxtnormsd <- zxtnormsd[i13]
for(i14 in 1:length(gx)){tgx <- gx[i14]
for(i15 in 1:length(cx)){tcx <- cx[i15]
for(i16 in 1:length(phix)){tphix <- phix[i16]
for(i17 in 1:length(nseasons)){tnseasons <- nseasons[i17]
for(i18 in 1:length(nsim)){tnsim <- nsim[i18]
for(i19 in 1:length(HPcut)){tHPcut<- HPcut[i19]
for(i20 in 1:length(pHScut)){tpHScut<- pHScut[i20]
for(i21 in 1:length(maY)){tmaY<- maY[i21]
for(i22 in 1:length(miY)){tmiY<- miY[i22]
for(i23 in 1:length(thetax)){tthetax<- thetax[i23]
for(i24 in 1:length(Ex)){tEx<-Ex[i24]


temp <- multisim(tpHSinit, tKx, tbetax, twxtnormm, twxtnormsd, thx, tmxtnormm, tmxtnormsd, taxtnormm, taxtnormsd,
                 trx, tzxtnormm, tzxtnormsd, tgx, tcx, tphix, tnseasons, tnsim, tHPcut, tpHScut, tmaY, tmiY,  tthetax,
                 tEx)$outfsum

outmpc[icomb,1:24] <- c(tpHSinit,tKx,tbetax,twxtnormm,twxtnormsd,thx,tmxtnormm, tmxtnormsd,taxtnormm,taxtnormsd,trx,tzxtnormm,tzxtnormsd,tgx,tcx, tphix,tnseasons,tnsim,tHPcut,tpHScut,tmaY,tmiY,tthetax, tEx)
outmpc[icomb,25:36] <- temp[1,]
outmpc[icomb,37:48] <- temp[2,]
outmpc[icomb,49:60] <- temp[3,]
outmpc[icomb,61:72] <- temp[4,]
outmpc[icomb,73:84] <- temp[5,]

icomb <- icomb + 1
}}}}}}}}}}}}}}}}}}}}}}}}

#---------------------------------------
# warning message

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

#---------------------------------------

}