nlmixr
Adding Covariances between random effects
You can simply add co-variances between two random effects by adding
the effects together in the model specification block, that is
eta.cl+eta.v ~. After that statement, you specify the lower
triangular matrix of the fit with c().
An example of this is the phenobarbitol data:
Model Specification
pheno <- function() {
ini({
tcl <- log(0.008) # typical value of clearance
tv <- log(0.6) # typical value of volume
## var(eta.cl)
eta.cl + eta.v ~ c(1,
0.01, 1) ## cov(eta.cl, eta.v), var(eta.v)
# interindividual variability on clearance and volume
add.err <- 0.1 # residual variability
})
model({
cl <- exp(tcl + eta.cl) # individual value of clearance
v <- exp(tv + eta.v) # individual value of volume
ke <- cl / v # elimination rate constant
d/dt(A1) = - ke * A1 # model differential equation
cp = A1 / v # concentration in plasma
cp ~ add(add.err) # define error model
})
}Fit with SAEM
fit := nlmixr(pheno, pheno_sd, "saem",
control=list(print=0),
table=list(cwres=TRUE, npde=TRUE))
print(fit)
#> ── nlmixr² SAEM OBJF by FOCEi approximation ──
#>
#> OBJF AIC BIC Log-likelihood Condition#(Cov) Condition#(Cor)
#> FOCEi 688.7367 985.6076 1003.868 -486.8038 20.63833 19.53368
#>
#> ── Time (sec $time): ──
#>
#> setup optimize covariance preprocess configure saem postprocess
#> elapsed 0.08438001 3.4232e-05 0.02600548 0.033 1.171 9.861 0.012
#> table compress
#> elapsed 3.617 0.084
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#>
#> Est. SE %RSE Back-transformed(95%CI) BSV(CV%)
#> tcl -5.00 0.0663 1.33 0.00671 (0.00589, 0.00764) 51.8
#> tv 0.347 0.0530 15.3 1.41 (1.27, 1.57) 41.6
#> add.err 2.78 6.95e-310 2.50e-308 2.78 (2.78, 2.78)
#> Shrink(SD)%
#> tcl 2.57
#> tv 1.30
#> add.err
#>
#> Covariance Type ($covMethod): sa
#> Some strong fixed parameter correlations exist ($cor) :
#> cor:tv,tcl cor:om.eta.cl,add.err
#> 0.903 -0.0986
#> cor:cov.eta.v.eta.cl,add.err cor:om.eta.v,add.err
#> 0.0345 -0.0645
#> cor:cov.eta.v.eta.cl,om.eta.cl cor:om.eta.v,om.eta.cl
#> 0.851 0.509
#> cor:om.eta.v,cov.eta.v.eta.cl
#> 0.816
#>
#>
#> Correlations in between subject variability (BSV) matrix:
#> cor:eta.v,eta.cl
#> 0.962
#>
#>
#> Full BSV covariance ($omega) or correlation ($omegaR; diagonals=SDs)
#> Distribution stats (mean/skewness/kurtosis/p-value) available in $shrink
#> Censoring ($censInformation): No censoring
#>
#> ── Fit Data (object is a modified tibble): ──
#> # A tibble: 155 × 26
#> ID TIME DV EPRED ERES NPDE NPD PDE PD PRED RES WRES
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 2 17.3 18.9 -1.62 -0.332 -0.100 0.37 0.46 17.5 -0.210 -0.0279
#> 2 1 112. 31 29.7 1.31 0.253 0.253 0.6 0.6 28.0 3.04 0.249
#> 3 2 2 9.7 11.4 -1.75 -0.664 -0.245 0.253 0.403 10.5 -0.806 -0.160
#> # ℹ 152 more rows
#> # ℹ 14 more variables: IPRED <dbl>, IRES <dbl>, IWRES <dbl>, CPRED <dbl>,
#> # CRES <dbl>, CWRES <dbl>, eta.cl <dbl>, eta.v <dbl>, A1 <dbl>, cl <dbl>,
#> # v <dbl>, ke <dbl>, tad <dbl>, dosenum <int>Basic Goodness of Fit Plots
plot(fit)
Those individual plots are not that great, it would be better to see
the actual curves; You can with augPred
plot(augPred(fit))
Two types of VPCs
library(ggplot2)
p1 <- vpcPlot(fit, show=list(obs_dv=TRUE));
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
p1 <- p1 + ylab("Concentrations")
## A prediction-corrected VPC
p2 <- vpcPlot(fit, pred_corr = TRUE, show=list(obs_dv=TRUE))
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
p2 <- p2 + ylab("Prediction-Corrected Concentrations")
library(patchwork)
p1 / p2