This describes the structure of prior models usable by
posologyr and illustrates how to define new models from
published population pharmacokinetic (ppk) models.
A posologyr prior ppk model is a named R list:
vignette("multiple_endpoints")
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).
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;
})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").
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.
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.
sqrt(exp(SD^2)-1), the standard deviation can be computed
back with sqrt(log((CV^2)+1))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.
The estimates of the parameters for the residual error model, either in a vector:
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.
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().
The names of every covariate defined in ppk_model, in a
character vector.
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"))mod_vancomyin_Goti2018
#> $ppk_model
#> rxode2 2.1.2 model named rx_274fc254fbfa90fbebccd323b4a60983 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"