library(Rsmlx)
library(ggplot2)
theme_set(theme_bw())

Read and plot the data

Let us first read and plot the phenobarbital data

phenobarbital <- list(
  dataFile = "data/phenobarbital_data.csv",
  headerTypes = c("id", "time", "amount", "contcov", "contcov", "observation"),
  administration = "iv")

phenobarbital.data <- mlxR::readDatamlx(data=phenobarbital)
ggplot(phenobarbital.data$y, aes(x=time, y=y, color=id)) + geom_point() + geom_line() + theme(legend.position = "none") 


Select the structural model

We use pkbuild to build a PK model for this data, using clearance as parameterization:

phenobarbital.clearance <- pkbuild(data=phenobarbital, new.dir="phenobarbital")
## 
## [1] "phenobarbital/pk_VCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
## 
## [1] "phenobarbital/pk_V1V2QCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...

According to (corrected) BIC, the best PK models for this data are the following ones

print(phenobarbital.clearance$bic)
##                       model  bic.cor
## 1     infusion_1cpt_VCl.txt 1039.733
## 2 infusion_2cpt_ClV1QV2.txt 1057.660


Select the statistical model

Once the structural model has been selected (infusion_1cpt_VCl.txt in this example), we can use buildmlx to build the statistical model

phenobarbital.res <- buildmlx(project=phenobarbital.clearance$project, covToTransform = "all")
## 
## 
## direction = "full" will be used for the covariate search
## 
## ____________________________________________
## Initialization:
## 
## Covariate model:
##    APGAR WEIGHT
## V      0      0
## Cl     0      0
## 
## Correlation model:
## list()
## 
## Residual error model:
##        CONC 
## "combined2" 
## Sampling of the conditional distribution using the initial model ... 
## 
## Estimated criteria (importanceSampling):
##    -2LL    s.e.     AIC    BICc     BIC 
## 1009.47    0.07 1021.47 1037.80 1033.94 
## ____________________________________________
## ____________________________________________
## Iteration 1:
## 
## Covariate model:
##    APGAR WEIGHT lAPGAR lWEIGHT
## V      0      1      0       0
## Cl     0      0      0       1
## 
## Correlation model:
## list()
## 
## Residual error model:
##       CONC 
## "constant" 
## 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 
## 872.58   0.06 886.58 904.02 901.12 
## ____________________________________________
## Iteration 2:
## No difference between two successive iterations
## ____________________________________________
## Final model:
## 
## Covariate model:
##    APGAR WEIGHT lAPGAR lWEIGHT
## V      0      1      0       0
## Cl     0      0      0       1
## 
## Correlation model:
## list()
## 
## Residual error model:
##       CONC 
## "constant" 
## 
## Estimated criteria (importanceSampling):
##   -2LL   s.e.    AIC   BICc    BIC 
## 872.58   0.06 886.58 904.02 901.12 
## ____________________________________________
## total time: 15.1s
## ____________________________________________


Using Monolix to fit the data with the selected model

The following Monolix project was created by buildmlx

print(phenobarbital.res$project)
## [1] "phenobarbital/pk_VCl_built.mlxtran"

We can load and run this project with Monolix (either using the MlxConnectors or the Monolix graphical user interface). Here are some individual fits obtained with Monolix:


Confidence intervals and statistical tests

90% confidence intervals for the population parameters can be computed using confintmlx with the estimated Fisher Information Matrix

pheno.fim <- confintmlx(project=phenobarbital.res$project)
## Estimation of the standard errors ...
print(pheno.fim$confint)
##                   estimate       lower       upper
## V_pop           1.38568091 1.327495982 1.446416119
## beta_V_cWEIGHT  0.52206417 0.459171104 0.584957238
## Cl_pop          0.00697066 0.006378209 0.007618142
## beta_Cl_lWEIGHT 1.19789900 0.948569615 1.447228379
## omega_V         0.15855692 0.126752232 0.198342032
## omega_Cl        0.24879614 0.179627300 0.344599718
## a               2.68247516 2.292205779 3.072744546

using the profile likelihood

