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

Read and plot the data

Let us first read and plot the theophylline data

theophylline <- list(
  dataFile = "data/theophylline.csv",
  headerTypes = c("id", "amount", "time", "observation", "contcov", "catcov"),
  administration = "oral")

theo.data <- mlxR::readDatamlx(data=theophylline)
ggplot(theo.data$y, aes(x=time, y=y, color=id)) + geom_point() + geom_line()


Select the structural model

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

theo.clearance <- pkbuild(data=theophylline, new.dir="theophylline")
## 
## [1] "theophylline/pk_kaVCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
## 
## [1] "theophylline/pk_Tk0VCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
## 
## [1] "theophylline/pk_TlagkaVCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
## 
## [1] "theophylline/pk_TlagTk0VCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
## 
## [1] "theophylline/pk_TlagkaV1V2QCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
## 
## [1] "theophylline/pk_TlagTk0V1V2QCl.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(theo.clearance$bic)
##                           model  bic.cor
## 1      oral1_1cpt_TlagkaVCl.txt 326.1454
## 2     oral0_1cpt_TlagTk0VCl.txt 334.3952
## 3         oral0_1cpt_Tk0VCl.txt 339.8382
## 4  oral1_2cpt_TlagkaClV1QV2.txt 344.6554
## 5 oral0_2cpt_TlagTk0ClV1QV2.txt 368.9396
## 6          oral1_1cpt_kaVCl.txt 383.8522


Select the statistical model

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

res.clearance <- buildmlx(project=theo.clearance$project, covToTransform = "all")
## 
## 
## direction = "full" will be used for the covariate search
## 
## ____________________________________________
## Initialization:
## 
## Covariate model:
##      SEX WEIGHT
## Tlag   0      0
## ka     0      0
## 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 
## 278.27   0.11 298.27 316.94 303.12 
## ____________________________________________
## ____________________________________________
## Iteration 1:
## 
## Covariate model:
##      SEX WEIGHT lWEIGHT
## Tlag   0      0       0
## ka     0      0       0
## V      0      0       1
## Cl     0      0       0
## 
## Correlation model:
## [[1]]
## [1] "V"  "Cl"
## 
## 
## Residual error model:
##           CONC 
## "proportional" 
## 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 
## 260.60   0.10 282.60 299.45 287.93 
## ____________________________________________
## Iteration 2:
## 
## Covariate model:
##      SEX WEIGHT lWEIGHT
## Tlag   0      0       0
## ka     0      1       0
## V      0      0       1
## Cl     0      0       0
## 
## Correlation model:
## [[1]]
## [1] "V"  "Cl"
## 
## 
## Residual error model:
##           CONC 
## "proportional" 
## 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 
## 257.06   0.09 281.06 298.39 286.88 
## ____________________________________________
## Iteration 3:
## No difference between two successive iterations
## ____________________________________________
## Final model:
## 
## Covariate model:
##      SEX WEIGHT lWEIGHT
## Tlag   0      0       0
## ka     0      1       0
## V      0      0       1
## Cl     0      0       0
## 
## Correlation model:
## [[1]]
## [1] "V"  "Cl"
## 
## 
## Residual error model:
##           CONC 
## "proportional" 
## 
## Estimated criteria (importanceSampling):
##   -2LL   s.e.    AIC   BICc    BIC 
## 257.06   0.09 281.06 298.39 286.88 
## ____________________________________________
## total time: 17.1s
## ____________________________________________

According BIC, the best statistical model assumes


Using rates instead of clearances

Let us select now the best model defined by rate constants instead of clearances

theo.rate <- pkbuild(data=theophylline, param="rate", new.dir="theophylline")
## 
## [1] "theophylline/pk_kaVk.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
## 
## [1] "theophylline/pk_Tk0Vk.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
## 
## [1] "theophylline/pk_TlagkaVk.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
## 
## [1] "theophylline/pk_TlagTk0Vk.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
## 
## [1] "theophylline/pk_TlagkaVk12k21k.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
## 
## [1] "theophylline/pk_TlagTk0Vk12k21k.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood...

The model is still a one compartment model with a \(Tlag\):

print(theo.rate$bic)
##                            model  bic.cor
## 1        oral1_1cpt_TlagkaVk.txt 317.0710
## 2       oral0_1cpt_TlagTk0Vk.txt 324.4401
## 3           oral0_1cpt_Tk0Vk.txt 330.7445
## 4  oral1_2cpt_TlagkaVkk12k21.txt 335.9315
## 5 oral0_2cpt_TlagTk0Vkk12k21.txt 359.5002
## 6            oral1_1cpt_kaVk.txt 376.0506

