Fit aggregate data using Iterative Reweighting with Monte Carlo updates
Source:R/fitIRMC.R
fitIRMC.RdfitIRMC 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 uses nloptr instead of optimx for optimization and includes additional
features like multiple chains and phase-based optimization.
Usage
fitIRMC(
opts,
obs,
maxiter = 100,
convcrit_nll = 1e-05,
single_dataframe = TRUE,
phase_fractions = c(0.2, 0.4, 0.2, 0.2),
max_worse_iterations = 10,
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.
- phase_fractions
Vector of fractions for each optimization phase. The sum should be 1. Default is c(0.2, 0.4, 0.2, 0.2).
- max_worse_iterations
Maximum number of consecutive worse iterations before skipping a phase. Default is 10.
- 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 iterationnll: Negative log-likelihood valuestime: Computation time for each iterationiter: Iteration numberchain: Chain number (if multiple chains are used)phase: Optimization phase number
Details
The function uses the Iterative Reweighting algorithm with Monte Carlo sampling for optimization. It includes several advanced features:
Multiple optimization phases with different convergence criteria
Chains with perturbed starting values
Phase-based optimization with automatic phase skipping
Convergence checking based on both likelihood and parameter stationarity
The optimization process is divided into phases, each with its own convergence criteria and settings. The algorithm can automatically skip phases if the optimization is not progressing.
Examples
# Load required libraries
library(admr)
library(rxode2)
#> rxode2 4.1.1 using 2 threads (see ?getRxThreads)
#> no cache: create with `rxCreateCache()`
library(nlmixr2)
#> ── Attaching packages ───────────────────────────────────────── nlmixr2 4.0.1 ──
#> ✔ lotri 1.0.2 ✔ nlmixr2extra 3.0.2
#> ✔ nlmixr2data 2.0.9 ✔ nlmixr2plot 3.0.3
#> ✔ nlmixr2est 4.1.1
#> ── Optional Packages Loaded/Ignored ─────────────────────────── nlmixr2 4.0.1 ──
#> ✖ babelmixr2 ✖ nonmem2rx
#> ✖ ggPMX ✖ posologyr
#> ✖ monolix2rx ✖ shinyMixR
#> ✖ nlmixr2lib ✖ xpose.nlmixr2
#> ✖ nlmixr2rpt
#> ── Conflicts ───────────────────────────────────────────── nlmixr2conflicts() ──
#> ✖ nlmixr2est::boxCox() masks rxode2::boxCox()
#> ✖ nlmixr2est::yeoJohnson() masks rxode2::yeoJohnson()
library(dplyr)
#>
#> Attaching package: ‘dplyr’
#> The following objects are masked from ‘package:stats’:
#>
#> filter, lag
#> The following objects are masked from ‘package:base’:
#>
#> intersect, setdiff, setequal, union
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,
f = predder
)
# Run optimization
result <- fitIRMC(opts, examplomycin_aggregated)
#> Chain 1:
#> Iter | NLL and Parameters (11 values)
#> --------------------------------------------------------------------------------
#> 1: 70.389 1.609 2.303 3.401 2.303 0.000 -2.408 -2.408 -2.408 -2.408 -2.408 -3.219
#>
#> ### Wide Search Phase ###
#> 2: 8980.099 1.592 4.303 2.776 3.112 2.000 -2.622 -1.782 -4.408 -3.384 -0.808 -2.831
#> 3: 2355.402 -0.408 2.303 2.738 2.995 0.000 -0.622 -2.060 -6.408 -5.012 -1.449 -2.096
#> 4: -1415.828 1.592 2.329 2.709 3.023 0.216 -2.622 -3.146 -5.235 -4.964 -2.509 -2.955
#> 5: -562.311 2.104 3.067 2.535 3.172 0.085 -0.622 -1.146 -7.235 -5.383 -1.219 -2.988
#> 6: -1593.680 1.541 2.335 2.465 3.198 0.025 -1.596 -3.146 -9.235 -4.394 -2.213 -3.035
#> 7: -1736.426 1.662 2.405 2.451 3.258 0.102 -2.199 -5.146 -9.255 -2.394 -1.363 -3.062
#> 8: -1561.368 1.578 2.348 2.428 3.081 -0.027 -1.742 -5.972 -10.276 -0.802 -1.462 -3.117
#> 9: -1772.619 1.644 2.321 2.420 3.325 0.047 -2.355 -6.735 -11.640 -1.977 -1.622 -3.131
#> 10: -1774.674 1.643 2.317 2.416 3.321 0.036 -2.298 -7.909 -11.692 -1.984 -1.624 -3.128
#> 11: -1770.642 1.632 2.376 2.417 3.344 0.105 -2.301 -6.470 -11.889 -2.083 -1.613 -3.128
#> 12: -1770.794 1.623 2.360 2.421 3.365 0.083 -2.355 -7.696 -11.958 -1.870 -1.715 -3.142
#> 13: -1769.597 1.624 2.351 2.425 3.357 0.081 -2.351 -8.703 -11.890 -1.806 -1.720 -3.134
#> 14: -1770.377 1.642 2.348 2.421 3.315 0.062 -2.268 -9.747 -11.914 -1.795 -1.566 -3.132
#> 15: -1767.866 1.619 2.334 2.419 3.305 0.036 -2.262 -8.551 -11.787 -1.806 -1.587 -3.134
#> 16: -1763.517 1.615 2.392 2.420 3.356 0.141 -2.342 -7.155 -11.733 -1.977 -1.585 -3.144
#> 17: -1765.093 1.645 2.388 2.414 3.272 0.092 -2.213 -8.509 -11.947 -1.602 -1.511 -3.146
#> 18: -1777.086 1.640 2.384 2.412 3.308 0.084 -2.364 -9.038 -12.193 -2.040 -1.609 -3.128
#> 19: -1776.439 1.640 2.381 2.409 3.297 0.080 -2.267 -10.058 -12.221 -1.886 -1.554 -3.135
#> 20: -1760.158 1.619 2.406 2.410 3.345 0.168 -2.280 -8.863 -12.065 -2.034 -1.525 -3.142
#>
#> ### Focussed Search Phase ###
#> 21: -1766.286 1.626 2.378 2.408 3.292 0.060 -2.158 -9.798 -12.527 -1.632 -1.509 -3.133
#> 22: -1779.490 1.638 2.376 2.407 3.312 0.076 -2.369 -10.053 -12.760 -2.028 -1.612 -3.133
#> 23: -1780.168 1.636 2.374 2.405 3.308 0.077 -2.366 -10.656 -12.834 -1.998 -1.610 -3.130
#> 24: -1780.697 1.636 2.372 2.404 3.306 0.071 -2.363 -11.304 -12.906 -2.061 -1.564 -3.137
#> 25: -1781.029 1.638 2.371 2.401 3.302 0.066 -2.300 -11.937 -12.993 -1.981 -1.574 -3.132
#> 26: -1782.014 1.638 2.370 2.400 3.298 0.062 -2.338 -12.552 -13.055 -1.988 -1.580 -3.136
#> 27: -1783.109 1.637 2.369 2.398 3.309 0.068 -2.342 -13.049 -13.083 -2.032 -1.595 -3.136
#> 28: -1783.642 1.637 2.368 2.397 3.308 0.071 -2.340 -13.725 -13.162 -2.064 -1.596 -3.136
#> 29: -1784.110 1.637 2.368 2.396 3.308 0.070 -2.339 -14.334 -13.237 -2.063 -1.629 -3.135
#> 30: -1784.720 1.638 2.368 2.394 3.306 0.065 -2.345 -14.934 -13.313 -2.064 -1.600 -3.138
#> 31: -1785.148 1.636 2.368 2.393 3.308 0.068 -2.348 -15.544 -13.392 -2.098 -1.600 -3.140
#> 32: -1784.573 1.636 2.368 2.395 3.306 0.067 -2.344 -16.189 -13.473 -2.033 -1.600 -3.140
#> 33: -1784.967 1.636 2.368 2.394 3.306 0.067 -2.374 -16.755 -13.543 -2.063 -1.602 -3.141
#> 34: -1785.241 1.637 2.367 2.393 3.308 0.067 -2.374 -17.362 -13.619 -2.095 -1.603 -3.142
#> 35: -1785.491 1.637 2.367 2.392 3.308 0.067 -2.374 -17.970 -13.694 -2.062 -1.602 -3.143
#> 36: -1785.725 1.637 2.367 2.392 3.308 0.066 -2.374 -18.579 -13.770 -2.030 -1.603 -3.143
#> 37: -1785.813 1.637 2.367 2.392 3.308 0.066 -2.374 -19.151 -13.841 -2.098 -1.603 -3.139
#> 38: -1785.844 1.637 2.367 2.392 3.308 0.066 -2.374 -19.760 -13.916 -2.065 -1.602 -3.139
#> 39: -1785.912 1.637 2.367 2.391 3.311 0.066 -2.374 -20.331 -13.987 -2.065 -1.602 -3.143
#> 40: -1785.960 1.637 2.367 2.391 3.311 0.067 -2.337 -20.901 -14.051 -2.033 -1.602 -3.143
#> 41: -1785.978 1.637 2.367 2.391 3.311 0.066 -2.337 -21.509 -14.126 -2.066 -1.602 -3.143
#> 42: -1786.004 1.637 2.367 2.391 3.311 0.066 -2.367 -22.078 -14.196 -2.036 -1.602 -3.143
#> 43: -1786.022 1.637 2.367 2.391 3.312 0.066 -2.367 -22.686 -14.272 -2.069 -1.602 -3.143
#> 44: -1786.034 1.637 2.367 2.391 3.312 0.066 -2.367 -23.294 -14.347 -2.036 -1.602 -3.143
#> 45: -1786.022 1.637 2.367 2.391 3.312 0.066 -2.367 -23.903 -14.423 -2.069 -1.602 -3.143
#> 46: -1786.041 1.637 2.367 2.391 3.312 0.066 -2.367 -24.511 -14.498 -2.036 -1.602 -3.143
#> 47: -1786.002 1.637 2.367 2.391 3.315 0.066 -2.367 -25.043 -14.498 -2.071 -1.602 -3.143
#> 48: -1786.025 1.637 2.367 2.391 3.315 0.066 -2.332 -25.577 -14.465 -2.040 -1.602 -3.143
#> 49: -1786.025 1.637 2.367 2.391 3.315 0.066 -2.332 -26.077 -14.465 -2.040 -1.602 -3.143
#> 50: -1785.972 1.637 2.367 2.391 3.317 0.066 -2.366 -26.503 -14.465 -2.078 -1.602 -3.143
#> 51: -1785.974 1.637 2.367 2.391 3.317 0.067 -2.366 -27.003 -14.465 -2.078 -1.602 -3.143
#> 52: -1786.032 1.637 2.367 2.391 3.317 0.067 -2.366 -27.003 -15.036 -2.007 -1.602 -3.143
#> 53: -1786.062 1.637 2.367 2.391 3.317 0.063 -2.363 -27.003 -15.568 -2.039 -1.602 -3.143
#> 54: -1786.071 1.637 2.367 2.391 3.317 0.063 -2.363 -27.003 -16.068 -2.039 -1.602 -3.143
#> 55: -1786.071 1.637 2.367 2.391 3.317 0.063 -2.363 -27.003 -16.568 -2.039 -1.602 -3.143
#> 56: -1786.071 1.637 2.367 2.391 3.317 0.063 -2.363 -27.003 -17.068 -2.039 -1.602 -3.143
#> 57: -1786.080 1.636 2.367 2.391 3.317 0.063 -2.363 -27.003 -17.567 -2.039 -1.602 -3.143
#> 58: -1786.088 1.636 2.367 2.391 3.317 0.063 -2.363 -27.003 -18.067 -2.039 -1.602 -3.143
#> 59: -1786.088 1.636 2.367 2.391 3.317 0.063 -2.363 -27.003 -18.567 -2.039 -1.602 -3.143
#> 60: -1786.094 1.636 2.367 2.391 3.317 0.063 -2.363 -27.003 -19.067 -2.039 -1.602 -3.143
#>
#> ### Fine-Tuning Phase ###
#> 61: -1786.078 1.636 2.367 2.391 3.317 0.063 -2.363 -27.236 -19.067 -2.056 -1.602 -3.143
#> 62: -1785.916 1.636 2.367 2.391 3.317 0.063 -2.363 -27.447 -19.049 -2.057 -1.622 -3.143
#> 63: -1785.880 1.636 2.367 2.391 3.317 0.063 -2.363 -27.662 -19.032 -2.075 -1.622 -3.143
#> 64: -1785.968 1.636 2.367 2.391 3.317 0.064 -2.363 -27.662 -19.337 -2.038 -1.622 -3.126
#> 65: -1786.066 1.636 2.367 2.391 3.317 0.064 -2.363 -27.662 -19.569 -2.038 -1.604 -3.126
#> 66: -1786.066 1.636 2.367 2.391 3.317 0.064 -2.363 -27.662 -19.819 -2.038 -1.604 -3.126
#> 67: -1786.066 1.636 2.367 2.391 3.317 0.064 -2.363 -27.662 -20.069 -2.038 -1.604 -3.126
#> 68: -1786.066 1.636 2.367 2.391 3.317 0.064 -2.363 -27.662 -20.319 -2.038 -1.604 -3.126
#> 69: -1786.063 1.636 2.367 2.391 3.317 0.062 -2.363 -27.662 -20.570 -2.038 -1.604 -3.126
#> 70: -1786.094 1.636 2.367 2.391 3.316 0.063 -2.362 -27.662 -20.841 -2.036 -1.604 -3.142
#>
#> ### Precision Phase ###
#> 71: -1786.114 1.635 2.367 2.391 3.316 0.064 -2.362 -27.013 -19.067 -2.039 -1.602 -3.143
#> 72: -1786.123 1.634 2.367 2.391 3.317 0.065 -2.363 -27.023 -19.067 -2.039 -1.602 -3.143
#> 73: -1786.113 1.635 2.367 2.391 3.317 0.065 -2.363 -27.033 -19.067 -2.039 -1.603 -3.143
#> 74: -1786.124 1.634 2.367 2.391 3.317 0.066 -2.363 -27.042 -19.068 -2.039 -1.603 -3.143
#> 75: -1786.114 1.635 2.367 2.391 3.316 0.067 -2.363 -27.052 -19.068 -2.039 -1.603 -3.143
#> 76: -1786.121 1.634 2.367 2.391 3.315 0.066 -2.363 -27.062 -19.068 -2.040 -1.603 -3.143
#> 77: -1786.112 1.635 2.367 2.391 3.316 0.067 -2.363 -27.071 -19.068 -2.040 -1.603 -3.143
#> 78: -1786.120 1.634 2.367 2.391 3.315 0.066 -2.363 -27.081 -19.068 -2.040 -1.604 -3.143
#> 79: -1786.111 1.635 2.367 2.391 3.316 0.067 -2.363 -27.091 -19.068 -2.040 -1.604 -3.143
#> 80: -1786.118 1.634 2.367 2.391 3.315 0.066 -2.363 -27.100 -19.068 -2.040 -1.604 -3.143
#> 81: -1786.108 1.635 2.367 2.391 3.316 0.067 -2.363 -27.110 -19.068 -2.041 -1.604 -3.143
#> 82: -1786.116 1.634 2.367 2.391 3.315 0.066 -2.363 -27.120 -19.068 -2.041 -1.604 -3.143
#> 83: -1786.108 1.635 2.367 2.391 3.316 0.067 -2.363 -27.130 -19.068 -2.041 -1.604 -3.143
#> 84: -1786.116 1.634 2.367 2.391 3.315 0.066 -2.363 -27.139 -19.068 -2.041 -1.604 -3.143
#>
#> Chain 1 Complete: Final NLL = -1786.124, Time Elapsed = 90.18 seconds
#>
print(result)
#> -- FitIRMC Summary --
#>
#> -- Objective Function and Information Criteria --
#> Log-likelihood: -1786.1235
#> AIC: 3583.25
#> BIC: 3640.61
#> Condition#(Cov): 29.39
#> Condition#(Cor): 35.43
#>
#> -- Timing Information --
#> Best Chain: 90.1788 seconds
#> All Chains: 90.1811 seconds
#> Covariance: 13.0886 seconds
#> Elapsed: 103.27 seconds
#>
#> -- Population Parameters --
#> # A tibble: 6 × 6
#> Parameter Est. SE `%RSE` `Back-transformed(95%CI)` `BSV(CV%)`
#> <chr> <dbl> <dbl> <dbl> <chr> <dbl>
#> 1 cl 1.63 0.0141 0.861 5.12 (4.98, 5.27) 30.7
#> 2 v1 2.37 0.0380 1.61 10.67 (9.90, 11.50) 0.000134
#> 3 v2 2.39 0.0206 0.863 10.92 (10.49, 11.38) 0.00724
#> 4 q 3.32 0.0282 0.849 27.57 (26.09, 29.13) 36.1
#> 5 ka 0.0660 0.0375 56.8 1.07 (0.99, 1.15) 44.9
#> 6 Residual Error 0.0432 NA NA 0.0432 NA
#>
#> -- Iteration Diagnostics --
#> Iter | NLL and Parameters
#> --------------------------------------------------------------------------------
#> 1: 70.389 1.609 2.303 3.401 2.303 0.000 -2.408 -2.408 -2.408 -2.408 -2.408 -3.219
#> 2: 8980.099 1.592 4.303 2.776 3.112 2.000 -2.622 -1.782 -4.408 -3.384 -0.808 -2.831
#> 3: 2355.402 -0.408 2.303 2.738 2.995 0.000 -0.622 -2.060 -6.408 -5.012 -1.449 -2.096
#> 4: -1415.828 1.592 2.329 2.709 3.023 0.216 -2.622 -3.146 -5.235 -4.964 -2.509 -2.955
#> 5: -562.311 2.104 3.067 2.535 3.172 0.085 -0.622 -1.146 -7.235 -5.383 -1.219 -2.988
#> ... (omitted iterations) ...
#> 80: -1786.118 1.634 2.367 2.391 3.315 0.066 -2.363 -27.100 -19.068 -2.040 -1.604 -3.143
#> 81: -1786.108 1.635 2.367 2.391 3.316 0.067 -2.363 -27.110 -19.068 -2.041 -1.604 -3.143
#> 82: -1786.116 1.634 2.367 2.391 3.315 0.066 -2.363 -27.120 -19.068 -2.041 -1.604 -3.143
#> 83: -1786.108 1.635 2.367 2.391 3.316 0.067 -2.363 -27.130 -19.068 -2.041 -1.604 -3.143
#> 84: -1786.116 1.634 2.367 2.391 3.315 0.066 -2.363 -27.139 -19.068 -2.041 -1.604 -3.143
#>