User defined models
Source:vignettes/posologyr_user_defined_models.Rmd
posologyr_user_defined_models.RmdIntroduction
This describes the structure of prior models usable by
posologyr and illustrates how to define new models from
published population pharmacokinetic (ppk) models.
Structure
A posologyr prior ppk model is a named R list:
- ppk_model
- A rxode2 model implementing the structural population pharmacokinetics 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 to implement 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::rxode() mini-language.
posologyr needs a structural model, defined with either
differential or algebraic equations, and an individual model.
The concentration in the central compartment must be named Cc.
The differential function d/dt(AUC) = Cc; is needed for
the optimisation 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;
})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] + sigma[2]*f #proportional model if sigma[1] == 0
return(g)
}Alternatively, the function can take the simulated concentrations,
and a matrix sigma of the estimates of the parameters for
the residual error model, as in the following example:
error_model <- function(f,sigma){
dv <- cbind(f,1)
g <- diag(dv%*%sigma%*%t(dv)) #sigma is the square matrix of the residual
return(sqrt(g)) #errors
}For multiple endpoint models, error_model must be a
named list with a function for each endpoint
vignette("multiple_endpoints").
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 standard deviation can be computed back withsqrt(log((CV^2)+1)) - Full covariance matrix: the easiest to reuse, but rarely seen in articles
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) | SD | Variance = SD^2 |
|---|---|---|---|
| BSV on CL | 39.8 | 0.383 | 0.147 |
| BSV on Vc | 81.6 | 0.714 | 0.510 |
| BSV on Vp | 57.1 | 0.531 | 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, either in a vector:
sigma = c(additive_a = 3.4, proportional_b = 0.227)in a matrix:
or in a named list, see
vignette("multiple_endpoints"):
sigma = list(
cp=c(additive_a = 0.144, proportional_b = 0.15),
pca=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"))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"