Nevertheless the BIC seems better with a model defined with rates (bic.cor=317.69) rather than clearances (bic.cor=326.22). In fact, this difference is due to different statistical models (no correlation between \(V\) and \(Cl\) implicitely assumes a correlation between \(V\) and \(k\) and vice versa).

Building now the statistical model leads to the same result

res.rate <- buildmlx(project=theo.rate$project, covToTransform = "all")
## 
## 
## direction = "full" will be used for the covariate search
## 
## ____________________________________________
## Initialization:
## 
## Covariate model:
##      SEX WEIGHT
## Tlag   0      0
## ka     0      0
## V      0      0
## k      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 
## 269.20   0.09 289.20 307.86 294.05 
## ____________________________________________
## ____________________________________________
## Iteration 1:
## 
## Covariate model:
##      SEX WEIGHT lWEIGHT
## Tlag   0      0       0
## ka     0      0       0
## V      0      0       1
## k      0      0       0
## 
## Correlation model:
## [[1]]
## [1] "V" "k"
## 
## 
## Residual error model:
##           CONC 
## "proportional" 
## 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 
## 260.28   0.11 282.28 299.13 287.61 
## ____________________________________________
## Iteration 2:
## No difference between two successive iterations
## ____________________________________________
## Final model:
## 
## Covariate model:
##      SEX WEIGHT lWEIGHT
## Tlag   0      0       0
## ka     0      0       0
## V      0      0       1
## k      0      0       0
## 
## Correlation model:
## [[1]]
## [1] "V" "k"
## 
## 
## Residual error model:
##           CONC 
## "proportional" 
## 
## Estimated criteria (importanceSampling):
##   -2LL   s.e.    AIC   BICc    BIC 
## 260.28   0.11 282.28 299.13 287.61 
## ____________________________________________
## total time: 6.4s
## ____________________________________________

Confidence intervals and statistical tests

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

theo.fim <- confintmlx(project=res.clearance$project, level=0.95)
## Estimation of the standard errors ...
print(theo.fim$confint)
##                    estimate         lower       upper
## Tlag_pop         0.11472295  0.0748118084  0.17592618
## ka_pop           2.44326619  1.6227891299  3.67857386
## beta_ka_cWEIGHT  0.04543698 -0.0003126989  0.09118666
## V_pop           33.29369058 30.8939617213 35.87982151
## beta_V_lWEIGHT   0.63320554  0.2861542139  0.98025686
## Cl_pop           2.77644551  2.3969582828  3.21601328
## omega_Tlag       0.60300214  0.3520120073  1.03295222
## omega_ka         0.67968643  0.4394940410  1.05114881
## omega_V          0.12321690  0.0779159941  0.19485606
## omega_Cl         0.25720177  0.1710223697  0.38680760
## corr_V_Cl        0.85125686  0.4825527417  0.96367922
## b                0.08399034  0.0706147657  0.09736592

0 belongs to the 95%CI for \(\beta_{ka, {\em weight}}\) which means that the statistical significance of the effect of WEIGHT on \(\ka\) is rather limited. This is confirmed by the statisical test performed using testmlx

theo.test <- testmlx(project=res.clearance$project)
print(theo.test$covariate$p.value.parameters)
##    parameter covariate p.value in.model
## 6     log.ka   cWEIGHT 0.09485     TRUE
## 11     log.V   lWEIGHT 0.01415     TRUE
## 3   log.Tlag   lWEIGHT 0.64450    FALSE
## 1   log.Tlag    WEIGHT 0.65570    FALSE
## 2   log.Tlag   cWEIGHT 0.65570    FALSE
## 4   log.Tlag       SEX 0.79970    FALSE
## 5     log.ka    WEIGHT 0.09485    FALSE
## 7     log.ka   lWEIGHT 0.10770    FALSE
## 8     log.ka       SEX 0.18540    FALSE
## 9      log.V    WEIGHT 0.01699    FALSE
## 10     log.V   cWEIGHT 0.01699    FALSE
## 12     log.V       SEX 0.09431    FALSE
## 15    log.Cl   lWEIGHT 0.39810    FALSE
## 13    log.Cl    WEIGHT 0.43530    FALSE
## 14    log.Cl   cWEIGHT 0.43530    FALSE
## 16    log.Cl       SEX 0.61330    FALSE