Fit aggregate data using the bobyqa algorithm with Monte Carlo updates

fitMC 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

fitMC(
  opts,
  obs,
  maxiter = 5000,
  convcrit_nll = 1e-05,
  single_dataframe = TRUE,
  chains = 1,
  perturbation = 0.1,
  seed = 1,
  use_grad = FALSE
)

Arguments

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.

maxiter

Maximum number of iterations for the optimization algorithm. Default is 100.

convcrit_nll

Convergence criterion for the negative log-likelihood. The algorithm stops when the relative change in negative log-likelihood is less than this value. Default is 1e-05.

single_dataframe

Logical indicating whether to use a single data frame (TRUE) or multiple data frames (FALSE). Default is TRUE.

chains

Number of chains to run. Default is 1.

perturbation

Perturbation factor for the initial parameter values of each chain. Default is 0.1.

seed

Random seed for reproducibility. Default is 1.

use_grad

Logical indicating whether to use gradient information in optimization.

Value

An object of class fit_admr_result containing:

• p: List of parameter estimates for each iteration

• nll: Negative log-likelihood values

• time: Computation time for each iteration

• iter: Iteration number

• chain: Chain number (not used in the standard MC algorithm)

• phase: Optimization phase number

Details

The function uses the Monte Carlo algorithm with Monte Carlo sampling for optimization. The optimization process is performed using the bobyqa algorithm, which is suitable for non-linear optimization problems. The function supports both single and multiple data frames, allowing flexibility in how the observed data is structured.

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
#>  
#>  
#> using C compiler: ‘gcc (Ubuntu 13.3.0-6ubuntu2~24.04) 13.3.0’

# 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 <- fitMC(opts, examplomycin_aggregated)
#> Error: Not a matrix.
print(result)
#> Error: object 'result' not found