Estimation: nlmixr2

Introduction

While the PKPDsim R package does not include functionality to perform population estimation within the package, it does offer a convenient translator to be able to use the nlmixr2 R package for parameter estimation. nlmixr2 allows fitting of population PK-PD models using common algorithms applied in pharmacometrics, such as FOCE and SAEM. For installation and usage of nlmixr2, please refer to the documentation in the nlmixr2 package and on the respective website.

Translation

Under the hood, nlmixr2 uses the rxode2 R package to perform PKPD model simulations, and which is very similar to PKPDsim. The model syntax for both rxode2 and nlmixr2 compared to PKPDsim are also very similar, making translation fairly straightforward. A model translator function has therefore been included in PKPDsim that allows translation of a PKPDsim model to nlmixr2 syntax:

## Error in get(paste0(generic, ".", class), envir = get_method_env()) : 
##   object 'type_sum.accel' not found

mod <- new_ode_model("pk_1cmt_iv")
f <- pkpdsim_to_nlmixr(
  model = mod,
  parameters = list(CL = 5, V = 50),
  omega = c(0.1, 0.05, 0.1),
  res_var = list(prop = 0.1, add = 0.1),
  log_transform = TRUE
)

The returned object f generated by pkpdsim_to_nlmixr() is an object that defines the required model and parameter definitions in nlmixr syntax:

f

## function () 
## {
##     ini({
##         logCL <- log(5)
##         logV <- log(50)
##         err_prop <- c(0, 0.1, 1)
##         err_add <- c(0, 1, 10)
##         eta_CL + eta_V ~ c(0.1, 0.05, 0.1)
##     })
##     model({
##         CL <- exp(logCL + eta_CL)
##         V <- exp(logV + eta_V)
##         d/dt(A1) = -(CL/V) * A1
##         y = A1/V
##         y ~ prop(err_prop) + add(err_add)
##     })
## }
## <environment: 0x55d9d8762048>

A full code example is included below.

Note: both nlmixr and the nlmixr-translator in PKPDsim are still under active development. Syntax and results may therefore change.

Example

library(tidyverse)
library(nlmixr2)

## Define parameters
par   <- list(CL = 5, V = 50)
ruv   <- list(prop = 0.1, add = 1)
omega <- c(0.1,
           0.05, 0.1)

## Simulate data
t_obs <- c(8, 23.5, 25, 71.5)
n <- 50
regimen <- new_regimen(
  amt = 1500, 
  n = 4,
  interval = 24
)
conc <- sim(
  ode = mod, 
  regimen = regimen, 
  parameters = par,
  omega = omega, 
  only_obs = TRUE, 
  t_obs = t_obs, 
  n = n,
  res_var = ruv
) %>%
  mutate(EVID = 0, AMT = 0, MDV = 0, CMT = 1) %>%
  select(ID = id, TIME = t, CMT, DV = y, AMT, EVID, MDV)
doses <- regimen_to_nm(
  regimen, 
  n_ind = n
)

## Combine observed data and dose data
simdat <- bind_rows(doses, conc) %>% 
  arrange(ID, TIME)

## Create nlmixr model object
f <- pkpdsim_to_nlmixr(
  model = mod,
  parameters = par,
  omega = c(0.1, 0.05, 0.1),
  res_var = list(prop = 0.1, add = 0.1),
  log_transform = TRUE
)

## Perform fit using nlmixr (SAEM)
fit <- nlmixr(
  object = f,
  data = simdat,
  est = "saem"
)

fit

## ── nlmixr² SAEM OBJF by FOCEi approximation ──
## 
##  Gaussian/Laplacian Likelihoods: AIC(fit) or fit$objf etc. 
##  FOCEi CWRES & Likelihoods: addCwres(fit) 
## 
## ── Time (sec fit$time): ──
## 
##            setup covariance   saem table compress   other
## elapsed 0.002713   0.017007 11.102 0.077    0.024 1.77728
## 
## ── Population Parameters (fit$parFixed or fit$parFixedDf): ──
## 
##           Est.     SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
## logCL     1.57 0.0399 2.54        4.81 (4.45, 5.2)     20.6      25.7% 
## logV      3.92 0.0548  1.4       50.2 (45.1, 55.9)     27.7      26.9% 
## err_prop 0.212                               0.212                     
## err_add    1.5                                 1.5                     
##  
##   Covariance Type (fit$covMethod): linFim
##   Correlations in between subject variability (BSV) matrix:
##     cor:eta_V,eta_CL 
##           0.689  
##  
## 
##   Full BSV covariance (fit$omega) or correlation (fit$omegaR; diagonals=SDs) 
##   Distribution stats (mean/skewness/kurtosis/p-value) available in fit$shrink 
##   Information about run found (fit$runInfo):
##    • 'f' has the following user-defined boundaries: err_prop, err_add which are ignored in 'saem' 
##   Censoring (fit$censInformation): No censoring
## 
## ── Fit Data (object fit is a modified tibble): ──
## # A tibble: 200 × 16
##   ID     TIME    DV  PRED   RES IPRED  IRES IWRES eta_CL  eta_V     y    A1
##   <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>  <dbl> <dbl> <dbl>
## 1 1       8   19.7  13.9   5.84 17.1  2.62  0.669 -0.239 -0.127 17.1   756.
## 2 1      23.5  4.85  3.15  1.70  4.54 0.310 0.174 -0.239 -0.127  4.54  201.
## 3 1      25   38.0  29.9   8.16 35.1  2.90  0.383 -0.239 -0.127 35.1  1554.
## # ℹ 197 more rows
## # ℹ 4 more variables: CL <dbl>, V <dbl>, tad <dbl>, dosenum <dbl>