robustSPM {MQMF}

robustSPM does a robustness test on quality of fit of an SPM

Description

robustSPM conducts a robustness test on the quality of fit of an SPM. This is done by using the original optimal model parameters or the original guessed parameter values, add random variation to each of them, and re-fit the model. This process needs to be repeated multiple times. This should enable an analysis of the stability of the modelling outcomes. If the optimum parameters are used then add more variation, if initial guesses are used you may need to select different starting points so that the random variation covers the parameter space reasonably well.

Usage

robustSPM(
  inpar,
  fish,
  N = 10,
  scaler = 40,
  verbose = FALSE,
  schaefer = TRUE,
  funk = simpspm,
  funkone = FALSE,
  steptol = 1e-06
)

Arguments

inpar

the parameter set with which to begin the trials

fish

the fisheries data: at least year, catch, and cpue

N

number of random trials to run; default = 10 = not enough

scaler

the divisor that sets the degree of normal random variation to add to the parameter values; default = 15 the smaller the value the more variable the outcome

verbose

progress and summary statistics to the screen? default = FALSE

schaefer

default = TRUE, which sets the analysis to the Schaefer model. setting it to FALSE applies the approximate Fox model

funk

the function used to generate the predicted cpue

funkone

defaults=FALSE; use negLL or negLL1, with FALSE robustSPM will use negLL, with TRUE it will use negLL1 which has a constraint on the first parameter to keep it > 0

steptol

is the steptol from nlm as used in fitSPM, the default value is 1e-06, as usual.

Value

a list containing the results from each run, the range of values across runs, and the median values.

Examples

  data(dataspm)
  param <- log(c(r=0.24,K=5174,Binit=2846,sigma=0.164))
  ans <- fitSPM(pars=param,fish=dataspm,schaefer=TRUE,maxiter=1000)
  out <- robustSPM(ans$estimate,dataspm,N=5,scaler=40,verbose=TRUE,
                   schaefer=TRUE) # N should be 50, 100, or more
  str(out)
  print(out$results)
  pairs(out$results[,c(6:8,11)]) # need a larger N!