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 <- 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... 
##           OFV           AIC           BIC          BICc standardError 
##  345.55005921  361.55005921  365.42931241  376.94223787    0.03606701 
## 
## [1] "theophylline/pk_Tk0VCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
##           OFV           AIC           BIC          BICc standardError 
##  301.53734841  317.53734841  321.41660161  332.92952708    0.02679921 
## 
## [1] "theophylline/pk_TlagkaVCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
##           OFV           AIC           BIC          BICc standardError 
##  278.22725507  298.22725507  303.07632156  316.89183212    0.09089036 
## 
## [1] "theophylline/pk_TlagTk0VCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
##           OFV           AIC           BIC          BICc standardError 
##      286.5850      306.5850      311.4340      325.2496        0.0643 
## 
## [1] "theophylline/pk_TlagkaV1V2QCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
##           OFV           AIC           BIC          BICc standardError 
##  278.02092244  306.02092244  312.80961554  331.23029628    0.09770697 
## 
## [1] "theophylline/pk_TlagTk0V1V2QCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
##           OFV           AIC           BIC          BICc standardError 
##  285.85827498  313.85827498  320.64696807  339.06764882    0.09648946

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

print(theo.clearance$res)
##                               model      OFV      AIC      BIC     BICc
## 1      lib:oral1_1cpt_TlagkaVCl.txt 278.2273 298.2273 303.0763 316.8918
## 2     lib:oral0_1cpt_TlagTk0VCl.txt 286.5850 306.5850 311.4340 325.2496
## 3  lib:oral1_2cpt_TlagkaClV1QV2.txt 278.0209 306.0209 312.8096 331.2303
## 4         lib:oral0_1cpt_Tk0VCl.txt 301.5373 317.5373 321.4166 332.9295
## 5 lib:oral0_2cpt_TlagTk0ClV1QV2.txt 285.8583 313.8583 320.6470 339.0676
## 6          lib:oral1_1cpt_kaVCl.txt 345.5501 361.5501 365.4293 376.9422


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")
## 
## --------------------------------------------------
## 
## Building:
##    -  The covariate model
##    -  The correlation model
##    -  The residual error model
##  
## __________________________________________________
## - - - Initialization - - -
## 
## Covariate model:
##      SEX WEIGHT
## Tlag   0      0
## ka     0      0
## V      0      0
## Cl     0      0
## 
## Correlation model:
## [1] "NULL"
## 
## Residual error model:
##        CONC 
## "combined2" 
## Sampling of the conditional distribution using the initial model ... 
## 
## Estimated criteria (importanceSampling):
##    AIC    BIC   BICc   s.e. 
## 298.23 303.08 316.89   0.09 
## 
## Estimation of the population parameters using the transformed covariates ... 
## Sampling of the conditional distribution using the the transformed covariates ... 
## __________________________________________________
## - - - Iteration 1 - - -
## 
## Covariate model:
##      SEX WEIGHT logtWEIGHT
## Tlag   0      0          0
## ka     0      0          1
## V      0      0          1
## Cl     0      0          0
## 
## Correlation model:
## [[1]]
## [1] "Cl" "V" 
## 
## 
## 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):
##    AIC    BIC   BICc   s.e. 
## 281.34 287.15 298.67   0.11 
## __________________________________________________
## - - - Iteration 2 - - -
## 
## Covariate model:
##      SEX WEIGHT logtWEIGHT
## Tlag   0      0          0
## ka     0      0          1
## V      0      0          1
## Cl     0      0          1
## 
## Correlation model:
## [[1]]
## [1] "Cl" "V" 
## 
## 
## 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):
##    AIC    BIC   BICc   s.e. 
## 282.69 289.00 300.51   0.09 
## __________________________________________________
## - - - Iteration 3 - - -
## 
## No difference between two successive iterations
## __________________________________________________
## - - - Further tests - - -
## _______________________
## Add parameters/covariates relationships:
##   parameter covariate p.value
## 4         V       SEX 0.02161
## 
## Run scenario for model 4 ... 
## Estimation of the population parameters... 
## _______________________
## Remove parameters/covariates relationships:
##          coefficient  p.value
## 1 beta_ka_logtWEIGHT 0.064731
## 3       beta_V_SEX_M 0.475142
## 
## Run scenario for model 5 ... 
## Estimation of the population parameters... 
## Sampling from the conditional distribution... 
## Estimation of the log-likelihood... 
## 
## Estimated criteria (importanceSampling):
##    AIC    BIC   BICc   s.e. 
## 282.75 288.08 299.60   0.10 
## _______________________
## Remove parameters/covariates relationships:
##          coefficient   p.value
## 1 beta_ka_logtWEIGHT 0.0646742
## 3 beta_ka_logtWEIGHT  0.067152
## 
## Run scenario for model 6 ... 
## Estimation of the population parameters... 
## Sampling from the conditional distribution... 
## Estimation of the log-likelihood... 
## 
## Estimated criteria (importanceSampling):
##    AIC    BIC   BICc   s.e. 
## 282.75 288.08 299.60   0.10 
## __________________________________________________
## 
## Final statistical model:
## 
## Covariate model:
##      logtWEIGHT SEX WEIGHT
## Tlag          0   0      0
## ka            1   0      0
## V             1   0      0
## Cl            0   0      0
## 
## Correlation model:
## [[1]]
## [1] "Cl" "V" 
## 
## 
## Residual error model:
##           CONC 
## "proportional" 
## 
## Estimated criteria (importanceSampling):
##    AIC    BIC   BICc   s.e. 
## 281.34 287.15 298.67   0.11 
## 
## total time: 88.7s
## __________________________________________________

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... 
##           OFV           AIC           BIC          BICc standardError 
##   337.7808318   353.7808318   357.6600850   369.1730104     0.0547343 
## 
## [1] "theophylline/pk_Tk0Vk.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
##           OFV           AIC           BIC          BICc standardError 
##  292.43897454  308.43897454  312.31822774  323.83115321    0.03424432 
## 
## [1] "theophylline/pk_TlagkaVk.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
##           OFV           AIC           BIC          BICc standardError 
##   269.2928020   289.2928020   294.1418685   307.9573791     0.1346861 
## 
## [1] "theophylline/pk_TlagTk0Vk.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
##           OFV           AIC           BIC          BICc standardError 
##  276.67070413  296.67070413  301.51977063  315.33528119    0.07301283 
## 
## [1] "theophylline/pk_TlagkaVk12k21k.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
##           OFV           AIC           BIC          BICc standardError 
##   260.3715933   288.3715933   295.1602864   313.5809671     0.1140874 
## 
## [1] "theophylline/pk_TlagTk0Vk12k21k.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
##           OFV           AIC           BIC          BICc standardError 
##   276.6522752   304.6522752   311.4409683   329.8616490     0.1150421

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

