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

Read and plot the data

Let us first read and plot the indomethacin data

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

indomethacin.data <- readDatamlx(data=indomethacin)
print(names(indomethacin.data))
## [1] "treatment" "y"         "id"        "N"
ggplot(indomethacin.data$y, aes(x=time, y=y, color=id)) + geom_point(size=1.5) + geom_line(size=1) +
  theme(legend.position = "none") + scale_y_continuous(trans='log10')


Select the structural model

We use pkbuild to select a PK model for this data, using both clearances and rates as parameterization, comparing linear and nonlinear elimination models:

indomethacin.res <- pkbuild(data=indomethacin, param="both", MM=TRUE, new.dir="indomethacin")
## 
## [1] "indomethacin/pk_VCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
##           OFV           AIC           BIC          BICc standardError 
## -69.626682699 -57.626682699 -58.876125884 -49.284544793   0.005498175 
## 
## [1] "indomethacin/pk_V1V2QCl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
##           OFV           AIC           BIC          BICc standardError 
## -177.29580239 -157.29580239 -159.37820770 -144.99083606    0.02466472 
## 
## [1] "indomethacin/pk_V1V2V3Q2Q3Cl.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
##           OFV           AIC           BIC          BICc standardError 
## -177.24126051 -149.24126051 -152.15662794 -132.97346576    0.06306874 
## 
## [1] "indomethacin/pk_V1V2QVmKm.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
##           OFV           AIC           BIC          BICc standardError 
## -172.27253510 -148.27253510 -150.77142147 -133.98615456    0.04532935 
## 
## [1] "indomethacin/pk_Vk.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
##           OFV           AIC           BIC          BICc standardError 
##  -69.69872326  -57.69872326  -58.94816645  -49.35658536    0.01171416 
## 
## [1] "indomethacin/pk_Vk12k21k.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
##           OFV           AIC           BIC          BICc standardError 
## -174.33641150 -154.33641150 -156.41881681 -142.03144518    0.04406726 
## 
## [1] "indomethacin/pk_Vk12k21k13k31k.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
##           OFV           AIC           BIC          BICc standardError 
## -174.72977760 -146.72977760 -149.64514503 -130.46198285    0.06434624 
## 
## [1] "indomethacin/pk_Vk12k21VmKm.mlxtran"
## Estimation of the population parameters...
## Estimation of the log-likelihood... 
##           OFV           AIC           BIC          BICc standardError 
## -169.27322536 -145.27322536 -147.77211173 -130.98684482    0.07600551

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

print(indomethacin.res$res)
##                                  model        OFV        AIC        BIC
## 1        lib:infusion_2cpt_ClV1QV2.txt -177.29580 -157.29580 -159.37821
## 2       lib:infusion_2cpt_Vkk12k21.txt -174.33641 -154.33641 -156.41882
## 3      lib:infusion_2cpt_V1QV2VmKm.txt -172.27254 -148.27254 -150.77142
## 4   lib:infusion_3cpt_ClV1Q2V2Q3V3.txt -177.24126 -149.24126 -152.15663
## 5    lib:infusion_2cpt_Vk12k21VmKm.txt -169.27323 -145.27323 -147.77211
## 6 lib:infusion_3cpt_Vkk12k21k13k31.txt -174.72978 -146.72978 -149.64515
## 7             lib:infusion_1cpt_Vk.txt  -69.69872  -57.69872  -58.94817
## 8            lib:infusion_1cpt_VCl.txt  -69.62668  -57.62668  -58.87613
##         BICc
## 1 -144.99084
## 2 -142.03145
## 3 -133.98615
## 4 -132.97347
## 5 -130.98684
## 6 -130.46198
## 7  -49.35659
## 8  -49.28454


Using Monolix to fit the data with the selected model

The following Monolix project was created by pkbuild

indomethacin.res$project
## [1] "indomethacin/pk_V1V2QCl.mlxtran"

We can load this project and plot some individual fits obtained with Monolix:

library(lixoftConnectors)
loadProject(indomethacin.res$project)
plotIndividualFits() + scale_y_continuous(trans='log10')