Fit aggregate data using Iterative Reweighting with Monte Carlo updates

timedIRMC implements the Iterative Reweighting (IRMC) algorithm for parameter estimation of aggregate data models, iterating over maximum likelihood updates with weighted Monte Carlo updates. This function is used to compare the performance of different implementations of aggregate data modeling.

Usage

timedIRMC(init, opts, obs, maxiter = 100, convcrit_nll = 5e-04, 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.
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 5e-04.
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 Iterative Reweighting algorithm with Monte Carlo sampling for optimization. At each iteration, it generates Monte Carlo samples and updates the parameter estimates using weighted importance sampling. 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 <- timedIRMC(opts$pt, opts, examplomycin_aggregated)
#> iteration 1, nll=58.6588992213695
#> iteration 2, nll=10120.0652052777
#> iteration 3, nll=976.217989450599
#> iteration 4, nll=1902.74575614608
#> iteration 5, nll=4.37702663109096
#> iteration 6, nll=-1281.31566755092
#> iteration 7, nll=-1737.5996228288
#> iteration 8, nll=-1834.58029289539
#> iteration 9, nll=-1845.11490679444
#> iteration 10, nll=-1845.11667765342
#> iteration 11, nll=-1845.15311166459
#> iteration 12, nll=-1845.13996335289
#> iteration 13, nll=-1845.14145074914
#> iteration 14, nll=-1845.14150897451
#> should break now due to no difference between OFV and appr OFV
print(result)
#> # A tibble: 14 × 5
#>    p               nll appr_nll time                 iter
#>    <list>        <dbl>    <dbl> <dttm>              <dbl>
#>  1 <dbl [11]>    58.7      58.7 2025-11-28 09:47:21     1
#>  2 <dbl [11]> 10120.    -1498.  2025-11-28 09:47:23     2
#>  3 <dbl [11]>   976.     -773.  2025-11-28 09:47:24     3
#>  4 <dbl [11]>  1903.    -1679.  2025-11-28 09:47:30     4
#>  5 <dbl [11]>     4.38  -1780.  2025-11-28 09:47:33     5
#>  6 <dbl [11]> -1281.    -1806.  2025-11-28 09:47:35     6
#>  7 <dbl [11]> -1738.    -1825.  2025-11-28 09:47:36     7
#>  8 <dbl [11]> -1835.    -1841.  2025-11-28 09:47:37     8
#>  9 <dbl [11]> -1845.    -1845.  2025-11-28 09:47:38     9
#> 10 <dbl [11]> -1845.    -1845.  2025-11-28 09:47:40    10
#> 11 <dbl [11]> -1845.    -1845.  2025-11-28 09:47:40    11
#> 12 <dbl [11]> -1845.    -1845.  2025-11-28 09:47:41    12
#> 13 <dbl [11]> -1845.    -1845.  2025-11-28 09:47:41    13
#> 14 <dbl [11]> -1845.    -1845.  2025-11-28 09:47:41    14