Fit several structural PK models and select the best one based on a corrected Bayesian Information Criterion for mixed effects models.
Models to compare can be defined by rate constants and/or clearances and can include or not nonlinear elimination models.
pkbuild <- function(data=NULL, project=NULL, stat=FALSE, param="clearance", new.dir=".",
MM=FALSE, level=NULL, settings.stat=NULL) Let us use the warfarin PK data in this example:
head(read.csv('data/warfarinPK.csv'))## id time y amount wt age sex
## 1 1 0 . 100 66.7 50 1
## 2 1 24 9.2 . 66.7 50 1
## 3 1 36 8.5 . 66.7 50 1
## 4 1 48 6.4 . 66.7 50 1
## 5 1 72 4.8 . 66.7 50 1
## 6 1 96 3.1 . 66.7 50 1warfarin is administrated orally:
library(Rsmlx)
warfarinPK <- list(
dataFile = "data/warfarinPK.csv",
headerTypes = c("id", "time", "observation", "amount", "contcov", "contcov", "catcov"),
administration = "oral")By default, PK models parameterized with clearance(s) are fitted.
warf.pk1 <- pkbuild(data=warfarinPK, new.dir="warfarinPK")##
## [1] "warfarinPK/pk_kaVCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_Tk0VCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagkaVCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagTk0VCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagTk0V1V2QCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagkaV1V2QCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...Comparison of these models is based on a corrected Bayesian Information Criterion for mixed effects models:
\[{\rm BIC}_{\rm cor} = -2 LL + \log(N)\times d_N + \log(n) \times d_n\] where \(N\) is the number of individuals and \(n\) the total number of observations. Then, \(d_N\) is the number of parameters associated to the variability of the individual parameters and \(d_n\) the number of parameters associated to the variability of the observations.
According to this criterion, the best PK model for this data is oral0_1cpt_TlagTk0VCl.txt.
print(warf.pk1)## $pop.ini
## Tlag Tk0 V Cl
## 0.4975827 2.6767010 8.1342567 0.1317972
##
## $project
## [1] "warfarinPK/pk_TlagTk0VCl.mlxtran"
##
## $model
## [1] "C:/ProgramData/Lixoft/MonolixSuite2019R1/factory/library/pk/oral0_1cpt_TlagTk0VCl.txt"
##
## $data
## $data$dataFile
## [1] "data/warfarinPK.csv"
##
## $data$headerTypes
## [1] "id" "time" "observation" "amount" "contcov"
## [6] "contcov" "catcov"
##
## $data$administration
## [1] "oral"
##
##
## $bic
## model bic.cor
## 1 oral0_1cpt_TlagTk0VCl.txt 713.0319
## 2 oral1_1cpt_TlagkaVCl.txt 716.8131
## 3 oral1_2cpt_TlagkaClV1QV2.txt 718.1846
## 4 oral0_2cpt_TlagTk0ClV1QV2.txt 720.6063
## 5 oral0_1cpt_Tk0VCl.txt 840.9993
## 6 oral1_1cpt_kaVCl.txt 919.0924
##
## $pop.est
## Tlag_pop Tk0_pop V_pop Cl_pop omega_Tlag omega_Tk0
## 0.83981762 1.33397176 8.01427029 0.13202203 0.59192544 0.62607546
## omega_V omega_Cl a b
## 0.22381333 0.29282128 0.33080886 0.07114829Models parameterized with rate constants can be used instead:
warf.pk2 <- pkbuild(data=warfarinPK, new.dir="warfarinPK", param="rate")##
## [1] "warfarinPK/pk_kaVk.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_Tk0Vk.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagkaVk.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagTk0Vk.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagTk0Vk12k21k.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagkaVk12k21k.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagkaVk12k21k13k31k.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...print(warf.pk2$bic)## model bic.cor
## 1 oral1_2cpt_TlagkaVkk12k21.txt 702.7145
## 2 oral0_2cpt_TlagTk0Vkk12k21.txt 706.4141
## 3 oral0_1cpt_TlagTk0Vk.txt 709.6980
## 4 oral1_1cpt_TlagkaVk.txt 715.5151
## 5 oral1_3cpt_TlagkaVkk12k21k13k31.txt 729.4074
## 6 oral0_1cpt_Tk0Vk.txt 836.8413
## 7 oral1_1cpt_kaVk.txt 918.5904print(warf.pk2$pop.est)## Tlag_pop ka_pop V_pop k_pop k12_pop k21_pop
## 0.84167886 1.25322089 7.35773190 0.01526672 0.00846664 0.04796776
## omega_Tlag omega_ka omega_V omega_k omega_k12 omega_k21
## 0.68084643 0.67159549 0.19227345 0.06875926 0.84023440 2.75131335
## a b
## 0.30642230 0.06354331or both:
warf.pk3 <- pkbuild(data=warfarinPK, new.dir="warfarinPK", param="both")##
## [1] "warfarinPK/pk_kaVCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_Tk0VCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagkaVCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagTk0VCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagTk0V1V2QCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagkaV1V2QCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_kaVk.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_Tk0Vk.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagkaVk.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagTk0Vk.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagTk0Vk12k21k.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagkaVk12k21k.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagkaVk12k21k13k31k.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...print(warf.pk3$bic)## model bic.cor
## 1 oral1_2cpt_TlagkaVkk12k21.txt 702.7145
## 2 oral0_2cpt_TlagTk0Vkk12k21.txt 706.4141
## 3 oral0_1cpt_TlagTk0Vk.txt 709.6980
## 4 oral0_1cpt_TlagTk0VCl.txt 713.0319
## 5 oral1_1cpt_TlagkaVk.txt 715.5151
## 6 oral1_1cpt_TlagkaVCl.txt 716.8131
## 7 oral1_2cpt_TlagkaClV1QV2.txt 718.1846
## 8 oral0_2cpt_TlagTk0ClV1QV2.txt 720.6063
## 9 oral1_3cpt_TlagkaVkk12k21k13k31.txt 729.4074
## 10 oral0_1cpt_Tk0Vk.txt 836.8413
## 11 oral0_1cpt_Tk0VCl.txt 840.9993
## 12 oral1_1cpt_kaVk.txt 918.5904
## 13 oral1_1cpt_kaVCl.txt 919.0924If a Monolix project has been already created using the data file, it can be used for providing data information
warf.pk4 <- pkbuild(project="projects/warfarinPK.mlxtran", new.dir="warfarinPK")##
## [1] "warfarinPK/pk_kaVCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_Tk0VCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagkaVCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagTk0VCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagTk0V1V2QCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagkaV1V2QCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...Of course, the results are the same as those obtained using the data file
print(warf.pk4$bic)## model bic.cor
## 1 oral0_1cpt_TlagTk0VCl.txt 703.7249
## 2 oral1_2cpt_TlagkaClV1QV2.txt 705.1475
## 3 oral1_1cpt_TlagkaVCl.txt 706.1825
## 4 oral0_2cpt_TlagTk0ClV1QV2.txt 711.8775
## 5 oral0_1cpt_Tk0VCl.txt 816.8908
## 6 oral1_1cpt_kaVCl.txt 895.0166warf.pk5 <- pkbuild(data=warfarinPK, stat=TRUE, new.dir="warfarinPK")##
## [1] "warfarinPK/pk_kaVCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_Tk0VCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagkaVCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagTk0VCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagTk0V1V2QCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagkaV1V2QCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
##
## direction = "full" will be used for the covariate search
##
## ____________________________________________
## Initialization:
##
## Covariate model:
## sex age wt
## Tlag 0 0 0
## Tk0 0 0 0
## V 0 0 0
## Cl 0 0 0
##
## Correlation model:
## list()
##
## Residual error model:
## y
## "combined2"
## Sampling of the conditional distribution using the initial model ...
##
## Estimated criteria (importanceSampling):
## -2LL s.e. AIC BICc BIC
## 657.94 0.11 677.94 704.86 692.60
## ____________________________________________
## ____________________________________________
## Iteration 1:
##
## Covariate model:
## sex age wt
## Tlag 0 0 0
## Tk0 0 0 0
## V 0 0 1
## Cl 0 0 1
##
## Correlation model:
## list()
##
## Residual error model:
## y
## "combined2"
## Run scenario for model 1 ...
## Estimation of the population parameters...
## Sampling from the conditional distribution...
## Estimation of the log-likelihood...
##
## Estimated criteria (importanceSampling):
## -2LL s.e. AIC BICc BIC
## 627.90 0.11 651.90 681.75 669.49
## ____________________________________________
## Iteration 2:
##
## Covariate model:
## sex age wt
## Tlag 0 0 0
## Tk0 0 0 0
## V 0 0 1
## Cl 0 1 1
##
## Correlation model:
## list()
##
## Residual error model:
## y
## "combined2"
## Run scenario for model 2 ...
## Estimation of the population parameters...
## Sampling from the conditional distribution...
## Estimation of the log-likelihood...
##
## Estimated criteria (importanceSampling):
## -2LL s.e. AIC BICc BIC
## 622.97 0.04 648.97 680.29 668.03
## ____________________________________________
## Iteration 3:
## No difference between two successive iterations
## ____________________________________________
## Final model:
##
## Covariate model:
## sex age wt
## Tlag 0 0 0
## Tk0 0 0 0
## V 0 0 1
## Cl 0 1 1
##
## Correlation model:
## list()
##
## Residual error model:
## y
## "combined2"
##
## Estimated criteria (importanceSampling):
## -2LL s.e. AIC BICc BIC
## 622.97 0.04 648.97 680.29 668.03
## ____________________________________________
## total time: 22.9s
## ____________________________________________print(warf.pk5)## $pop.ini
## Tlag Tk0 V Cl
## 0.4975827 2.6767010 8.1342567 0.1317972
##
## $project
## [1] "warfarinPK/pk_TlagTk0VCl_built.mlxtran"
##
## $model
## [1] "C:/ProgramData/Lixoft/MonolixSuite2019R1/factory/library/pk/oral0_1cpt_TlagTk0VCl.txt"
##
## $data
## $data$dataFile
## [1] "data/warfarinPK.csv"
##
## $data$headerTypes
## [1] "id" "time" "observation" "amount" "contcov"
## [6] "contcov" "catcov"
##
## $data$administration
## [1] "oral"
##
##
## $bic
## model bic.cor
## 1 oral0_1cpt_TlagTk0VCl.txt 713.0319
## 2 oral1_1cpt_TlagkaVCl.txt 716.8131
## 3 oral1_2cpt_TlagkaClV1QV2.txt 718.1846
## 4 oral0_2cpt_TlagTk0ClV1QV2.txt 720.6063
## 5 oral0_1cpt_Tk0VCl.txt 840.9993
## 6 oral1_1cpt_kaVCl.txt 919.0924
##
## $pop.est
## Tlag_pop Tk0_pop V_pop beta_V_cwt Cl_pop
## 0.800605379 1.523056962 7.978402734 0.013493997 0.132219081
## beta_Cl_cage beta_Cl_cwt omega_Tlag omega_Tk0 omega_V
## 0.009411299 0.009340206 0.552894204 0.537567272 0.141679307
## omega_Cl a b
## 0.253304510 0.329102249 0.070179168
##
## $covariate.model
## $covariate.model$Tlag
## sex age wt
## FALSE FALSE FALSE
##
## $covariate.model$Tk0
## sex age wt
## FALSE FALSE FALSE
##
## $covariate.model$V
## sex age wt
## FALSE FALSE TRUE
##
## $covariate.model$Cl
## sex age wt
## FALSE TRUE TRUE
##
##
## $correlation.model
## list()
##
## $error.model
## y
## "combined2"Define the settings of buildmlx in a list
warf.pk6 <- pkbuild(data=warfarinPK, stat=TRUE, new.dir="warfarinPK",
settings.stat =list(covToTransform = c("wt"),
criterion ="AIC"))##
## [1] "warfarinPK/pk_kaVCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_Tk0VCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagkaVCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagTk0VCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagTk0V1V2QCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
## [1] "warfarinPK/pk_TlagkaV1V2QCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...
##
##
## direction = "full" will be used for the covariate search
##
## ____________________________________________
## Initialization:
##
## Covariate model:
## sex age wt
## Tlag 0 0 0
## Tk0 0 0 0
## V 0 0 0
## Cl 0 0 0
##
## Correlation model:
## list()
##
## Residual error model:
## y
## "combined2"
## Sampling of the conditional distribution using the initial model ...
##
## Estimated criteria (importanceSampling):
## -2LL s.e. AIC BICc BIC
## 657.94 0.11 677.94 704.86 692.60
## ____________________________________________
## ____________________________________________
## Iteration 1:
##
## Covariate model:
## sex age wt lwt
## Tlag 0 0 0 0
## Tk0 0 0 0 0
## V 0 0 0 1
## Cl 0 1 0 1
##
## Correlation model:
## list()
##
## Residual error model:
## y
## "combined2"
## Run scenario for model 1 ...
## Estimation of the population parameters...
## Sampling from the conditional distribution...
## Estimation of the log-likelihood...
##
## Estimated criteria (importanceSampling):
## -2LL s.e. AIC BICc BIC
## 621.35 0.08 647.35 678.67 666.40
## ____________________________________________
## Iteration 2:
## No difference between two successive iterations
## ____________________________________________
## Final model:
##
## Covariate model:
## sex age wt lwt
## Tlag 0 0 0 0
## Tk0 0 0 0 0
## V 0 0 0 1
## Cl 0 1 0 1
##
## Correlation model:
## list()
##
## Residual error model:
## y
## "combined2"
##
## Estimated criteria (importanceSampling):
## -2LL s.e. AIC BICc BIC
## 621.35 0.08 647.35 678.67 666.40
## ____________________________________________
## total time: 18.6s
## ____________________________________________