Fit aggregate data using the bobyqa algorithm with Monte Carlo updates

timedbobyqa implements the bobyqa algorithm for parameter estimation of aggregate data models, iterating over maximum likelihood updates with Monte Carlo updates. Each iteration creates new Monte Carlo samples and updates the parameter values. This function is used to compare the performance of different implementations of aggregate data modeling.

Usage

timedbobyqa(init, opts, obs, nomap = TRUE)

Arguments

init: Initial parameter values for optimization. These should be transformed parameters as generated by opts$pt.
opts: A list of model options generated by genopts(). Contains settings for the model, including the prediction function, time points, parameter structure, and simulation settings.
obs: Observed data in aggregate form (mean and covariance) or as a matrix of raw data.
nomap: Logical indicating whether to use multiple models (FALSE) or a single model (TRUE). Default is TRUE.

Value

A data frame containing:

p: List of parameter estimates for each iteration
nll: Negative log-likelihood values
time: Computation time for each iteration
iter: Iteration number

Details

The function uses the bobyqa algorithm from the optimx package for optimization. It performs Monte Carlo sampling at each iteration to update the parameter estimates. The algorithm continues until convergence or until the maximum number of iterations is reached.

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
predder <- function(time, theta_i, dose = 100) {
  n_individuals <- nrow(theta_i)
  if (is.null(n_individuals)) n_individuals <- 1

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

  out <- rxSolve(rxModel, params = theta_i, events = ev, cores = 0)
  cp_matrix <- matrix(out$cp, nrow = n_individuals, ncol = length(time),
                      byrow = TRUE)
  return(cp_matrix)
}

# Create options
opts <- genopts(
  time = c(.1, .25, .5, 1, 2, 3, 5, 8, 12),
  p = list(
    beta = c(cl = 5, v1 = 10, v2 = 30, q = 10, ka = 1),
    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),
    Sigma_prop = 0.04
  ),
  nsim = 2500,
  n = 500,
  fo_appr = FALSE,
  omega_expansion = 1.2,
  f = predder
)

# Run optimization
result <- timedbobyqa(opts$pt, opts, examplomycin_aggregated)
#> Error: Not a matrix.
print(result)
#> Error: object 'result' not found