library(Rsmlx)
initializeMlxConnectors(software = "monolix")
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 <- 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"
## [INFO] Results have been successfully loaded
## 
## [1] "phenobarbital/pk_V1V2QCl.mlxtran"
## [INFO] Results have been successfully loaded

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 1038.771
## 2 infusion_2cpt_ClV1QV2.txt 1058.571


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")
## [INFO] Results have been successfully loaded
## 
## 
## direction = "full" will be used for the covariate search
## 
## ____________________________________________
## Initialization:
## 
## Covariate model:
## [1] 0 0
## 
## Correlation model:
## list()
## 
## Residual error model:
##        CONC 
## "combined2" 
## 
## Estimated criteria (importanceSampling):
##    -2LL    s.e.     AIC     BIC 
## 1008.51    0.06 1020.51 1032.98 
## ____________________________________________
## ____________________________________________
## Iteration 1:
## 
## Covariate model:
##    WEIGHT lWEIGHT
## V       1       0
## Cl      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    BIC 
## 872.08   0.06 886.08 900.63 
## ____________________________________________
## Iteration 2:
## No difference between two successive iterations
## ____________________________________________
## Final model:
## 
## Covariate model:
##    WEIGHT lWEIGHT
## V       1       0
## Cl      0       1
## 
## Correlation model:
## list()
## 
## Residual error model:
##       CONC 
## "constant" 
## 
## Estimated criteria (importanceSampling):
##   -2LL   s.e.    AIC    BIC 
## 872.08   0.06 886.08 900.63 
## ____________________________________________
## total time: 6.4s
## ____________________________________________


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 population parameters... 
## [INFO] Results have been successfully loaded
## Estimation of the standard errors ... 
## [INFO] The linearization methods is more relevant around the conditional mode, we encourage you to make this individual mode task in your scenario
print(pheno.fim$confint)
##                    estimate       lower       upper
## V_pop           1.377913577 1.317480341 1.441118904
## beta_V_cWEIGHT  0.534840156 0.469385808 0.600294503
## Cl_pop          0.007096061 0.006567173 0.007667543
## beta_Cl_lWEIGHT 1.098755289 0.879262976 1.318247602
## omega_Cl        0.194226056 0.130896995 0.288194248
## omega_V         0.170148377 0.137605056 0.210388127
## a               2.712560877 2.323560262 3.101561493

using the profile likelihood

pheno.prl <- confintmlx(project=phenobarbital.res$project, method="proflike")
## [INFO] Results have been successfully loaded
## [INFO] Results have been successfully loaded
## /**********************************************************************/ 
##  LL search on V_pop
## Upper bound search / Iteration 2 / V_pop = 1.682
## Upper bound search / Iteration 3 / V_pop = 1.455
## Upper bound search / Iteration 4 / V_pop = 1.448
## Lower bound search / Iteration 2 / V_pop = 1.128
## Lower bound search / Iteration 3 / V_pop = 1.309
## Lower bound search / Iteration 4 / V_pop = 1.331
## Lower bound search / Iteration 5 / V_pop = 1.329
## parameter  V_pop 
## Value  1.377 
## CI =  [1.329 , 1.448]  
## diff. =  [-0.048 , 0.07] 
## rel. diff. =  [-3.483 , 5.095] 
## /**********************************************************************/ 
##  LL search on beta_V_cWEIGHT
## Upper bound search / Iteration 2 / beta_V_cWEIGHT = 0.734
## Upper bound search / Iteration 3 / beta_V_cWEIGHT = 0.605
## Upper bound search / Iteration 4 / beta_V_cWEIGHT = 0.595
## Upper bound search / Iteration 5 / beta_V_cWEIGHT = 0.589
## Upper bound search / Iteration 6 / beta_V_cWEIGHT = 0.594
## Upper bound search / Iteration 7 / beta_V_cWEIGHT = 0.589
## Lower bound search / Iteration 2 / beta_V_cWEIGHT = 0.334
## Lower bound search / Iteration 3 / beta_V_cWEIGHT = 0.464
## Lower bound search / Iteration 4 / beta_V_cWEIGHT = 0.475
## Lower bound search / Iteration 5 / beta_V_cWEIGHT = 0.477
## parameter  beta_V_cWEIGHT 
## Value  0.534 
## CI =  [0.475 , 0.589]  
## diff. =  [-0.06 , 0.055] 
## rel. diff. =  [-11.054 , 10.306] 
## /**********************************************************************/ 
##  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
## Upper bound search / Iteration 5 / 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
## parameter  Cl_pop 
## Value  0.007 
## CI =  [0.006 , 0.007]  
## diff. =  [-0.001 , 0] 
## rel. diff. =  [-8.495 , 7.121] 
## /**********************************************************************/ 
##  LL search on beta_Cl_lWEIGHT
## Upper bound search / Iteration 2 / beta_Cl_lWEIGHT = 1.298
## Upper bound search / Iteration 3 / beta_Cl_lWEIGHT = 1.319
## Lower bound search / Iteration 2 / beta_Cl_lWEIGHT = 0.898
## Lower bound search / Iteration 3 / beta_Cl_lWEIGHT = 0.903
## parameter  beta_Cl_lWEIGHT 
## Value  1.098 
## CI =  [0.903 , 1.319]  
## diff. =  [-0.195 , 0.221] 
## rel. diff. =  [-17.744 , 20.12] 
## /**********************************************************************/ 
##  LL search on omega_Cl
## Upper bound search / Iteration 2 / omega_Cl = 0.237
## Upper bound search / Iteration 3 / omega_Cl = 0.282
## Upper bound search / Iteration 4 / omega_Cl = 0.344
## Upper bound search / Iteration 5 / omega_Cl = 0.291
## Lower bound search / Iteration 2 / omega_Cl = 0.159
## Lower bound search / Iteration 3 / omega_Cl = 0.145
## Lower bound search / Iteration 4 / omega_Cl = 0.119
## Lower bound search / Iteration 5 / omega_Cl = 0.133
## parameter  omega_Cl 
## Value  0.194 
## CI =  [0.133 , 0.291]  
## diff. =  [-0.061 , 0.097] 
## rel. diff. =  [-31.391 , 50.334] 
## /**********************************************************************/ 
##  LL search on omega_V
## Upper bound search / Iteration 2 / omega_V = 0.207
## Upper bound search / Iteration 3 / omega_V = 0.206
## Lower bound search / Iteration 2 / omega_V = 0.139
## Lower bound search / Iteration 3 / omega_V = 0.138
## parameter  omega_V 
## Value  0.17 
## CI =  ]-Inf , 0.207]  
## diff. =  ]-Inf , 0.037] 
## rel. diff. =  ]-Inf , 22.14] 
## /**********************************************************************/ 
##  LL search on a
## Upper bound search / Iteration 2 / a = 3.313
## Upper bound search / Iteration 3 / a = 3.139
## Upper bound search / Iteration 4 / a = 3.134
## Lower bound search / Iteration 2 / a = 2.22
## Lower bound search / Iteration 3 / a = 2.396
## Lower bound search / Iteration 4 / a = 2.403