pheno.prl <- confintmlx(project=phenobarbital.res$project, method="proflike")
## /**********************************************************************/ 
##  LL search on V_pop
## Upper bound search / Iteration 2 / V_pop = 1.692
## Upper bound search / Iteration 3 / V_pop = 1.461
## Upper bound search / Iteration 4 / V_pop = 1.452
## Lower bound search / Iteration 2 / V_pop = 1.134
## Lower bound search / Iteration 3 / V_pop = 1.315
## Lower bound search / Iteration 4 / V_pop = 1.318
## parameter  V_pop 
## Value  1.385 
## CI =  [1.315 , 1.452]  
## diff. =  [-0.071 , 0.066] 
## rel. diff. =  [-5.097 , 4.832] 
## /**********************************************************************/ 
##  LL search on beta_V_cWEIGHT
## Upper bound search / Iteration 2 / beta_V_cWEIGHT = 0.722
## Upper bound search / Iteration 3 / beta_V_cWEIGHT = 0.599
## Upper bound search / Iteration 4 / beta_V_cWEIGHT = 0.608
## Upper bound search / Iteration 5 / beta_V_cWEIGHT = 0.6
## Upper bound search / Iteration 6 / beta_V_cWEIGHT = 0.606
## Upper bound search / Iteration 7 / beta_V_cWEIGHT = 0.602
## Lower bound search / Iteration 2 / beta_V_cWEIGHT = 0.322
## Lower bound search / Iteration 3 / beta_V_cWEIGHT = 0.453
## Lower bound search / Iteration 4 / beta_V_cWEIGHT = 0.464
## Lower bound search / Iteration 5 / beta_V_cWEIGHT = 0.461
## parameter  beta_V_cWEIGHT 
## Value  0.522 
## CI =  [0.461 , 0.606]  
## diff. =  [-0.061 , 0.084] 
## rel. diff. =  [-11.528 , 16.13] 
## /**********************************************************************/ 
##  LL search on Cl_pop
## Upper bound search / Iteration 2 / Cl_pop = 0.008
## Upper bound search / Iteration 3 / Cl_pop = 0.007
## Upper bound search / Iteration 4 / Cl_pop = 0.007
## Lower bound search / Iteration 2 / Cl_pop = 0.005
## Lower bound search / Iteration 3 / Cl_pop = 0.006
## Lower bound search / Iteration 4 / Cl_pop = 0.006
## Lower bound search / Iteration 5 / Cl_pop = 0.006
## Lower bound search / Iteration 6 / Cl_pop = 0.006
## parameter  Cl_pop 
## Value  0.006 
## CI =  [0.006 , 0.007]  
## diff. =  [-0.001 , 0] 
## rel. diff. =  [-8.389 , 9.648] 
## /**********************************************************************/ 
##  LL search on beta_Cl_lWEIGHT
## Upper bound search / Iteration 2 / beta_Cl_lWEIGHT = 1.397
## Upper bound search / Iteration 3 / beta_Cl_lWEIGHT = 1.364
## Upper bound search / Iteration 4 / beta_Cl_lWEIGHT = 1.355
## Lower bound search / Iteration 2 / beta_Cl_lWEIGHT = 0.997
## Lower bound search / Iteration 3 / beta_Cl_lWEIGHT = -0.447
## Lower bound search / Iteration 4 / beta_Cl_lWEIGHT = 0.926
## Lower bound search / Iteration 5 / beta_Cl_lWEIGHT = 0.811
## Lower bound search / Iteration 6 / beta_Cl_lWEIGHT = 0.887
## Lower bound search / Iteration 7 / beta_Cl_lWEIGHT = 0.872
## Lower bound search / Iteration 8 / beta_Cl_lWEIGHT = 0.828
## Lower bound search / Iteration 9 / beta_Cl_lWEIGHT = 0.855
## parameter  beta_Cl_lWEIGHT 
## Value  1.197 
## CI =  [0.855 , 1.364]  
## diff. =  [-0.343 , 0.166] 
## rel. diff. =  [-28.582 , 13.93] 
## /**********************************************************************/ 
##  LL search on omega_Cl
## Upper bound search / Iteration 2 / omega_Cl = 0.303
## Upper bound search / Iteration 3 / omega_Cl = 0.311
## Lower bound search / Iteration 2 / omega_Cl = 0.203
## Lower bound search / Iteration 3 / omega_Cl = 0.048
## Lower bound search / Iteration 4 / omega_Cl = 0.168
## Lower bound search / Iteration 5 / omega_Cl = 0.146
## Lower bound search / Iteration 6 / omega_Cl = 0.062
## Lower bound search / Iteration 7 / omega_Cl = 0.138
## Lower bound search / Iteration 8 / omega_Cl = 0.141
## parameter  omega_Cl 
## Value  0.248 
## CI =  [0.141 , 0.311]  
## diff. =  [-0.108 , 0.062] 
## rel. diff. =  [-43.045 , 25.07] 
## /**********************************************************************/ 
##  LL search on omega_V
## Upper bound search / Iteration 2 / omega_V = 0.193
## Upper bound search / Iteration 3 / omega_V = 0.236
## Upper bound search / Iteration 4 / omega_V = 0.211
## Upper bound search / Iteration 5 / omega_V = 0.228
## Upper bound search / Iteration 6 / omega_V = 0.214
## Upper bound search / Iteration 7 / omega_V = 0.223
## Upper bound search / Iteration 8 / omega_V = 0.215
## Lower bound search / Iteration 2 / omega_V = 0.129
## parameter  omega_V 
## Value  0.158 
## CI =  [0.129 , 0.215]  
## diff. =  [-0.029 , 0.056] 
## rel. diff. =  [-18.127 , 35.918] 
## /**********************************************************************/ 
##  LL search on a
## Upper bound search / Iteration 2 / a = 3.276
## Upper bound search / Iteration 3 / a = 3.099
## Upper bound search / Iteration 4 / a = 3.126
## Lower bound search / Iteration 2 / a = 2.196
## Lower bound search / Iteration 3 / a = 2.358
## Lower bound search / Iteration 4 / a = 2.363

