log_lik: Calculate log likelihood

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

Calculate log likelihood

Description

Calculates the sum of log likelihoods of each observation given the model parameterization.

Current implementations enable log likelihood calculations for:

1. continuous data, considering a normal distribution around the prediction for each datapoint,

2. count data, considering a multinomial distribution for data reporting the number of survivors over time.

The log likelihood calculation for count data 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

log_lik(
  obs,
  pred,
  data_type = c("continuous", "count"),
  log_scale = FALSE,
  npars = NULL
)

Arguments

obs	numeric vector of observed values
pred	numeric vector of predicted values
data_type	determines the if likelihood profiling is conducted for "continuous" (default) or "count" data
log_scale	FALSE (default), option to calculate the log likelihood on a log scale (i.e., observations and predictions are log transformed during calculation)
npars	named numeric vector of parameters that the model was calibrated on, required for "continuous" data type, optional for "count".

Value

the log likelihood value

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

# simple example for continuous data #####
# observations
obs <- c(12, 38, 92, 176, 176, 627, 1283, 2640)
# intercept, a, and slope, b, of a Poisson regression fitted through obs
pars <- c(a = 2, b = 0.73)
# predictions with the Poisson regression
pred <- c(15.43, 32.15, 66.99, 139.57, 290.82, 605.94, 1262.52, 2630.58)
# example plot
plot(seq(1:length(obs)), obs)
lines(seq(1:length(obs)), pred)
log_lik(
  obs = obs,
  pred = pred,
  npars = length(pars),
)

# example with count data and GUTS model #####
library(dplyr)
# observational data
dt <- ringtest_c %>% filter(replicate == "E")
myexposure <- dt %>% select(time, conc)
obs <- dt %>%
  mutate(S=Nsurv / max(Nsurv)) %>%
  select(time, S)
# GUTS model
GUTS_RED_IT() %>%
  set_param(c(hb = 0)) %>%
  set_exposure(myexposure) -> myscenario
# fit
fit <- calibrate(
  x = myscenario,
  par = c(kd=1.2, alpha=9.2, beta=4.3),
  data = obs,
  output = "S")
# update
myscenario <- myscenario %>% set_param(fit$par)
# simulate
pred <- myscenario %>% simulate()
pred <- pred$S #* max(obs$S)
obs <- obs$S
# calc likelihood
log_lik(obs,
    pred,
    data_type = "count")

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