lik_profile: Likelihood profiling

In cvasi: Calibration, Validation, and Simulation of TKTD Models

Likelihood profiling

Description

The aim of the function is two-fold: 1) estimate a 95% confidence around each parameter of a calibrated model, and 2) see if perhaps a local minimum was found rather than a global minimum. To achieve this, the likelihood profiling goes through every parameter one by one. For each parameter, the model is sequentially refit with the parameter value set to increasingly lower and higher values, and the likelihood of the model given the data calculated (using log_lik()). The likelihood is then compared to the likelihood of the original model (using a likelihood ratio). This leads to the development of a likelihood profile, from which a plot a 95% confidence interval for the parameter is derived.

The idea of the function is a variable stepwise algorithm: When the likelihood ratio changes very little (less than l_crit_min), the stepsize is increased (up to a maximum, specified by f_step_max). When the lik. ratio changes too much (more than l_crit_max), the algorithm tries again with a smaller stepsize (also bound to a minimum: f_step_min). Note that the stepsize is used as a fraction of the parameter value that is tried. To prevent very small stepsizes when the value goes towards zero (as can be the case for effect thresholds), an absolute minimum stepsize (f_step_abs), which is specified as a fraction of the best parameter value (Xhat) (unless it is zero, then the algorithm takes something small).

Note that the likelihood of the model given the data can be calculated across all datasets provided in the calibration set x, or calculated separately for each individual dataset before being combined into one likelihood (by adjusting the optional parameter individual). The latter has the advantage that different datasets can be given different weights in the likelihood calculation (using the "weight" slot of the caliset objects,x). Further, for continuous data (e.g. biomass), the likelihood considers the variance (standard deviation) in the log likelihood calculation, which can vary between datasets when the likelihood is calculated for each dataset separately before combining into an overall likelihood. The latter could be relevant when factors might lead to variability between datasets (e.g. different labs, different animal culture,...)

To conduct the likelihood calculations on separate datasets, the parameter individual which by default is 'FALSE' can be set to 'TRUE'. Then, then log likelihoods are calculated for each dataset individually (or in subgroups, using the "tag" names of the caliset object, if provided, to group datasets with the same "tag" before calculating the log likelihood). Subsequently, the log likelihoods for the subsets are combined into an overall likelihood (considering the set weights provided in the "weight" slot of the caliset object). Note that for each set only 1 weight can be provided (i.e. not individual weights for each datapoint within the set), and that set with the same tag should have identical weight.

The function was inspired by a MatLab BYOM v.6.8 procedure, created by Tjalling Jager. For details, please refer to BYOM (http://debtox.info/byom.html) as well as Jager (2021).

Usage

lik_profile(
  x,
  par,
  output,
  data = NULL,
  bounds = NULL,
  refit = TRUE,
  type = c("coarse", "fine"),
  individual = FALSE,
  break_prof = FALSE,
  log_scale = FALSE,
  data_type = c("continuous", "count"),
  ...
)

Arguments

x	either a single scenario or a list of caliset objects
par	named vector - parameters (names and values) to be profiled
output	character vector, name of output column of simulate() that is used in calibration
data	only needed if x is a scenario
bounds	optional list of lists (including lower and upper bound): uses defaults in x object, but can be overwritten here (e.g. bounds <- list(k_resp = list(0,10), k_phot_max = list(0,30)) )
refit	if TRUE (default), refit if a better minimum is found
type	"fine" or "coarse" (default) likelihood profiling
individual	if FALSE (default), the log likelihood is calculated across the whole dataset. Alternatively, if TRUE, log likelihoods are calculated for each (group of) set(s) individually.
break_prof	If TRUE, then stop the profiling if a better optimum is located. Default is FALSE.
log_scale	FALSE (default), option to calculate the log likelihood on a log scale (i.e., observations and predictions are log transformed during calculation)
data_type	Character argument, "continuous" (default) or "count", to specify the data type for the log likelihood calculations.
...	additional parameters passed on to calibrate() and simulate(). To avoid parameter confusion, use argument method to select optimization algorithms of calibrate() and argument ode_method to select numerical integration schemes of package deSolve.

Value

A list containing, for each parameter profiled, the likelihood profiling results as a dataframe; the 95% confidence interval; the original parameter value; the likelihood plot object; and the recalibrated parameter values (in case a lower optimum was found)

References

Jager T, 2021: Robust Likelihood-Based Optimization and Uncertainty Analysis of Toxicokinetic-Toxicodynamic Models. Integrated Environmental Assessment and Management 17:388-397. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1002/ieam.4333")}

Examples

# Example with Lemna model - physiological params
library(dplyr)

# observations - control run
obs <- schmitt2013 %>%
  filter(trial == "T0")

# update metsulfuron
myscenario <- metsulfuron %>%
  set_param(c(k_phot_fix = TRUE, Emax = 1)) %>%
  set_init(c(BM = 0.0012)) %>%
  set_noexposure() %>%
  set_bounds(list(k_phot_max=c(0, 1)))

fit <- calibrate(
  x = myscenario,
  par = c(k_phot_max = 1),
  data = obs,
  output = "FrondNo",
  method = "Brent"
)

# Likelihood profiling

res <- lik_profile(
  x = myscenario,
  data = obs,
  output = "FrondNo",
  par = fit$par,
  bounds = list(
    k_phot_max = list(0, 30)
  ),
  refit = FALSE,
  type = "fine",
  method = "Brent"
)
# plot
plot(res)

cvasi documentation built on Sept. 11, 2025, 5:11 p.m.