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:
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: 0x55dd00641f90>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 = T,
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) %>%
dplyr::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 = T
)
## 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.001745 0.016009 9.159 0.063 0.022 1.415246
##
## ── Population Parameters (fit$parFixed or fit$parFixedDf): ──
##
## Est. SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
## logCL 1.57 0.0382 2.44 4.79 (4.44, 5.16) 17.1 31.6%
## logV 3.92 0.055 1.4 50.6 (45.4, 56.4) 25.2 31.6%
## err_prop 0.23 0.23
## err_add 1.62 1.62
##
## Covariance Type (fit$covMethod): linFim
## Correlations in between subject variability (BSV) matrix:
## cor:eta_V,eta_CL
## 0.761
##
##
## Full BSV covariance (fit$omega) or correlation (fit$omegaR; diagonals=SDs)
## Distribution stats (mean/skewness/kurtosis/p-value) available in fit$shrink
## 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.82 16.0 3.67 0.911 -0.156 -0.109 16.0 728.
## 2 1 23.5 4.85 3.20 1.64 3.96 0.889 0.478 -0.156 -0.109 3.96 180.
## 3 1 25 38.0 29.7 8.28 33.7 4.38 0.553 -0.156 -0.109 33.7 1528.
## # ℹ 197 more rows
## # ℹ 4 more variables: CL <dbl>, V <dbl>, tad <dbl>, dosenum <dbl>