## parameter  a 
## Value  2.712 
## CI =  [2.403 , 3.139]  
## diff. =  [-0.31 , 0.426] 
## rel. diff. =  [-11.396 , 15.726]
## /**********************************************************************/ 
## parameter  V_pop 
## Value  1.377 
## CI =  [1.329 , 1.448]  
## diff. =  [-0.048 , 0.07] 
## rel. diff. =  [-3.483 , 5.095] 
## /**********************************************************************/ 
## parameter  beta_V_cWEIGHT 
## Value  0.534 
## CI =  [0.475 , 0.589]  
## diff. =  [-0.06 , 0.055] 
## rel. diff. =  [-11.054 , 10.306] 
## /**********************************************************************/ 
## parameter  Cl_pop 
## Value  0.007 
## CI =  [0.006 , 0.007]  
## diff. =  [-0.001 , 0] 
## rel. diff. =  [-8.495 , 7.121] 
## /**********************************************************************/ 
## parameter  beta_Cl_lWEIGHT 
## Value  1.098 
## CI =  [0.903 , 1.319]  
## diff. =  [-0.195 , 0.221] 
## rel. diff. =  [-17.744 , 20.12] 
## /**********************************************************************/ 
## parameter  omega_Cl 
## Value  0.194 
## CI =  [0.133 , 0.291]  
## diff. =  [-0.061 , 0.097] 
## rel. diff. =  [-31.391 , 50.334] 
## /**********************************************************************/ 
## parameter  omega_V 
## Value  0.17 
## CI =  ]-Inf , 0.207]  
## diff. =  ]-Inf , 0.037] 
## rel. diff. =  ]-Inf , 22.14] 
## /**********************************************************************/ 
## parameter  a 
## Value  2.712 
## CI =  [2.403 , 3.139]  
## diff. =  [-0.31 , 0.426] 
## rel. diff. =  [-11.396 , 15.726]
print(pheno.prl$confint)
##                    estimate       lower       upper
## V_pop           1.377914000 1.329921871 1.448120255
## beta_V_cWEIGHT  0.534840200 0.475719293 0.589964210
## Cl_pop          0.007096061 0.006493293 0.007601384
## beta_Cl_lWEIGHT 1.098755000 0.903792842 1.319829871
## omega_Cl        0.194226100 0.133258341 0.291988004
## omega_V         0.170148400        -Inf 0.207819725
## a               2.712561000 2.403445647 3.139163045

or non parametric bootstrapping