## parameter  a 
## Value  2.682 
## CI =  [2.363 , 3.126]  
## diff. =  [-0.319 , 0.444] 
## rel. diff. =  [-11.877 , 16.56]
## /**********************************************************************/ 
## parameter  V_pop 
## Value  1.385 
## CI =  [1.315 , 1.452]  
## diff. =  [-0.071 , 0.066] 
## rel. diff. =  [-5.097 , 4.832] 
## /**********************************************************************/ 
## parameter  beta_V_cWEIGHT 
## Value  0.522 
## CI =  [0.461 , 0.606]  
## diff. =  [-0.061 , 0.084] 
## rel. diff. =  [-11.528 , 16.13] 
## /**********************************************************************/ 
## parameter  Cl_pop 
## Value  0.006 
## CI =  [0.006 , 0.007]  
## diff. =  [-0.001 , 0] 
## rel. diff. =  [-8.389 , 9.648] 
## /**********************************************************************/ 
## parameter  beta_Cl_lWEIGHT 
## Value  1.197 
## CI =  [0.855 , 1.364]  
## diff. =  [-0.343 , 0.166] 
## rel. diff. =  [-28.582 , 13.93] 
## /**********************************************************************/ 
## parameter  omega_Cl 
## Value  0.248 
## CI =  [0.141 , 0.311]  
## diff. =  [-0.108 , 0.062] 
## rel. diff. =  [-43.045 , 25.07] 
## /**********************************************************************/ 
## parameter  omega_V 
## Value  0.158 
## CI =  [0.129 , 0.215]  
## diff. =  [-0.029 , 0.056] 
## rel. diff. =  [-18.127 , 35.918] 
## /**********************************************************************/ 
## parameter  a 
## Value  2.682 
## CI =  [2.363 , 3.126]  
## diff. =  [-0.319 , 0.444] 
## rel. diff. =  [-11.877 , 16.56]
print(pheno.prl$confint)
##                   estimate       lower      upper
## V_pop           1.38568000 1.315064085 1.45264201
## beta_V_cWEIGHT  0.52206420 0.461880898 0.60627685
## Cl_pop          0.00697066 0.006385915 0.00764319
## beta_Cl_lWEIGHT 1.19790000 0.855519366 1.36477072
## omega_Cl        0.24879610 0.141703342 0.31116996
## omega_V         0.15855690 0.129815410 0.21550782
## a               2.68248000 2.363883711 3.12671918

or non parametric bootstrapping

