We still use the same theophylline project, where the parameters and the Fisher information matrix have been estimated with Monolix.

We will see how to use the estimated Fisher information matrix for taking into account the uncertainty on the parameter estimates.

project.file <- 'monolixRuns/theophylline_project.mlxtran'

The estimated population parameters are used when Simulx is used with a Monolix project

res1  <- simulx(project = project.file)

print(res1$population, digits=3)
##       ka_pop        V_pop  beta_V_lw70       Cl_pop beta_Cl_lw70 
##       1.5190      33.1129       0.7487       2.8080       0.5141 
##     omega_ka      omega_V     omega_Cl            b 
##       0.6445       0.0944       0.2363       0.1591

and also when several replicates are simulated

res2 <- simulx(project = project.file, 
               nrep    = 4,
               settings=list(seed=123456))
print(res2$population, digits=3)
##       ka_pop        V_pop  beta_V_lw70       Cl_pop beta_Cl_lw70 
##       1.5190      33.1129       0.7487       2.8080       0.5141 
##     omega_ka      omega_V     omega_Cl            b 
##       0.6445       0.0944       0.2363       0.1591
head(res2$y1)
##   rep id time       y1
## 1   1  1 0.25 1.152985
## 2   1  1 0.57 1.951585
## 3   1  1 1.12 4.780470
## 4   1  1 2.02 5.756897
## 5   1  1 3.82 6.256080
## 6   1  1 5.10 5.559285
print(ggplot(data=res2$y1,aes(x=time, y=y1, colour=id)) + 
        geom_point() +  geom_line() + facet_wrap(~ rep, ncol=2))

when npop is defined, then the population parameters are drawn from the distribution of the estimates (derived from the Fisher information matrix)

res3a <- simulx(project = project.file, 
               npop    = 4,
               settings=list(seed=12345))
print(res3a$population, digits=3)
##   pop ka_pop  V_pop beta_V_lw70 Cl_pop beta_Cl_lw70 omega_ka  omega_V
## 1   1 1.6936 33.953     0.68654 2.8840      0.06257  0.77730 0.053053
## 2   2 1.7306 30.973     0.45263 2.9053      0.35734  0.91398 0.188201
## 3   3 1.4863 33.876     0.73265 2.6594      1.06971  0.54807 0.077733
## 4   4 1.3837 32.701     1.23773 2.9724      0.73322  0.45926 0.118682
##   omega_Cl       b
## 1  0.27503 0.18715
## 2  0.21804 0.17334
## 3  0.28871 0.16370
## 4  0.38167 0.16235
head(res3a$y1)
##   pop id time       y1
## 1   1  1 0.25 5.215778
## 2   1  1 0.57 7.948964
## 3   1  1 1.12 8.085736
## 4   1  1 2.02 7.973053
## 5   1  1 3.82 6.075647
## 6   1  1 5.10 5.265295
print(ggplot(data=res3a$y1,aes(x=time, y=y1, colour=id)) + 
        geom_point() +  geom_line() + facet_wrap(~ pop, ncol=2))

If the FIM has been estimated by stochastic approximation, then this matrix is used by default (which is equivalent to set fim="sa").

If the FIM has been estimated by linearization of the model, this matrix can be used by setting fim="lin",

res3b <- simulx(project = project.file, 
                npop    = 4,
                fim     = "lin",
                settings=list(seed=12345))
print(res3b$population, digits=3)
##   pop ka_pop  V_pop beta_V_lw70 Cl_pop beta_Cl_lw70 omega_ka  omega_V
## 1   1 1.6938 33.984     0.67644 2.8857     0.056955  0.76653 0.049636
## 2   2 1.7308 30.957     0.45475 2.9054     0.365272  0.89098 0.192128
## 3   3 1.4863 33.886     0.72612 2.6604     1.067738  0.55839 0.078254
## 4   4 1.3836 32.681     1.26403 2.9684     0.739155  0.45616 0.122249
##   omega_Cl       b
## 1  0.27033 0.18663
## 2  0.22820 0.17469
## 3  0.28263 0.16305
## 4  0.38399 0.16342

Drawing population parameters can also be done by first simulating the population parameters, and then using this data frame with Simulx as parameter.

Remark: the column pop is mandatory.

pop4 <- simpopmlx(n=4, project=project.file)
print(pop4, digits=3)
##   pop ka_pop V_pop beta_V_lw70 Cl_pop beta_Cl_lw70 omega_ka omega_V
## 1   1   1.42  34.5       1.010   2.93        0.587    0.762  0.1020
## 2   2   1.02  29.9       1.087   2.54        0.623    0.793  0.0578
## 3   3   2.05  32.1       0.337   2.69        0.111    1.083  0.1176
## 4   4   1.53  34.3       0.927   3.20        0.788    0.372  0.1690
##   omega_Cl     b
## 1    0.210 0.166
## 2    0.154 0.142
## 3    0.284 0.157
## 4    0.334 0.148
res4 <- simulx(project   = project.file, 
               parameter = pop4)
print(res4$population, digits=3)
##   pop ka_pop V_pop beta_V_lw70 Cl_pop beta_Cl_lw70 omega_ka omega_V
## 1   1   1.42  34.5       1.010   2.93        0.587    0.762  0.1020
## 2   2   1.02  29.9       1.087   2.54        0.623    0.793  0.0578
## 3   3   2.05  32.1       0.337   2.69        0.111    1.083  0.1176
## 4   4   1.53  34.3       0.927   3.20        0.788    0.372  0.1690
##   omega_Cl     b
## 1    0.210 0.166
## 2    0.154 0.142
## 3    0.284 0.157
## 4    0.334 0.148

it is possible to combine the simulation of several populations and several replicates with a given group size,

res5a <- simulx(project = project.file, 
                npop    = 4,
                nrep    = 3,
                group   = list(size=20))
head(res5a$y1)
##   pop rep id time       y1
## 1   1   1  1 0.25 1.698415
## 2   1   1  1 0.57 3.414193
## 3   1   1  1 1.12 5.356088
## 4   1   1  1 2.02 8.914687
## 5   1   1  1 3.82 7.924648
## 6   1   1  1 5.10 6.974387

Set settings$disp.iter=TRUE to display the iteration numbers

res5b <- simulx(project = project.file, 
                npop    = 4,
                nrep    = 3,
                group   = list(size=20),
                settings=list(disp.iter=TRUE))
## 
## population:  1 
## replicate:  1 
## replicate:  2 
## replicate:  3 
## 
## population:  2 
## replicate:  1 
## replicate:  2 
## replicate:  3 
## 
## population:  3 
## replicate:  1 
## replicate:  2 
## replicate:  3 
## 
## population:  4 
## replicate:  1 
## replicate:  2 
## replicate:  3

If the results are stored in a datafile, then a column pop is automatically added,

g1 = list(size=100, treatment=list(amount=5, time=0))
g2 = list(size=100, treatment=list(amount=10, time=0))
res5c <- simulx(project = project.file, 
                npop    = 4,
                group   = list(g1, g2),
                result.file = "theo5c.csv")
head(read.table("theo5c.csv", header=T, sep=","))                
##   pop id time       y group amount
## 1   1  1 0.00       .     1      5
## 2   1  1 0.25 0.04247     1      .
## 3   1  1 0.57 0.08591     1      .
## 4   1  1 1.12 0.09003     1      .
## 5   1  1 2.02 0.11324     1      .
## 6   1  1 3.82 0.11569     1      .