nlmixr2/rxode2 mu referencing 3.0/4.0
By Matthew Fidler in rxode2 nlmixr2
December 11, 2024
This month, I will talk about about a new iteration of mu-referencing in nlmixr2, which I call mu3 and mu4.
What is mu referencing in nlmixr2 – Review from another post
From the last blog post about mu-referencing, I will give a brief overview of mu-referencing and what mu-2 referencing is and how it is expanded a bit more.
mu-referencing is combining a fixed effect, random effect and possibly a covariate in the form:
\[ \theta_\mathsf{pop}+\eta_\mathsf{individual}+\theta_\mathsf{covariate}\times \mathsf{DataCovariate} \]
Often they are placed in exponentials for these to be log-normally distributed like:
\[ \exp\left(\theta_\mathsf{pop}+\eta_\mathsf{individual}+\theta_\mathsf{covariate}\times \mathsf{DataCovariate}\right) \]
In optimization routines like saem, these are switched out with a single parameter during optimization classically called \(\phi\) in both NONMEM and Monolix.
Once the best \(\phi\) is found, then the population, individual and covariates can be found by linear regression of the individual \(\phi\) values versus the information in the optimization.
This linear model adds important stability and speed when determining these parameters in the mu-expression. (It also adds rules like they must not be time-varying for instance).
In old versions of rxode2 and nlmixr2 the \(\mathsf{DataCovariate}\) had to be in the dataset itself. This classic weight covariate adjustment:
\[ \exp\left(\theta_{Cl}+\eta_{Cl}\right)\times\left(\frac{WT}{70}\right)^{3/4} \]
would have to written:
\[ \exp\left(\theta_{Cl}+\eta_{Cl} + 3/4\times\log\left(\frac{WT}{70}\right)\right) \]
If you wanted to estimate the population parameter \(3/4\) to see if it approaches the correct value you could with:
\[ \exp\left(\theta_{Cl}+\eta_{Cl} + \theta_{Cl, \textsf{cov}}\times\mathsf{DataCovariate}\right) \]
Where
\[ \mathsf{DataCovariate} =\log\left(\frac{WT}{70}\right) \]
This is easy enough to do and adds stabilization.
However, with mu referencing 2.0 you can simply use an additive expression to setup mu-referencing:
\[ \exp\left(\theta_{Cl}+\eta_{Cl} + \frac{WT}{70}+\theta_{Cl, \textsf{cov}}\times\log\left(WT/70\right)\right) \]
This is a bit more convenient than creating a column in the dataset that does this conversion and less user-based intervention to make nlmixr2 use linear models when it can.
Checking for mu2 referencing
In rxode2 / nlmixr2 you can check to see if your version of nlmixr2 supports the mu2 referencing by evaluating the functional form:
one.compartment <- function() {
ini({
tka <- log(1.57); label("Ka")
tcl <- log(2.72); label("Cl")
tv <- log(31.5); label("V")
wt.cl <- 0.75; label("WT on CL")
eta.ka ~ 0.6
eta.cl ~ 0.3
eta.v ~ 0.1
add.sd <- 0.7
})
# and a model block with the error specification and model specification
model({
ka <- exp(tka + eta.ka)
cl <- exp(tcl + eta.cl + wt.cl*log(WT/70))
v <- exp(tv + eta.v)
d/dt(depot) <- -ka * depot
d/dt(center) <- ka * depot - cl / v * center
cp <- center / v
cp ~ add(add.sd)
})
}
one <- one.compartment()
print(one)## ── rxode2-based free-form 2-cmt ODE model ──────────────────────────────────────
## ── Initalization: ──
## Fixed Effects ($theta):
## tka tcl tv wt.cl add.sd
## 0.4510756 1.0006319 3.4499875 0.7500000 0.7000000
##
## Omega ($omega):
## eta.ka eta.cl eta.v
## eta.ka 0.6 0.0 0.0
## eta.cl 0.0 0.3 0.0
## eta.v 0.0 0.0 0.1
##
## States ($state or $stateDf):
## Compartment Number Compartment Name
## 1 1 depot
## 2 2 center
## ── μ-referencing ($muRefTable): ──
## theta eta level covariates
## 1 tka eta.ka id
## 2 tcl eta.cl id log(0.0142857142857143 * WT)*wt.cl
## 3 tv eta.v id
##
## ── Model (Normalized Syntax): ──
## function() {
## ini({
## tka <- 0.451075619360217
## label("Ka")
## tcl <- 1.00063188030791
## label("Cl")
## tv <- 3.44998754583159
## label("V")
## wt.cl <- 0.75
## label("WT on CL")
## add.sd <- c(0, 0.7)
## eta.ka ~ 0.6
## eta.cl ~ 0.3
## eta.v ~ 0.1
## })
## model({
## ka <- exp(tka + eta.ka)
## cl <- exp(tcl + eta.cl + wt.cl * log(WT/70))
## v <- exp(tv + eta.v)
## d/dt(depot) <- -ka * depot
## d/dt(center) <- ka * depot - cl/v * center
## cp <- center/v
## cp ~ add(add.sd)
## })
## }If the mu2 referencing is supported it will show the following mu reference table:
── μ-referencing ($muRefTable): ──
theta eta level covariates
1 tka eta.ka id
2 tcl eta.cl id log(0.0142857142857143 * WT)*wt.cl
3 tv eta.v idHere it shows the transformation that is used when creating the transformed data used for mu-based covariate modeling. In this case, we have log(0.0142857142857143 * WT). It is a bit different than what is written because it is prepossessed by symengine and looks at the derivative with respect to the covariate parameter wt.cl.
In general, I am excited by this new feature in nlmixr2 because it adds a new level of simplicity to user-based models and more often detects mu-referenced code when it may not have been detected in the past.
What is new with mu referencing 3.0
The mu referencing 3.0 allows you to specify the mu referencing anywhere that makes sense in the model. For example:
one.compartment <- function() {
ini({
tka <- log(1.57); label("Ka")
tcl <- log(2.72); label("Cl")
tv <- log(31.5); label("V")
wt.cl <- 0.75; label("WT on CL")
eta.ka ~ 0.6
eta.cl ~ 0.3
eta.v ~ 0.1
add.sd <- 0.7
})
# and a model block with the error specification and model specification
model({
wt70 <- log(WT/70)
ka <- exp(tka + eta.ka)
cl <- exp(tcl + eta.cl + wt.cl*wt70)
v <- exp(tv + eta.v)
d/dt(depot) <- -ka * depot
d/dt(center) <- ka * depot - cl / v * center
cp <- center / v
cp ~ add(add.sd)
})
}
one <- one.compartment()
f <- nlmixr2(one, theo_sd, "saem", control=saemControl(print=0))## ℹ loading model to look for μ₃/μ₄ references## → loading into symengine environment...## → pruning branches (`if`/`else`) of saem model...## ✔ done## ℹ done## ℹ μ₃ item: wt70## → loading into symengine environment...## → pruning branches (`if`/`else`) of saem model...## ✔ done## → finding duplicate expressions in saem model...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → optimizing duplicate expressions in saem model...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## ✔ done## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’## ℹ calculate uninformed etas## ℹ done## Calculating covariance matrix## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → loading into symengine environment...## → pruning branches (`if`/`else`) of saem model...## ✔ done## → finding duplicate expressions in saem predOnly model 0...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → finding duplicate expressions in saem predOnly model 1...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → optimizing duplicate expressions in saem predOnly model 1...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → finding duplicate expressions in saem predOnly model 2...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## ✔ done## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’## → Calculating residuals/tables## ✔ done## → compress origData in nlmixr2 object, save 7040## → compress phiM in nlmixr2 object, save 63104## → compress parHistData in nlmixr2 object, save 14840## → compress saem0 in nlmixr2 object, save 30440print(f)## ── nlmixr² SAEM OBJF by FOCEi approximation ──
##
## Gaussian/Laplacian Likelihoods: AIC() or $objf etc.
## FOCEi CWRES & Likelihoods: addCwres()
##
## ── Time (sec $time): ──
##
## setup covariance saem table compress other
## elapsed 0.001598 0.011004 5.108 0.057 0.021 2.379398
##
## ── Population Parameters ($parFixed or $parFixedDf): ──
##
## Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
## tka Ka 0.462 0.194 42 1.59 (1.09, 2.32) 70.8 0.188%
## tcl Cl 1.02 0.0832 8.16 2.77 (2.35, 3.26) 26.5 5.12%
## tv V 3.45 0.0447 1.29 31.7 (29, 34.6) 13.1 9.33%
## wt.cl WT on CL 0.564 0.62 110 0.564 (-0.65, 1.78)
## add.sd 0.696 0.696
##
## Covariance Type ($covMethod): linFim
## No correlations in between subject variability (BSV) matrix
## Full BSV covariance ($omega) or correlation ($omegaR; diagonals=SDs)
## Distribution stats (mean/skewness/kurtosis/p-value) available in $shrink
## Information about run found ($runInfo):
## • μ₃ item: wt70
## Censoring ($censInformation): No censoring
##
## ── Fit Data (object is a modified tibble): ──
## # A tibble: 132 × 21
## ID TIME DV PRED RES IPRED IRES IWRES eta.ka eta.cl eta.v cp
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 0 0.74 0 0.74 0 0.74 1.06 0.0875 -0.549 -0.0849 0
## 2 1 0.25 2.84 3.27 -0.430 3.84 -0.998 -1.43 0.0875 -0.549 -0.0849 3.84
## 3 1 0.57 6.57 5.84 0.733 6.77 -0.202 -0.290 0.0875 -0.549 -0.0849 6.77
## # ℹ 129 more rows
## # ℹ 9 more variables: depot <dbl>, center <dbl>, wt70 <dbl>, ka <dbl>,
## # cl <dbl>, v <dbl>, tad <dbl>, dosenum <dbl>, WT <dbl>Note that the run info shows that nlmixr2 shows that the mu-3 referencing was detected:
f$runInfo## [1] "μ₃ item: wt70"This is only calculated when running nlmixr2
What is new with mu referencing 4.0
Another new feature in nlmixr2 is assigning string variables. Internally, rxode2 converts this to a numeric and the converts to a factor in the output.
For nlmixr2, this means string expressions can be converted to mu referenced expressions too. Here is an example:
one.compartment <- function() {
ini({
tka <- log(1.57); label("Ka")
tcl <- log(2.72); label("Cl")
tv <- log(31.5); label("V")
wt.cl <- 0.75; label("WT on CL")
eta.ka ~ 0.6
eta.cl ~ 0.3
eta.v ~ 0.1
add.sd <- 0.7
})
# and a model block with the error specification and model specification
model({
if (WT < 70) {
wt70 <- "Less than 70"
} else {
wt70 <- "Greater than 70"
}
ka <- exp(tka + eta.ka)
cl <- exp(tcl + eta.cl + wt.cl*(wt70 == "Less than 70"))
v <- exp(tv + eta.v)
d/dt(depot) <- -ka * depot
d/dt(center) <- ka * depot - cl / v * center
cp <- center / v
cp ~ add(add.sd)
})
}
one <- one.compartment()
f <- nlmixr2(one, theo_sd, "saem", control=saemControl(print=0))## ℹ loading model to look for μ₃/μ₄ references## → loading into symengine environment...## → pruning branches (`if`/`else`) of saem model...## ✔ done## ℹ done## ℹ μ₄ item: (wt70 == "Less than 70")## → loading into symengine environment...## → pruning branches (`if`/`else`) of saem model...## ✔ done## → finding duplicate expressions in saem model...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → optimizing duplicate expressions in saem model...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## ✔ done## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’## ℹ calculate uninformed etas## ℹ done## Calculating covariance matrix## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → loading into symengine environment...## → pruning branches (`if`/`else`) of saem model...## ✔ done## → finding duplicate expressions in saem predOnly model 0...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → finding duplicate expressions in saem predOnly model 1...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → optimizing duplicate expressions in saem predOnly model 1...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → finding duplicate expressions in saem predOnly model 2...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## → optimizing duplicate expressions in saem predOnly model 2...## [====|====|====|====|====|====|====|====|====|====] 0:00:00## ✔ done## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’## → Calculating residuals/tables## ✔ done## → compress origData in nlmixr2 object, save 6600## → compress phiM in nlmixr2 object, save 61240## → compress parHistData in nlmixr2 object, save 14864## → compress saem0 in nlmixr2 object, save 30480print(f)## ── nlmixr² SAEM OBJF by FOCEi approximation ──
##
## Gaussian/Laplacian Likelihoods: AIC() or $objf etc.
## FOCEi CWRES & Likelihoods: addCwres()
##
## ── Time (sec $time): ──
##
## setup covariance saem table compress other
## elapsed 0.001763 0.008005 5.204 0.058 0.019 1.907232
##
## ── Population Parameters ($parFixed or $parFixedDf): ──
##
## Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%)
## tka Ka 0.465 0.196 42 1.59 (1.09, 2.34) 71.4
## tcl Cl 1.05 0.108 10.3 2.84 (2.3, 3.51) 26.5
## tv V 3.46 0.046 1.33 31.7 (29, 34.7) 13.6
## wt.cl WT on CL -0.0911 0.166 182 -0.0911 (-0.416, 0.234)
## add.sd 0.697 0.697
## Shrink(SD)%
## tka 0.174%
## tcl 4.75%
## tv 10.8%
## wt.cl
## add.sd
##
## Covariance Type ($covMethod): linFim
## No correlations in between subject variability (BSV) matrix
## Full BSV covariance ($omega) or correlation ($omegaR; diagonals=SDs)
## Distribution stats (mean/skewness/kurtosis/p-value) available in $shrink
## Information about run found ($runInfo):
## • μ₄ item: (wt70 == "Less than 70")
## Censoring ($censInformation): No censoring
##
## ── Fit Data (object is a modified tibble): ──
## # A tibble: 132 × 21
## ID TIME DV PRED RES IPRED IRES IWRES eta.ka eta.cl eta.v cp
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 0 0.74 0 0.74 0 0.74 1.06 0.0689 -0.499 -0.0907 0
## 2 1 0.25 2.84 3.28 -0.436 3.81 -0.968 -1.39 0.0689 -0.499 -0.0907 3.81
## 3 1 0.57 6.57 5.85 0.721 6.74 -0.170 -0.244 0.0689 -0.499 -0.0907 6.74
## # ℹ 129 more rows
## # ℹ 9 more variables: depot <dbl>, center <dbl>, wt70 <fct>, ka <dbl>,
## # cl <dbl>, v <dbl>, tad <dbl>, dosenum <dbl>, WT <dbl>As above, you can see this mu-referenced expression by the $runInfo:
f$runInfo## [1] "μ₄ item: (wt70 == \"Less than 70\")"Why the Pokemon icon?
Why the icon for mu2 referencing? Every time I hear mu2 I think of the ultimate genetically modified Pokemon mewtwo (which is a link to where the image comes from).
Since this was expanded from the original article I used it as the base, even though now we are at mu3 and mu4 referencing.
Sicne this article came out in December, the image was expanded with holiday lights.