Estimate the posterior distribution of individual parameters by MCMC
Source:R/param_estim.R
poso_estim_mcmc.RdEstimates the posterior distribution of individual parameters by Markov Chain Monte Carlo (using a Metropolis-Hastings algorithm)
Arguments
- dat
Dataframe. An individual subject dataset following the structure of NONMEM/rxode2 event records.
- prior_model
A
posologyrprior population pharmacokinetics model, a list of six objects.- return_model
A boolean. Returns a rxode2 model using the estimated ETAs if set to
TRUE.- burn_in
Number of burn-in iterations for the Metropolis-Hastings algorithm.
- n_iter
Total number of iterations (following the burn-in iterations) for each Markov chain of the Metropolis-Hastings algorithm.
- n_chains
Number of Markov chains
- nocb
A boolean. for time-varying covariates: the next observation carried backward (nocb) interpolation style, similar to NONMEM. If
FALSE, the last observation carried forward (locf) style will be used. Defaults toFALSE.- control
A list of parameters controlling the Metropolis-Hastings algorithm.
Value
If return_model is set to FALSE, a list of one element: a
dataframe $eta of ETAs from the posterior distribution, estimated by
Markov Chain Monte Carlo.
If return_model is set to TRUE, a list of the dataframe of the posterior
distribution of ETA, and a rxode2 model using the estimated distributions of
ETAs.
References
Comets E, Lavenu A, Lavielle M. Parameter estimation in nonlinear mixed effect models using saemix, an R implementation of the SAEM algorithm. Journal of Statistical Software 80, 3 (2017), 1-41.
Examples
# model
mod_run001 <- function() {
ini({
THETA_Cl <- 4.0
THETA_Vc <- 70.0
THETA_Ka <- 1.0
ETA_Cl ~ 0.2
ETA_Vc ~ 0.2
ETA_Ka ~ 0.2
prop.sd <- sqrt(0.05)
})
model({
TVCl <- THETA_Cl
TVVc <- THETA_Vc
TVKa <- THETA_Ka
Cl <- TVCl*exp(ETA_Cl)
Vc <- TVVc*exp(ETA_Vc)
Ka <- TVKa*exp(ETA_Ka)
K20 <- Cl/Vc
Cc <- centr/Vc
d/dt(depot) = -Ka*depot
d/dt(centr) = Ka*depot - K20*centr
Cc ~ prop(prop.sd)
})
}
# df_patient01: event table for Patient01, following a 30 minutes intravenous
# infusion
df_patient01 <- data.frame(ID=1,
TIME=c(0.0,1.0,14.0),
DV=c(NA,25.0,5.5),
AMT=c(2000,0,0),
EVID=c(1,0,0),
DUR=c(0.5,NA,NA))
# estimate the posterior distribution of population parameters
poso_estim_mcmc(dat=df_patient01,prior_model=mod_run001,
n_iter=50,n_chains=2)
#>
#>
#>
#>
#> $eta
#> ETA_Cl ETA_Vc ETA_Ka
#> 1 0.40129734 -0.221078172 0.930425874
#> 2 0.46219626 -0.408861774 0.919479583
#> 3 0.46733199 -0.303257501 0.231586070
#> 4 0.70605689 -0.389072743 0.231586070
#> 5 0.71226862 -0.460682597 0.220690179
#> 6 0.45848168 -0.156826054 0.655923830
#> 7 0.71578116 -0.356073917 0.443066290
#> 8 0.39429314 -0.250227766 0.029447638
#> 9 0.79934008 -0.300123505 0.029447638
#> 10 0.64891523 -0.602467643 -0.257596645
#> 11 0.68534557 -0.393584606 0.471011853
#> 12 0.41497705 -0.469426296 0.200708751
#> 13 -0.28560298 -0.336731063 0.574758130
#> 14 0.69972327 -0.384404332 -0.094838457
#> 15 0.56697115 -0.427648209 -0.154894296
#> 16 0.52583412 -0.193851708 0.029564358
#> 17 0.73647568 -0.145748104 0.240492264
#> 18 0.72223572 -0.174057317 0.765403502
#> 19 0.61003356 -0.349027845 0.290769918
#> 20 0.69904470 -0.047834356 0.422258849
#> 21 0.43859699 -0.612412417 -0.098503394
#> 22 0.69096316 -0.051944909 0.914823122
#> 23 0.54729672 0.005972018 0.337291445
#> 24 0.27761493 -0.237833012 0.271761486
#> 25 0.64326906 0.011240352 0.623424475
#> 26 0.61126409 -0.364932526 -0.290915542
#> 27 0.36850878 -0.351410700 0.388243848
#> 28 0.27661006 -0.581498663 -0.230915305
#> 29 -0.23732985 -0.526465862 0.407033222
#> 30 0.19616325 -0.645638552 -0.494479393
#> 31 0.47330453 -0.034944362 0.470257171
#> 32 0.33639400 -0.515573239 -0.050276895
#> 33 0.78830319 -0.387354276 -0.031447965
#> 34 0.48064119 -0.343826881 0.244052299
#> 35 0.60189998 0.034441740 0.469481783
#> 36 0.77837261 -0.215748962 0.185718893
#> 37 0.75914818 -0.176417195 0.302207225
#> 38 0.68090244 -0.346820687 0.036191074
#> 39 0.54490321 -0.629963624 0.473465262
#> 40 0.49828937 -0.283363262 0.356249587
#> 41 0.54257518 -0.494014971 0.374892911
#> 42 0.62410466 -0.260718197 0.409061732
#> 43 0.15855684 -0.878222982 -0.131889899
#> 44 0.75140183 -0.141857343 0.326987451
#> 45 0.07179731 -0.032165060 0.268610671
#> 46 0.41210387 -0.258656564 0.462389655
#> 47 0.74462060 -0.382517639 0.673070802
#> 48 0.78287745 -0.073199344 0.521618599
#> 49 0.73952962 0.231836592 0.907839064
#> 50 0.53905726 -0.042012025 0.483275077
#> 51 0.00000000 0.000000000 0.000000000
#> 52 0.42534184 -0.403708105 1.056572451
#> 53 0.55618191 -0.264477715 0.046032939
#> 54 0.45520633 -0.679464662 -0.016242431
#> 55 0.44009565 -0.409873762 0.340896578
#> 56 0.57429692 -0.387726602 0.613547280
#> 57 0.54488941 -0.287829000 0.140149483
#> 58 0.48402318 -0.798975605 0.161344579
#> 59 0.56660040 -0.479785996 0.138614462
#> 60 0.61758410 -0.297867848 0.278458135
#> 61 0.73780864 0.023043168 0.551507628
#> 62 0.68002525 -0.319894513 0.234239231
#> 63 0.89503107 -0.037166313 0.600292939
#> 64 0.78780453 -0.263742002 0.175122402
#> 65 0.51360658 0.052946634 0.892954576
#> 66 0.64420451 -0.318996298 0.650988522
#> 67 0.68813422 -0.251065357 0.415733557
#> 68 0.78325406 0.094584973 0.439275927
#> 69 0.48416155 -0.480703520 0.124244451
#> 70 0.48645250 -0.663751410 0.007279203
#> 71 0.42371624 -0.511424673 0.239185177
#> 72 0.79998492 -0.327110133 -0.287008745
#> 73 0.55606168 -0.300343887 -0.103813643
#> 74 -0.18397797 -0.957175830 0.033660836
#> 75 0.33043554 -0.828096146 -0.054257139
#> 76 0.18819311 -0.839629729 0.143333722
#> 77 0.08262417 -0.414383051 -0.225163689
#> 78 0.37230513 -0.375817501 0.666112298
#> 79 0.51699274 -0.934285571 -0.299695700
#> 80 0.08783468 -0.800477076 -0.147109306
#> 81 0.68133590 -0.368815872 0.040840687
#> 82 0.65031284 -0.323466120 0.398707403
#> 83 0.52147855 -0.586734782 -0.127577976
#> 84 0.71718247 -0.171249178 0.372917693
#> 85 0.52524631 -0.377795316 0.640937425
#> 86 0.80589852 0.156365752 0.993932092
#> 87 0.43878864 -0.746336825 0.356184103
#> 88 0.13803478 -0.812773685 0.083051161
#> 89 0.49241332 -0.249973773 0.507466133
#> 90 0.81523875 -0.116686814 0.532386442
#> 91 0.28613945 -0.434134521 -0.055910865
#> 92 0.16400899 -0.611578714 -0.265540778
#> 93 0.49782289 -0.500626081 0.300302428
#> 94 0.69449874 -0.069173401 0.542000876
#> 95 0.44901310 -0.630662806 -0.542320029
#> 96 0.39487244 0.142613375 1.402640547
#> 97 0.64429477 -0.461460431 0.136638918
#> 98 0.57768550 -0.319792195 0.097921632
#> 99 0.56412639 -0.496782084 -0.041628981
#> 100 0.32575969 -0.106694594 0.762696118
#>
#> $model
#> ── Solved rxode2 object ──
#> ── Parameters ($params): ──
#> # A tibble: 100 × 8
#> sim.id THETA_Cl THETA_Vc THETA_Ka prop.sd ETA_Cl ETA_Vc ETA_Ka
#> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 4 70 1 0.224 0.401 -0.221 0.930
#> 2 2 4 70 1 0.224 0.462 -0.409 0.919
#> 3 3 4 70 1 0.224 0.467 -0.303 0.232
#> 4 4 4 70 1 0.224 0.706 -0.389 0.232
#> 5 5 4 70 1 0.224 0.712 -0.461 0.221
#> 6 6 4 70 1 0.224 0.458 -0.157 0.656
#> 7 7 4 70 1 0.224 0.716 -0.356 0.443
#> 8 8 4 70 1 0.224 0.394 -0.250 0.0294
#> 9 9 4 70 1 0.224 0.799 -0.300 0.0294
#> 10 10 4 70 1 0.224 0.649 -0.602 -0.258
#> # ℹ 90 more rows
#> ── Initial Conditions ($inits): ──
#> depot centr AUC
#> 0 0 0
#>
#> Simulation without uncertainty in parameters, omega, or sigma matricies
#>
#> ── First part of data (object): ──
#> # A tibble: 200 × 14
#> sim.id time TVCl TVVc TVKa Cl Vc Ka K20 rxCc Cc depot
#> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 1 4 70 1 5.98 56.1 2.54 0.106 28.4 28.4 3.19e+2
#> 2 1 14 4 70 1 5.98 56.1 2.54 0.106 8.61 8.61 7.99e-9
#> 3 2 1 4 70 1 6.35 46.5 2.51 0.137 33.7 33.7 3.25e+2
#> 4 2 14 4 70 1 6.35 46.5 2.51 0.137 6.96 6.96 -2.81e-9
#> 5 3 1 4 70 1 6.38 51.7 1.26 0.123 22.2 22.2 7.90e+2
#> 6 3 14 4 70 1 6.38 51.7 1.26 0.123 7.85 7.85 6.03e-5
#> # ℹ 194 more rows
#> # ℹ 2 more variables: centr <dbl>, AUC <dbl>
#>