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 algoritm takes something small).

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"),
  break_prof = FALSE,
  ...
)

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
break_prof	if 'TRUE' (default), stop the profiling if a better optimum is located
...	additional parameters passed on to stats::optim() and calibrate()

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)

# exposure - control run
exp <- Schmitt2013 %>%
  filter(ID == "T0") %>%
  select(time=t, conc)

# observations - control run
obs <- Schmitt2013 %>%
  filter(ID == "T0") %>%
   select(t, BM=obs)

# update metsulfuron
myscenario <- metsulfuron %>%
  set_param(c(k_phot_fix = TRUE,Emax = 1)) %>%
  set_init(c(BM = 12)) %>%
  set_exposure(exp)

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

# Likelihood profiling

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

cvasi documentation built on Sept. 23, 2024, 9:08 a.m.