Classic posologyr models
Source:vignettes/articles/classic_posologyr_models.Rmd
classic_posologyr_models.RmdIntroduction
Originally, posologyr models were R lists;
posologyr still uses this format internally, but it’s often
more useful to use the rxode2 syntax, or even to import a
NONMEM model using the nonmem2rx package.
This article describes the structure of classical
posologyr models. It illustrates how to define them from
published population models.
Structure
A posologyr model is a named R list of the following
items:
- ppk_model
-
A
rxode2model implementing the structural population model with the individual model (i.e. the model of inter-individual variability) and the covariates - error_model
-
A function of the residual error model, alternatively a named list of
functions for multiple endpoints model
vignette("multiple_endpoints") - theta
- A named vector of the population estimates of the fixed effects parameters (called THETAs, following NONMEM terminology)
- omega
- A named square variance-covariance matrix of the population parameters inter-individual variability
- sigma
- The estimates of the parameters of the residual error model
- pi_matrix
- Optional. A named square variance-covariance matrix of the population parameters inter-occasion variability
- covariates
- A character vector of the covariates of the model
Definition of a prior model through an example
The model in this example is a two-compartment ppk model of vancomycin derived from a retrospective study with a cohort of over 1,800 patients (doi:10.1097/FTD.0000000000000490).
ppk_model
A model defined in the rxode2 mini-language.
posologyr needs a structural model, defined with either
differential or algebraic equations, and an individual model. Depending
on model type, naming conventions are more or less strict:
Single endpoint model (e.g. most pharmacokinetic models)
- The concentration in the central compartment must be named Cc.
Multiple endpoints model (eg. PK-PD, parent-metabolite, blood-urine…)
- The names of the endpoints are flexible, but they must be consistent with the names of the error models and their parameters.
The differential function d/dt(AUC) = Cc; is needed for
the optimization function poso_dose_auc().
ppk_model = rxode2::rxode({
centr(0) = 0;
TVCl = THETA_Cl*(CLCREAT/120)^0.8*(0.7^DIAL);
TVVc = THETA_Vc*(WT/70) *(0.5^DIAL);
TVVp = THETA_Vp;
TVQ = THETA_Q;
Cl = TVCl*exp(ETA_Cl);
Vc = TVVc*exp(ETA_Vc);
Vp = TVVp*exp(ETA_Vp);
Q = TVQ;
ke = Cl/Vc;
k12 = Q/Vc;
k21 = Q/Vp;
Cc = centr/Vc;
d/dt(centr) = - ke*centr - k12*centr + k21*periph;
d/dt(periph) = + k12*centr - k21*periph;
d/dt(AUC) = Cc;
})
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’error_model
A function of the residual error model, taking two arguments: the
simulated concentrations, and a vector sigma of the
estimates of the parameters for the residual error model.
error_model <- function(f,sigma){ #additive model if sigma[2] == 0
g <- sigma[1]^2 + (sigma[2]^2)*(f^2) #proportional model if sigma[1] == 0
return(sqrt(g))
}For multiple endpoints models, error_model must be a
named list with a function for each endpoint
vignette("multiple_endpoints").
error_model = list(
first_endpoint = function(f,sigma){
g <- sigma[1]^2 + (sigma[2]^2)*(f^2)
return(sqrt(g))
},
second_endpoint = function(f,sigma){
g <- sigma[1]^2 + (sigma[2]^2)*(f^2)
return(sqrt(g))
}
)To obtain individual estimations for multiple endpoints, consistency in the naming convention must be maintained across the following:
- The dataset (using the column DVID).
- The residual error models (stored in a named list).
- The standard deviation of the residual error models (stored in a named list, see sigma).
As many residual error models as desired can be defined. Each model
defined in the named list error_models must have its counterpart in the
named list sigma, and the names must match those defined in
the DVID column of the dataset.
theta
The estimations of the parameters for the fixed effects of the model
(THETA), in a named vector. The names must match the names used in
ppk_model.
theta = c(THETA_Cl=4.5, THETA_Vc=58.4, THETA_Vp=38.4, THETA_Q=6.5)omega
The variance-covariance matrix of the random effects (ETA) for the
individual model. A symmetric matrix. The names must match the names
used in ppk_model. An easy way to define it is using
lotri::lotri().
The estimates of the variances of the random effects can be given under different parameterizations depending on the authors.
- Standard deviation (SD): the square root of the variance, as returned by Monolix
- Coefficient of variation (CV): calculated as
sqrt(exp(SD^2)-1), the variance can be computed back withlog((CV^2)+1) - Full covariance matrix: the easiest to reuse, but less common in the literature
In the case of the vancomycin
model, the estimates of between subject variability (BSV) are given
as CV%. They must be converted to variances prior to their inclusion in
omega.
| Parameter | CV% (from the article) | Variance = SD^2 |
|---|---|---|
| BSV on CL | 39.8 | 0.147 |
| BSV on Vc | 81.6 | 0.510 |
| BSV on Vp | 57.1 | 0.282 |
omega = lotri::lotri({ETA_Cl + ETA_Vc + ETA_Vp + ETA_Q ~
c(0.147,
0 , 0.510 ,
0 , 0 , 0.282,
0 , 0 , 0 , 0)})The estimates of covariance (off-diagonal) are sometimes given as
coefficients of correlation between ETAs. The covariance between ETA_a
and ETA_b can be computed with the following product:
standard_deviation(ETA_a) * standard_deviation(ETA_b) * correlation(ETA_a and ETA_b).
In this example, all covariances are equal to zero.
sigma
The estimates of the parameters for the residual error model on the standard deviation scale, either in a vector:
sigma = c(additive_a = 3.4, proportional_b = 0.227)in a matrix:
or in a named list (for multiple endpoints):
sigma = list(
first_endpoint=c(additive_a = 0.144, proportional_b = 0.15),
second_endpoint=c(additive_a = 3.91, proportional_b = 0.0)
)depending on the residual error model.
pi_matrix
Optional: only needed for models with inter-occasion variability
(IOV). The variance-covariance matrix of the random effects (KAPPA) for
the IOV. As for the omega matrix, the names must match the
names used in ppk_model. An easy way to define it is using
lotri::lotri().
covariates
The names of every covariate defined in ppk_model, in a
character vector.
covariates = c("CLCREAT","WT","DIAL")Full model
The posologyr model is the list of all these objects.
Note: this model does not include inter-occasion variability, so the
pi_matrix is omitted.
mod_vancomyin_Goti2018 <- list(
ppk_model = rxode2::rxode({
centr(0) = 0;
TVCl = THETA_Cl*(CLCREAT/120)^0.8*(0.7^DIAL);
TVVc = THETA_Vc*(WT/70) *(0.5^DIAL);
TVVp = THETA_Vp;
TVQ = THETA_Q;
Cl = TVCl*exp(ETA_Cl);
Vc = TVVc*exp(ETA_Vc);
Vp = TVVp*exp(ETA_Vp);
Q = TVQ;
ke = Cl/Vc;
k12 = Q/Vc;
k21 = Q/Vp;
Cc = centr/Vc;
d/dt(centr) = - ke*centr - k12*centr + k21*periph;
d/dt(periph) = + k12*centr - k21*periph;
d/dt(AUC) = Cc;
}),
error_model = function(f,sigma){
g <- sigma[1] + sigma[2]*f
return(g)
},
theta = c(THETA_Cl=4.5, THETA_Vc=58.4, THETA_Vp=38.4,THETA_Q=6.5),
omega = lotri::lotri({ETA_Cl + ETA_Vc + ETA_Vp + ETA_Q ~
c(0.147,
0 , 0.510 ,
0 , 0 , 0.282,
0 , 0 , 0 , 0)}),
sigma = c(additive_a = 3.4, proportional_b = 0.227),
covariates = c("CLCREAT","WT","DIAL"))
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’Resulting R object
mod_vancomyin_Goti2018
#> $ppk_model
#> rxode2 NA model named rx_bce7176cfa18100af62e71ae38d3cd48 model (✔ ready).
#> $state: centr, periph, AUC
#> $params: THETA_Cl, CLCREAT, DIAL, THETA_Vc, WT, THETA_Vp, THETA_Q, ETA_Cl, ETA_Vc, ETA_Vp
#> $lhs: TVCl, TVVc, TVVp, TVQ, Cl, Vc, Vp, Q, ke, k12, k21, Cc
#>
#> $error_model
#> function(f,sigma){
#> g <- sigma[1] + sigma[2]*f
#> return(g)
#> }
#>
#> $theta
#> THETA_Cl THETA_Vc THETA_Vp THETA_Q
#> 4.5 58.4 38.4 6.5
#>
#> $omega
#> ETA_Cl ETA_Vc ETA_Vp ETA_Q
#> ETA_Cl 0.147 0.00 0.000 0
#> ETA_Vc 0.000 0.51 0.000 0
#> ETA_Vp 0.000 0.00 0.282 0
#> ETA_Q 0.000 0.00 0.000 0
#>
#> $sigma
#> additive_a proportional_b
#> 3.400 0.227
#>
#> $covariates
#> [1] "CLCREAT" "WT" "DIAL"