Generate a fitting function for optimization

genfitfunc generates a fitting function for optimization that computes the negative log-likelihood of the model given the observed data and current parameter estimates. The function handles both raw data and aggregate data (mean and covariance) formats.

Usage

genfitfunc(opts, obs)

Arguments

opts

A list of model options generated by genopts(). Contains settings for the model, including:

• f: Prediction function for the model

• time: Vector of observation times

• p: List of parameter values and structure

• nsim: Number of Monte Carlo samples

• n: Number of individuals

• pt: Parameter values on transformed scale

• ptrans: Parameter transformation function

obs

Observed data in one of two formats: 1. Aggregate form: List with elements: - E: Vector of means for each time point - V: Covariance matrix of observations 2. Raw data matrix: - Rows: Individual observations (nsim) - Columns: Time points

Value

A function with signature function(p, givedetails = FALSE, opts_overrides) where:

• p: Parameter values (transformed or untransformed)

• givedetails: If TRUE, returns additional attributes: - EV: Expected values (mean and covariance) - obs: Observed data - nllfun: Negative log-likelihood function - opts: Model options used

• opts_overrides: Optional list to override model options

• Returns: Negative log-likelihood value

Details

Algorithm Steps:

• 1. Data Preparation:

     • Validates input data format

     • Converts raw data to aggregate form if needed

     • Uses expected data if observations not provided

  2. Parameter Processing:

     • Handles both transformed and untransformed parameters

     • Updates model options with current parameters

  3. Likelihood Computation:

     • Generates expected values using gen_pop_EV

     • Computes covariance matrix inverse

     • Calculates negative log-likelihood

  4. Error Handling:

     • Checks for matrix inversion problems

     • Validates data dimensions

     • Ensures proper parameter transformations

Mathematical Details:

• Negative Log-Likelihood: $$-2\log L = n\log|V| + (y - \mu)^T V^{-1} (y - \mu)$$ where: - n: Number of individuals - V: Model-predicted covariance matrix - y: Observed means - μ: Model-predicted means

• Parameter Transformations: - Log transform for positive parameters - Logit transform for bounded parameters - Identity for unrestricted parameters

Examples

# Load required libraries
library(admr)
library(rxode2)
library(nlmixr2)
library(dplyr)
library(tidyr)
library(mnorm)


# Load and prepare data
data(examplomycin)
examplomycin_wide <- examplomycin %>%
  filter(EVID != 101) %>%
  dplyr::select(ID, TIME, DV) %>%
  pivot_wider(names_from = TIME, values_from = DV) %>%
  dplyr::select(-c(1))

# Create aggregated data
examplomycin_aggregated <- examplomycin_wide %>%
  admr::meancov()

# Define RxODE model
rxModel <- function(){
model({
  # Parameters
  ke = cl / v1             # Elimination rate constant
  k12 = q / v1             # Rate constant for central to peripheral transfer
  k21 = q / v2             # Rate constant for peripheral to central transfer

  # Differential equations
  d/dt(depot)    = -ka * depot
  d/dt(central)  = ka * depot - ke * central - k12 * central + k21 * peripheral
  d/dt(peripheral) = k12 * central - k21 * peripheral

  # Concentration in central compartment
  cp = central / v1
})
}

rxModel <- rxode2(rxModel)
#>  
#>  
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
rxModel <- rxModel$simulationModel
#>  
#>  

# Define prediction function for a two-compartment model
predder <- function(time, theta_i, dose = 100) {
  n_individuals <- nrow(theta_i)
  if (is.null(n_individuals)) n_individuals <- 1

  # Create event table for dosing and sampling
  ev <- eventTable(amount.units="mg", time.units="hours")
  ev$add.dosing(dose = dose, nbr.doses = 1, start.time = 0)
  ev$add.sampling(time)

  # Solve ODE system
  out <- rxSolve(rxModel, params = theta_i, events = ev, cores = 0)

  # Return matrix of predictions
  matrix(out$cp, nrow = n_individuals, ncol = length(time), byrow = TRUE)
}

# Create options for a two-compartment model
opts <- genopts(
  f = predder,
  time = c(.1, .25, .5, 1, 2, 3, 5, 8, 12),
  p = list(
    # Population parameters (fixed effects)
    beta = c(cl = 5,    # Clearance (L/h)
            v1 = 10,    # Central volume (L)
            v2 = 30,    # Peripheral volume (L)
            q = 10,     # Inter-compartmental clearance (L/h)
            ka = 1),    # Absorption rate (1/h)

    # Between-subject variability (30% CV on all parameters)
    Omega = matrix(c(0.09, 0, 0, 0, 0,
                    0, 0.09, 0, 0, 0,
                    0, 0, 0.09, 0, 0,
                    0, 0, 0, 0.09, 0,
                    0, 0, 0, 0, 0.09), nrow = 5, ncol = 5),

    # Residual error (20% CV)
    Sigma_prop = 0.04
  ),
  nsim = 2500,  # Number of Monte Carlo samples
  n = 500,      # Number of individuals
  fo_appr = FALSE  # Use Monte Carlo approximation
)

# Generate objective function for optimization
objfun <- genfitfunc(opts, examplomycin_aggregated)

# Test the objective function with initial parameters
init_params <- opts$p
nll <- objfun(init_params)