Skip to contents

Random Effect Covariances

2023-05-30

Source: vignettes/addingCovariances.Rmd
addingCovariances.Rmd

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:

## Load phenobarbitol data
library(nlmixr2)

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))
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00

print(fit)
#> ── nlmixr² SAEM OBJF by FOCEi approximation ──
#> 
#>           OBJF      AIC      BIC Log-likelihood Condition#(Cov) Condition#(Cor)
#> FOCEi 688.6792 985.5502 1003.811      -486.7751        7.570439        6.604927
#> 
#> ── Time (sec $time): ──
#> 
#>            setup optimize covariance   saem table compress
#> elapsed 0.003262    6e-06   0.033009 15.337 3.511    0.029
#> 
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#> 
#>                          Parameter  Est.     SE %RSE    Back-transformed(95%CI)
#> tcl     typical value of clearance -4.99 0.0743 1.49 0.00677 (0.00586, 0.00784)
#> tv         typical value of volume 0.346 0.0538 15.5          1.41 (1.27, 1.57)
#> add.err       residual variability  2.83                                   2.83
#>         BSV(CV%) Shrink(SD)%
#> tcl         52.4      1.50% 
#> tv          41.0      1.09% 
#> add.err                     
#>  
#>   Covariance Type ($covMethod): linFim
#>   Correlations in between subject variability (BSV) matrix:
#>     cor:eta.v,eta.cl 
#>           0.989  
#>  
#> 
#>   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.8 -1.53 -0.422  0.0167 0.337 0.507  17.5 -0.222 -0.0297
#> 2 1      112.  31    29.9  1.15  0.394  0.394  0.653 0.653  27.9  3.11   0.254 
#> 3 2        2    9.7  10.9 -1.18 -0.994 -0.193  0.16  0.423  10.5 -0.813 -0.162 
#> # ℹ 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 <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)
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