Function predRMSE

Let’s start with the function for the RMSE, predRMSE. We can define it as a simple function by using the ellipsis ... as the (only) input argument, and inside this function we can pass it on to the function redictionError. Via the ellipsis ... we can pass input arguments to the predictionError function without having to explicitly specify them in our definition of predRMSE. Obviously, we could explicitly do so, but in this case this is a short and convenient way to define our function. Thus: the predRMSE function takes input via the ..., which it passes on to predictionError, which returns the prediction errors. Then, we compute the RMSE from these errors and return it:

# Define a function that computes the Root Mean Squared Error, RMSE
# NOTE: nothing needs to be filled in here, so this function can be used as-is
predRMSE <- function(...) { # NOTE: these ... are ellipsis, so DO NOT fill in something here!
  # Get prediction errors
  err <- predictionError(...) # NOTE: these ... are ellipsis, so DO NOT fill in something here!
  
  # Compute the measure of model fit (here the Root Mean Squared Error)
  RMSE <- sqrt(mean(err^2))
  
  # Return the measure of model fit
  return(RMSE)
}

Let’s test whether this function works:

predRMSE(p = ParmsLogistic(),
         y0 = InitLogistic(),
         fn = RatesLogistic,
         name = "N",
         obs_value = surveyData$N,
         obs_times = surveyData$t) # optionally add values for out_time

## [1] 40.91814

Indeed, the value is similar to what we computed above (as it should be)!

Transformations

For the transformations and back-transformations, we have coded a simple function for you to facilitate these steps. This function allows for 3 types of transformations, which are among the most common types of transformations used (but there are more):

• “none”: no transformation, we simply leave the parameter value as-is. This should be the case for model parameters that can have any value $> -\infty$ and $< \infty$
• “log”: logarithmic transformation: this is a good transformation for parameters whose support is $>= 0$ (e.g. the $r$ and $K$ parameters in our logistic growth model);
• “logit”: logit transformation: this is a convenient transformation for parameters whose support is $>= 0$ and $<= 1$ (e.g. fractions or probabilities).

The function is called transformParameters and takes as input:

• parm: the vector with parameters, e.g., a vector as returned by the ParmsLogistic function;
• trans: a character vector with transformations: it should be of the same length as the parameter vector, and possible values for each parameter are “none”, “log” and “logit”.
• rev: a logical value (i.e., TRUE or FALSE), specifying whether or not we and the back-transformation. If FALSE, we transform the parameter from its own scope to unconstrained scope; if TRUE we back-transform from the unconstrained scale to the parameter scale:

# Define a function that allows for transformation and back-transformation
# NOTE: nothing needs to be filled in here, so this function can be used as-is
transformParameters <- function(parm, trans, rev = FALSE) {
  # Check transformation vector
  trans <- match.arg(trans, choices = c("none","logit","log"), several.ok = TRUE)
  
  # Get transformation function per parameter
  if(rev == TRUE) {
    # Back-transform from real scope to parameter scope
    transfun <- sapply(trans, switch,
                       "none" = mean,
                       "logit" = plogis,
                       "log" = exp)
  }else {
    # Transform from parameter scope to real scope
    transfun <- sapply(trans, switch,
                       "none" = mean,
                       "logit" = qlogis,
                       "log" = log)
  }
  
  # Apply transformation function
  y <- mapply(function(x, f){f(x)}, parm, transfun)
  
  # Return
  return(y)
}

Let’s test this function for our logistic growth model, for parameters $r$ and $K$, both via a “log” transform, and with values as defined by function ParmsLogistic:

prm_t <- transformParameters(parm = ParmsLogistic(), trans = c("log","log"), rev = FALSE)
prm_t

##         r         K 
## -1.203973  4.605170

We can now easily back-transform the transformed parameter values (as stored in object prm_t) with the same function using the argument rev = TRUE:

prm_t_back <- transformParameters(parm = prm_t, trans = c("log","log"), rev = TRUE)
prm_t_back

##     r     K 
##   0.3 100.0

Indeed, after transformation and back-transformation, the parameter values are back to their original value: $r = 0.3$ and $K = 100$.

Optimization function

We have predefined a function that allows you to do the optimization when supplied with the needed data. The function needs all the inputs that the predictionError function needs (see them described above), but also 2 additional function inputs:

• “calibrateWhich”: a character vector with the names of the parameters that are to be calibrated. This vector needs to be at least of length 2 (otherwise the optim function will not work), but you can choose to only select a subset of the model’s parameters for calibration;
• “transformations”: a character vector with the transformation for each model parameter (see above in the description of function transformParameters).

Optionally, you can supply additional inputs to the optim function via the three dots ellipsis argument ... (e.g. when the calibration does not converge within the default 500 iterations, you can then pass a higher number of iteration as input to optim via argument control, e.g. adding control = list(maxit = 1000) increases the number of iterations to 1000):