print(theo.rate$res)
##                                model      OFV      AIC      BIC     BICc
## 1        lib:oral1_1cpt_TlagkaVk.txt 269.2928 289.2928 294.1419 307.9574
## 2  lib:oral1_2cpt_TlagkaVkk12k21.txt 260.3716 288.3716 295.1603 313.5810
## 3       lib:oral0_1cpt_TlagTk0Vk.txt 276.6707 296.6707 301.5198 315.3353
## 4           lib:oral0_1cpt_Tk0Vk.txt 292.4390 308.4390 312.3182 323.8312
## 5 lib:oral0_2cpt_TlagTk0Vkk12k21.txt 276.6523 304.6523 311.4410 329.8616
## 6            lib:oral1_1cpt_kaVk.txt 337.7808 353.7808 357.6601 369.1730

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

Building now the statistical model leads to the same result

res.rate <- buildmlx(project=theo.rate$project, covToTransform = "all")
## 
## --------------------------------------------------
## 
## Building:
##    -  The covariate model
##    -  The correlation model
##    -  The residual error model
##  
## __________________________________________________
## - - - Initialization - - -
## 
## Covariate model:
##      SEX WEIGHT
## Tlag   0      0
## ka     0      0
## V      0      0
## k      0      0
## 
## Correlation model:
## [1] "NULL"
## 
## Residual error model:
##        CONC 
## "combined2" 
## Sampling of the conditional distribution using the initial model ... 
## 
## Estimated criteria (importanceSampling):
##    AIC    BIC   BICc   s.e. 
## 289.29 294.14 307.96   0.13 
## 
## Estimation of the population parameters using the transformed covariates ... 
## Sampling of the conditional distribution using the the transformed covariates ... 
## __________________________________________________
## - - - Iteration 1 - - -
## 
## Covariate model:
##      SEX WEIGHT logtWEIGHT
## Tlag   0      0          0
## ka     0      0          1
## V      0      0          1
## k      0      0          0
## 
## Correlation model:
## [[1]]
## [1] "k" "V"
## 
## 
## 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):
##    AIC    BIC   BICc   s.e. 
## 281.09 286.91 298.42   0.11 
## __________________________________________________
## - - - Iteration 2 - - -
## 
## No difference between two successive iterations
## __________________________________________________
## - - - Further tests - - -
## _______________________
## Add parameters/covariates relationships:
##   parameter covariate p.value
## 4         V       SEX 0.02609
## 
## Run scenario for model 3 ... 
## Estimation of the population parameters... 
## _______________________
## Remove parameters/covariates relationships:
##          coefficient   p.value
## 1 beta_ka_logtWEIGHT 0.0641787
## 3       beta_V_SEX_M  0.405099
## 
## Run scenario for model 4 ... 
## Estimation of the population parameters... 
## Sampling from the conditional distribution... 
## Estimation of the log-likelihood... 
## 
## Estimated criteria (importanceSampling):
##    AIC    BIC   BICc   s.e. 
## 282.22 287.55 299.07   0.10 
## _______________________
## Remove parameters/covariates relationships:
##          coefficient   p.value
## 1 beta_ka_logtWEIGHT 0.0650827
## 3 beta_ka_logtWEIGHT 0.0641303
## 
## Run scenario for model 5 ... 
## Estimation of the population parameters... 
## Sampling from the conditional distribution... 
## Estimation of the log-likelihood... 
## 
## Estimated criteria (importanceSampling):
##    AIC    BIC   BICc   s.e. 
## 282.22 287.55 299.07   0.10 
## __________________________________________________
## 
## Final statistical model:
## 
## Covariate model:
##      logtWEIGHT SEX WEIGHT
## Tlag          0   0      0
## ka            1   0      0
## V             1   0      0
## k             0   0      0
## 
## Correlation model:
## [[1]]
## [1] "k" "V"
## 
## 
## Residual error model:
##           CONC 
## "proportional" 
## 
## Estimated criteria (importanceSampling):
##    AIC    BIC   BICc   s.e. 
## 281.09 286.91 298.42   0.11 
## 
## total time: 35.9s
## __________________________________________________

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)
print(theo.fim$confint)
##                       estimate       lower       upper
## Tlag_pop            0.11631132  0.07848992  0.17235746
## ka_pop              2.51388195  1.68150864  3.75829320
## beta_ka_logtWEIGHT  2.95943632 -0.18014571  6.09901835
## V_pop              33.30467284 30.91644263 35.87738881
## beta_V_logtWEIGHT   0.62156942  0.27002129  0.97311755
## Cl_pop              2.77872701  2.39962018  3.21772748
## omega_Tlag          0.58568795  0.35493024  0.96647265
## omega_ka            0.68732363  0.45197105  1.04523015
## omega_V             0.12264342  0.07886405  0.19072579
## omega_Cl            0.25781764  0.17357407  0.38294851
## corr_V_Cl           0.84225218  0.54226640  0.95178229
## b                   0.08370576  0.07033487  0.09707664

0 belongs to the 95%CI for \(\beta_{ka, {\rm 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.ttest    p.lrt p.wald_lin   p.wald_SA in.model
## 6     log.ka logtWEIGHT 0.09853 0.063700 0.06467420 0.067152000     TRUE
## 9      log.V logtWEIGHT 0.01508 0.006409 0.00052945 0.000718617     TRUE
## 12    log.Cl logtWEIGHT 0.39160 0.336000         NA          NA    FALSE
## 11    log.Cl     WEIGHT 0.42870 0.374100         NA          NA    FALSE
## 10    log.Cl        SEX 0.59820 0.554000         NA          NA    FALSE
## 5     log.ka     WEIGHT 0.08634 0.054240         NA          NA    FALSE
## 4     log.ka        SEX 0.17970 0.132000         NA          NA    FALSE
## 3   log.Tlag logtWEIGHT 0.68260 0.646000         NA          NA    FALSE
## 2   log.Tlag     WEIGHT 0.69320 0.657700         NA          NA    FALSE
## 1   log.Tlag        SEX 0.82350 0.802200         NA          NA    FALSE
## 8      log.V     WEIGHT 0.01800 0.007965         NA          NA    FALSE
## 7      log.V        SEX 0.10000 0.064890         NA          NA    FALSE