nlmixr

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 Phenobarbital data:

## Load Phenobarb data
library(nlmixr)

## These options cache the models and the model simulations in R
## To run the actual models on your system, take the save options off.
options(nlmixr.save=TRUE,
        nlmixr.save.dir=system.file(package="nlmixr"));

invisible(memoise::cache_filesystem(file.path(system.file(package="nlmixr"),"memo")))

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")

print(fit)
#> ── nlmixr SAEM(ODE); OBJF by Gaussian Quadrature (n.nodes=3, n.sd=1.6) fit  
#>                OBJF      AIC      BIC Log-likelihood Condition Number
#> gauss3_1.6 721.8749 1018.746 1037.006      -503.3729         7.462652
#> 
#> ── Time (sec; $time): ───────────────────────────────────────────────────── 
#>           saem    setup table covariance logLik    other
#> elapsed 22.104 0.226164 0.019      0.023  0.066 0.281836
#> 
#> ── Population Parameters ($parFixed or $parFixedDf): ────────────────────── 
#>                          Parameter  Est.     SE %RSE
#> tcl     typical value of clearance    -5  0.075  1.5
#> tv         typical value of volume 0.347 0.0536 15.4
#> add.err       residual variability  2.84            
#>           Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
#> tcl     0.00674 (0.00582, 0.0078)     52.8      1.68% 
#> tv              1.42 (1.27, 1.57)     40.8      1.18% 
#> add.err                      2.84                      
#> 
#>   Covariance Type ($covMethod): linFim
#>   Correlations in between subject variability (BSV) matrix:
#>     cor__eta.v.eta.cl 
#>              0.986   
#> 
#>   Full BSV covariance ($omega) or correlation ($omegaR; diagonals=SDs) 
#>   Distribution stats (mean/skewness/kurtosis/p-value) available in $shrink 
#> 
#> ── Fit Data (object is a modified tibble): ──────────────────────────────── 
#> # A tibble: 155 x 16
#>   ID     TIME    DV  EVID  PRED    RES IPRED  IRES  IWRES  eta.cl   eta.v
#>   <fct> <dbl> <dbl> <int> <dbl>  <dbl> <dbl> <dbl>  <dbl>   <dbl>   <dbl>
#> 1 1        2   17.3     0  17.5 -0.198  18.4 -1.13 -0.397 -0.0731 -0.0516
#> 2 1      112.  31       0  27.9  3.08   29.6  1.39  0.489 -0.0731 -0.0516
#> 3 2        2    9.7     0  10.5 -0.799  12.5 -2.76 -0.973 -0.215  -0.171 
#> # … with 152 more rows, and 5 more variables: cl <dbl>, v <dbl>, ke <dbl>,
#> #   cp <dbl>, A1 <dbl>

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

Two types of VPCs

library(ggplot2)
## A traditional VPC
p1 <- vpc(fit, show=list(obs_dv=TRUE)) + ylab("Concentrations")

## A prediction-corrected VPC
p2 <- vpc(fit, pred_corr = TRUE, show=list(obs_dv=TRUE)) + ylab("Prediction-Corrected Concentrations")

library(gridExtra)
grid.arrange(p1,p2)