Compute expected population-level mean and covariance

gen_pop_EV computes the expected population-level mean and covariance of the model predictions using either first-order approximation (FO) or Monte Carlo (MC) methods. The function is a core component of the aggregate data modeling framework, used to generate expected values for model fitting and optimization.

Usage

gen_pop_EV(opts)

Arguments

opts

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

• Prediction function (f)

• Time points for prediction

• Population parameters (beta)

• Between-subject variability (Omega)

• Method selection (fo_appr)

• Number of Monte Carlo samples (nsim)

• Omega expansion factor

Value

A list containing:

• E: Expected mean of the model predictions at each time point

• V: Expected covariance matrix of the predictions across time points

Details

The function implements two methods for computing population expectations:

1. First-Order (FO) Approximation (when opts$fo_appr is TRUE):

   • Uses linearization around the random effects (η)

   • Steps: a. Generate random effects samples b. Compute importance sampling weights based on η distribution c. For each individual:

     • Compute FOCE approximation of mean and covariance

     • Apply importance sampling weights d. Combine weighted results across individuals

   • Advantages:

     • Computationally efficient for complex models

     • Good accuracy for nearly linear systems

   • Used when:

     • Number of Monte Carlo samples is small (nsim < 10)

     • Model is approximately linear in random effects

2. Monte Carlo (MC) Approximation (when opts$fo_appr is FALSE):

   • Uses direct simulation to estimate expectations

   • Steps: a. Generate random effects samples b. Compute full model predictions for each sample c. Calculate empirical mean and covariance d. Apply importance sampling if omega expansion is used

   • Advantages:

     • No linearization assumptions

     • More accurate for highly nonlinear models

   • Used when:

     • Sufficient Monte Carlo samples available

     • High accuracy required

     • Model is highly nonlinear

The choice between methods depends on:

• Model complexity and nonlinearity

• Required accuracy

• Computational resources

• Number of available Monte Carlo samples

Examples

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


# 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 (30% CV)
    Sigma_prop = 0.04
  ),
  nsim = 2500,  # Number of Monte Carlo samples
  n = 2500,      # Number of individuals
  fo_appr = FALSE  # Use Monte Carlo approximation
)

# Generate population expectations
ev <- gen_pop_EV(opts)