nlmixr

nlmixr

This is an example of a complex model that can be estimated.

In the example below, a target-mediated drug disposition PK model for nimotuzumab is illustrated (Rodríguez-Vera et al. 2015).

Model Schematic

Model Schematic

nlmixr model

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.examples"));

## Note to get the datasets you can use:
## remotes::install_github("nlmixrdevelopment/nlmixr.examples")
library(nlmixr.examples) ## For nlmixr.examples data

library(xpose)
#> Loading required package: ggplot2
#> 
#> Attaching package: 'xpose'
#> The following object is masked from 'package:nlmixr':
#> 
#>     vpc
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(xpose.nlmixr)
#> 
#> Attaching package: 'xpose.nlmixr'
#> The following object is masked from 'package:nlmixr':
#> 
#>     vpc
library(ggplot2)

nimo <- function() {
    ini({
        ## Note that the UI can take expressions
        ## Also note that these initial estimates should be provided on the log-scale
        tcl <- log(0.001)
        tv1 <- log(1.45)
        tQ <- log(0.004)
        tv2 <- log(44)
        tkss <- log(12)
        tkint <- log(0.3)
        tksyn <- log(1)
        tkdeg <- log(7)
        ## Initial estimates should be high for SAEM ETAs
        eta.cl  ~ 2
        eta.v1  ~ 2
        eta.kss ~ 2
        ##  Also true for additive error (also ignored in SAEM)
        add.err <- 10
    })
    model({
        cl <- exp(tcl + eta.cl)
        v1 <- exp(tv1 + eta.v1)
        Q  <- exp(tQ)
        v2 <- exp(tv2)
        kss <- exp(tkss + eta.kss)
        kint <- exp(tkint)
        ksyn <- exp(tksyn)
        kdeg <- exp(tkdeg)

        k <- cl/v1
        k12 <- Q/v1
        k21 <- Q/v2

        eff(0) <- ksyn/kdeg ##initializing compartment

        ## Concentration is calculated
        conc = 0.5*(central/v1-eff-kss)+0.5*sqrt((central/v1-eff-kss)**2+4*kss*central/v1)

        d/dt(central)  = -(k+k12)*conc*v1+k21*peripheral-kint*eff*conc*v1/(kss+conc)
        d/dt(peripheral) = k12*conc*v1-k21*peripheral  ##Free Drug second compartment amount
        d/dt(eff) = ksyn - kdeg*eff - (kint-kdeg)*conc*eff/(kss+conc)

        IPRED=log(conc)

        IPRED ~ add(add.err)
    })
}

Fit

fit <- nlmixr(nimo, nimoData, est="saem", control=list(print=0))
#> Loading model already run (/home/matt/R/x86_64-pc-linux-gnu-library/3.5/nlmixr.examples/nlmixr-nimo-nimoData-saem-71a71aacdcb2cc5c08a2b950a446cefd.rds)

Goodness of fit Plots

## Add cwres/npde after fit (can also use nlmixr(...,table=list(cwres=TRUE,npde=TRUE)))
fit  <- fit %>% addCwres() %>% addNpde();
#> Calculating residuals/tables
#> done.
#> Warning in .collectWarnings(do.call(foceiFit.data.frame0, call, envir =
#> parent.frame(1))): Calculated residuals like IPRED are masked by nlmixr
#> calculated values
#> Compiling NPDE model...
#> done


##Goodness-of-fit plots

plot(fit); ## Standard nlmixr plots


################################################################################
## Xpose plots; Need to print otherwise running a script won't
## show xpose plots
################################################################################
xpdb <- xpose_data_nlmixr(fit) ## first convert to nlmixr object

print(dv_vs_pred(xpdb) +
      ylab("Observed Nimotuzumab Concentrations (ug/mL)") +
      xlab("Population Predicted Nimotuzumab Concentrations (ug/mL)"))


print(dv_vs_ipred(xpdb) +
      ylab("Observed Nimotuzumab Concentrations (ug/mL)") +
      xlab("Individual Predicted Nimotuzumab Concentrations (ug/mL)"))


print(res_vs_pred(xpdb) +
      ylab("Conditional Weighted Residuals") +
      xlab("Population Predicted Nimotuzumab Concentrations (ug/mL)"))


print(res_vs_idv(xpdb) +
      ylab("Conditional Weighted Residuals") +
      xlab("Time (h)"))


print(absval_res_vs_idv(xpdb, res = 'IWRES') +
      ylab("Individual Weighted Residuals") +
      xlab("Time (h)"))


print(absval_res_vs_pred(xpdb, res = 'IWRES') +
      ylab("Individual Weighted Residuals") +
      xlab("Population Predicted Nimotuzumab Concentrations (ug/mL)"))


print(ind_plots(xpdb, nrow=3, ncol=4) +
      ylab("Predicted and Observed Nimotuzumab concentrations (ug/mL)") +
      xlab("Time (h)"))


print(res_distrib(xpdb) +
     ylab("Density") +
     xlab("Conditional Weighted Residuals"))


################################################################################
##Visual Predictive Checks
################################################################################
vpc.ui(fit,n=500,stratify=c("dos"), show=list(obs_dv=T),
       bins = c(-0.5,0,25,75,100,200,400,600,750,900,1100,1200,1400,1600,1900,2150,2300),
       ylab = "Nimotuzumab Concentrations (ug/mL)", xlab = "Time (h)")
#> Compiling VPC model...done
#> done (3.22 sec)
#> Warning in data.frame(..., check.names = FALSE): row names were found from
#> a short variable and have been discarded
#>   $rxsim = original simulated data
#>   $sim = merge simulated data
#>   $obs = observed data
#>   $gg = vpc ggplot
#> use vpc(...) to change plot options
#> plotting the object now


vpc.ui(fit,n=500, show=list(obs_dv=T),
       bins = c(-0.5,0,25,75,100,200,400,600,750,900,1100,1200,1400,1600,1900,2150,2300),
       ylab = "Nimotuzumab Concentrations (ug/mL)", xlab = "Time (h)")
#> Compiling VPC model...
#> done
#> done (3.25 sec)
#>   $rxsim = original simulated data
#>   $sim = merge simulated data
#>   $obs = observed data
#>   $gg = vpc ggplot
#> use vpc(...) to change plot options
#> plotting the object now