pheno.boot <- confintmlx(project=phenobarbital.res$project, method="bootstrap")
## Generating data sets with initial data set resampling...
## Generating projects with bootstrap data sets...
## Project 1/100 => Estimating the population parameters 
## Project 2/100 => Estimating the population parameters 
## Project 3/100 => Estimating the population parameters 
## Project 4/100 => Estimating the population parameters 
## Project 5/100 => Estimating the population parameters 
## Project 6/100 => Estimating the population parameters 
## Project 7/100 => Estimating the population parameters 
## Project 8/100 => Estimating the population parameters 
## Project 9/100 => Estimating the population parameters 
## Project 10/100 => Estimating the population parameters 
## Project 11/100 => Estimating the population parameters 
## Project 12/100 => Estimating the population parameters 
## Project 13/100 => Estimating the population parameters 
## Project 14/100 => Estimating the population parameters 
## Project 15/100 => Estimating the population parameters 
## Project 16/100 => Estimating the population parameters 
## Project 17/100 => Estimating the population parameters 
## Project 18/100 => Estimating the population parameters 
## Project 19/100 => Estimating the population parameters 
## Project 20/100 => Estimating the population parameters 
## Project 21/100 => Estimating the population parameters 
## Project 22/100 => Estimating the population parameters 
## Project 23/100 => Estimating the population parameters 
## Project 24/100 => Estimating the population parameters 
## Project 25/100 => Estimating the population parameters 
## Project 26/100 => Estimating the population parameters 
## Project 27/100 => Estimating the population parameters 
## Project 28/100 => Estimating the population parameters 
## Project 29/100 => Estimating the population parameters 
## Project 30/100 => Estimating the population parameters 
## Project 31/100 => Estimating the population parameters 
## Project 32/100 => Estimating the population parameters 
## Project 33/100 => Estimating the population parameters 
## Project 34/100 => Estimating the population parameters 
## Project 35/100 => Estimating the population parameters 
## Project 36/100 => Estimating the population parameters 
## Project 37/100 => Estimating the population parameters 
## Project 38/100 => Estimating the population parameters 
## Project 39/100 => Estimating the population parameters 
## Project 40/100 => Estimating the population parameters 
## Project 41/100 => Estimating the population parameters 
## Project 42/100 => Estimating the population parameters 
## Project 43/100 => Estimating the population parameters 
## Project 44/100 => Estimating the population parameters 
## Project 45/100 => Estimating the population parameters 
## Project 46/100 => Estimating the population parameters 
## Project 47/100 => Estimating the population parameters 
## Project 48/100 => Estimating the population parameters 
## Project 49/100 => Estimating the population parameters 
## Project 50/100 => Estimating the population parameters 
## Project 51/100 => Estimating the population parameters 
## Project 52/100 => Estimating the population parameters 
## Project 53/100 => Estimating the population parameters 
## Project 54/100 => Estimating the population parameters 
## Project 55/100 => Estimating the population parameters 
## Project 56/100 => Estimating the population parameters 
## Project 57/100 => Estimating the population parameters 
## Project 58/100 => Estimating the population parameters 
## Project 59/100 => Estimating the population parameters 
## Project 60/100 => Estimating the population parameters 
## Project 61/100 => Estimating the population parameters 
## Project 62/100 => Estimating the population parameters 
## Project 63/100 => Estimating the population parameters 
## Project 64/100 => Estimating the population parameters 
## Project 65/100 => Estimating the population parameters 
## Project 66/100 => Estimating the population parameters 
## Project 67/100 => Estimating the population parameters 
## Project 68/100 => Estimating the population parameters 
## Project 69/100 => Estimating the population parameters 
## Project 70/100 => Estimating the population parameters 
## Project 71/100 => Estimating the population parameters 
## Project 72/100 => Estimating the population parameters 
## Project 73/100 => Estimating the population parameters 
## Project 74/100 => Estimating the population parameters 
## Project 75/100 => Estimating the population parameters 
## Project 76/100 => Estimating the population parameters 
## Project 77/100 => Estimating the population parameters 
## Project 78/100 => Estimating the population parameters 
## Project 79/100 => Estimating the population parameters 
## Project 80/100 => Estimating the population parameters 
## Project 81/100 => Estimating the population parameters 
## Project 82/100 => Estimating the population parameters 
## Project 83/100 => Estimating the population parameters 
## Project 84/100 => Estimating the population parameters 
## Project 85/100 => Estimating the population parameters 
## Project 86/100 => Estimating the population parameters 
## Project 87/100 => Estimating the population parameters 
## Project 88/100 => Estimating the population parameters 
## Project 89/100 => Estimating the population parameters 
## Project 90/100 => Estimating the population parameters 
## Project 91/100 => Estimating the population parameters 
## Project 92/100 => Estimating the population parameters 
## Project 93/100 => Estimating the population parameters 
## Project 94/100 => Estimating the population parameters 
## Project 95/100 => Estimating the population parameters 
## Project 96/100 => Estimating the population parameters 
## Project 97/100 => Estimating the population parameters 
## Project 98/100 => Estimating the population parameters 
## Project 99/100 => Estimating the population parameters 
## Project 100/100 => Estimating the population parameters
print(pheno.boot$confint)
##                   estimate       lower       upper
## V_pop           1.38568091 1.342239240 1.438484557
## beta_V_cWEIGHT  0.52206417 0.481044952 0.611446067
## Cl_pop          0.00697066 0.006421249 0.007448236
## beta_Cl_lWEIGHT 1.19789900 0.780601331 1.411766575
## omega_V         0.15855692 0.124206298 0.205383578
## omega_Cl        0.24879614 0.088396813 0.335811728
## a               2.68247516 2.311560665 3.126748665