# Define a function that performs the calibration of specified parameters
# NOTE: nothing needs to be filled in here, so this function can be used as-is
calibrateParameters <- function(par,  # parameters
                                init, # initial state values
                                fun,  # function that returns the rates of change
                                stateCol, # name of the state variable
                                obs,  # observed values of the state variable
                                t,    # timepoints of the observed state values
                                calibrateWhich, # names of the parameters to calibrate
                                transformations, # transformation per calibration parameter
                                times = seq(from = 0, # times for the numerical integration
                                            to = max(t), # default: 0 - max(t)
                                            by = 1), # in steps of 1
                                ... # NOTE: these ... are ellipsis, so DO NOT fill in something here!
                                    # these are the optional extra arguments past on to optim()
                                ) {
  # check names of parameters that need calibration
  calibrateWhich <- match.arg(calibrateWhich, choices = names(par), several.ok = TRUE)
  
  # Create vector with transformations for all parameters (set to "none")
  tranforms <- rep("none", length(par))
  names(tranforms) <- names(par)
  
  # Overwrite the elements for parameters that need calibration
  tranforms[calibrateWhich] <- transformations
  
  # Transform parameters
  par <- transformParameters(parm = par, trans = tranforms, rev = FALSE)
  
  # Get parameters to be estimated
  par_fit <- par[calibrateWhich]
  
  # Specify the cost function: the function that will be optimized for the parameters to be estimated
  costFunction <- function(parm_optim, parm_all, 
                           init, fun, stateCol, obs, t,
                           times = t,
                           optimNames = calibrateWhich,
                           transf = tranforms,
                           ... # NOTE: these ... are ellipsis, so DO NOT fill in something here!
                           ) {
    # Gather parameters (those that are and are not to be estimated)
    par <- parm_all
    par[optimNames] <- parm_optim[optimNames]
    
    # Back-transform
    par <- transformParameters(parm = par, trans = transf, rev = TRUE)
    
    # Get RMSE
    y <- predRMSE(p = par, y0 = init, fn = fun, 
                  name = stateCol, obs_values = obs, obs_times = t, 
                  out_time = times)

    # Return
    return(y)
  }
  
  # Fit using optim
  fit <- optim(par = par_fit,
               fn =  costFunction,
               parm_all = par, 
               init = init, 
               fun = fun,
               stateCol = stateCol,
               obs = obs,
               t = t,
               times = times,
               ...) # NOTE: these ... are ellipsis, so DO NOT fill in something here!
                    # These are the optional additional arguments passed on to optim()
  
  # Gather parameters: all parameters updated with those that are estimated
  y <- par
  y[calibrateWhich] <- fit$par[calibrateWhich]
  
  # Back-transform
  y <- transformParameters(parm = y, trans = tranforms, rev = TRUE)
  
  # put all (constant and estimated) parameters back in $par
  fit$par <- y
  
  # Add the calibrateWhich information
  fit$calibrated <- calibrateWhich
  
  # return 'fit'object:
  return(fit)
}

We have now defined all the necessary components needed for calibration, so let’s calibrate the model parameters $r$ and $K$ based on the data as defined in surveyData, and check the output:

fit <- calibrateParameters(par = ParmsLogistic(),
                           init = InitLogistic(),
                           fun = RatesLogistic,
                           stateCol = "N",
                           obs = surveyData$N,
                           t = surveyData$t,
                           times = seq(from = 0, to = 50, by = 1),
                           calibrateWhich = c("r","K"),
                           transformations = c("log","log"))
fit

## $par
##          r          K 
##  0.3291443 51.1938939 
## 
## $value
## [1] 6.262373
## 
## $counts
## function gradient 
##       65       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
## 
## $calibrated
## [1] "r" "K"

What the function returns is a list with elements:

• “par”: the calibrated parameter values. This includes all model parameters: those that are kept constant as well as those that have been calibrated. Here the calibrated parameters are approximately $r = 0.33$ and $K = 51.2$;
• “value”: this is the value of the cost-function at the values as given in “par”. Here, this is thus the RMSE of the model predictions given the parameter values in “par”;
• “counts”: information about the number of iterations that the optim function needed;
• “convergence”: information about the convergence of the optim function: a value of 0 means that the estimation successfully converged! For other values, see ?optim for the help-documentation of the optim function;
• “message”: a character string giving any additional information returned by the optimizer, or NULL;
• “calibrated”: a character string with the model parameters that have been calibrated (thus the vector passed as input to argument calibrateWhich).

Using the calibrated parameter values, we can predict the state’s dynamics:

statesFitted <- ode(y = InitLogistic(),
                    times = surveyData$t,
                    parms = fit$par,
                    func = RatesLogistic,
                    method = "ode45")
# Plot the results
plot(surveyData$t, surveyData$N, xlab = "time", ylab = "N")
lines(statesFitted[,"time"], statesFitted[,"N"], col = 2, lwd = 2)
# For clarity, let's add a legend:
legend(x="topleft", bty = "n", lty = c(NA,1), pch = c(1,NA), legend = c("observation","prediction"))

And, we can compute the RMSE for the calibrated model:

predRMSE(p = fit$par,
         y0 = InitLogistic(),
         fn = RatesLogistic,
         name = "N",
         obs_values = surveyData$N,
         obs_times = surveyData$t) # optionally add values for out_time

## [1] 6.262373

Indeed, the model performs much better now! The RMSE decreased from 40.92 to 6.26!

Note that this last step is not a proper evaluation of the calibrated model! Namely, for that we would need an independent dataset.