pheno.boot <- confintmlx(project=phenobarbital.res$project, method="bootstrap")
## [INFO] Results have been successfully loaded
## [INFO] Results have been successfully loaded
## [INFO] Results have been successfully loaded
## Generating data sets...
## Generating projects with bootstrap data sets...
## Project 1/100 => Running SAEM 
## Project 2/100 => Running SAEM 
## Project 3/100 => Running SAEM 
## Project 4/100 => Running SAEM 
## Project 5/100 => Running SAEM 
## Project 6/100 => Running SAEM 
## Project 7/100 => Running SAEM 
## Project 8/100 => Running SAEM 
## Project 9/100 => Running SAEM 
## Project 10/100 => Running SAEM 
## Project 11/100 => Running SAEM 
## Project 12/100 => Running SAEM 
## Project 13/100 => Running SAEM 
## Project 14/100 => Running SAEM 
## Project 15/100 => Running SAEM 
## Project 16/100 => Running SAEM 
## Project 17/100 => Running SAEM 
## Project 18/100 => Running SAEM 
## Project 19/100 => Running SAEM 
## Project 20/100 => Running SAEM 
## Project 21/100 => Running SAEM 
## Project 22/100 => Running SAEM 
## Project 23/100 => Running SAEM 
## Project 24/100 => Running SAEM 
## Project 25/100 => Running SAEM 
## Project 26/100 => Running SAEM 
## Project 27/100 => Running SAEM 
## Project 28/100 => Running SAEM 
## Project 29/100 => Running SAEM 
## Project 30/100 => Running SAEM 
## Project 31/100 => Running SAEM 
## Project 32/100 => Running SAEM 
## Project 33/100 => Running SAEM 
## Project 34/100 => Running SAEM 
## Project 35/100 => Running SAEM 
## Project 36/100 => Running SAEM 
## Project 37/100 => Running SAEM 
## Project 38/100 => Running SAEM 
## Project 39/100 => Running SAEM 
## Project 40/100 => Running SAEM 
## Project 41/100 => Running SAEM 
## Project 42/100 => Running SAEM 
## Project 43/100 => Running SAEM 
## Project 44/100 => Running SAEM 
## Project 45/100 => Running SAEM 
## Project 46/100 => Running SAEM 
## Project 47/100 => Running SAEM 
## Project 48/100 => Running SAEM 
## Project 49/100 => Running SAEM 
## Project 50/100 => Running SAEM 
## Project 51/100 => Running SAEM 
## Project 52/100 => Running SAEM 
## Project 53/100 => Running SAEM 
## Project 54/100 => Running SAEM 
## Project 55/100 => Running SAEM 
## Project 56/100 => Running SAEM 
## Project 57/100 => Running SAEM 
## Project 58/100 => Running SAEM 
## Project 59/100 => Running SAEM 
## Project 60/100 => Running SAEM 
## Project 61/100 => Running SAEM 
## Project 62/100 => Running SAEM 
## Project 63/100 => Running SAEM 
## Project 64/100 => Running SAEM 
## Project 65/100 => Running SAEM 
## Project 66/100 => Running SAEM 
## Project 67/100 => Running SAEM 
## Project 68/100 => Running SAEM 
## Project 69/100 => Running SAEM 
## Project 70/100 => Running SAEM 
## Project 71/100 => Running SAEM 
## Project 72/100 => Running SAEM 
## Project 73/100 => Running SAEM 
## Project 74/100 => Running SAEM 
## Project 75/100 => Running SAEM 
## Project 76/100 => Running SAEM 
## Project 77/100 => Running SAEM 
## Project 78/100 => Running SAEM 
## Project 79/100 => Running SAEM 
## Project 80/100 => Running SAEM 
## Project 81/100 => Running SAEM 
## Project 82/100 => Running SAEM 
## Project 83/100 => Running SAEM 
## Project 84/100 => Running SAEM 
## Project 85/100 => Running SAEM 
## Project 86/100 => Running SAEM 
## Project 87/100 => Running SAEM 
## Project 88/100 => Running SAEM 
## Project 89/100 => Running SAEM 
## Project 90/100 => Running SAEM 
## Project 91/100 => Running SAEM 
## Project 92/100 => Running SAEM 
## Project 93/100 => Running SAEM 
## Project 94/100 => Running SAEM 
## Project 95/100 => Running SAEM 
## Project 96/100 => Running SAEM 
## Project 97/100 => Running SAEM 
## Project 98/100 => Running SAEM 
## Project 99/100 => Running SAEM 
## Project 100/100 => Running SAEM 
## [INFO] Results have been successfully loaded
print(pheno.boot$confint)
##                    estimate       lower       upper
## V_pop           1.377913577 1.344912113 1.449963187
## beta_V_cWEIGHT  0.534840156 0.473530703 0.614754523
## Cl_pop          0.007096061 0.006383202 0.007577012
## beta_Cl_lWEIGHT 1.098755289 0.688229307 1.347975875
## omega_Cl        0.194226056 0.090782311 0.315073483
## omega_V         0.170148377 0.123893502 0.200468772
## a               2.712560877 2.331803524 3.139335949