<!–html_preserve–><!–/html_preserve–>

# Simulation of Variability with RxODE

Variability in model parameters can be simulated by creating a matrix of parameter values for use in the simulation. In the example below, 40% variability in clearance is simulated.

``````set.seed(42);
library(RxODE)
mod <- RxODE({
eff(0) = 1
C2 = centr/V2;
C3 = peri/V3;
CL =  TCl*exp(eta.Cl) ## This is coded as a variable in the model
d/dt(depot) =-KA*depot;
d/dt(centr) = KA*depot - CL*C2 - Q*C2 + Q*C3;
d/dt(peri)  =                    Q*C2 - Q*C3;
d/dt(eff)  = Kin - Kout*(1-C2/(EC50+C2))*eff;
})

theta <- c(KA=2.94E-01, TCl=1.86E+01, V2=4.02E+01,  # central
Q=1.05E+01, V3=2.97E+02,                # peripheral
Kin=1, Kout=1, EC50=200)                # effects
``````

Each subproblem can be simulated by using the rxSolve function to run the simulation for each set of parameters of in the parameter matrix.

``````## the column names of the omega matrix need to match the parameters specified by RxODE
omega <- matrix(0.4^2,dimnames=list(NULL,c("eta.Cl")))

ev <- eventTable(amount.units="mg", time.units="hours") %>%

sim  <- rxSolve(mod,theta,ev,omega=omega,nSub=100)

library(ggplot2)

ggplot(sim,aes(time,centr,color=factor(sim.id))) + geom_line(size=1) + coord_trans(y = "log10") + ylab("Central Concentration") +
xlab("Time (hr)") + guides(color=FALSE)
``````

``````ggplot(sim,aes(time,eff,color=factor(sim.id))) + geom_line(size=1) + coord_trans(y = "log10") + ylab("Effect") +
xlab("Time (hr)") + guides(color=FALSE)
``````

It is now straightforward to perform calculations and generate plots with the simulated data. Below, the 5th, 50th, and 95th percentiles of the simulated data are plotted.

``````library(dplyr)

p <- c(0.05, 0.5, 0.95);
s <-sim %>% group_by(time) %>%
do(data.frame(p=p, eff=quantile(.\$eff, probs=p),
eff.n = length(.\$eff), eff.avg = mean(.\$eff),
centr=quantile(.\$centr, probs=p),
centr.n=length(.\$centr),centr.avg = mean(.\$centr))) %>%
mutate(Percentile=factor(sprintf("%d%%",p*100),levels=c("5%","50%","95%")))

ggplot(s,aes(time,centr,color=Percentile)) + geom_line(size=1) + coord_trans(y = "log10") + ylab("Central Concentration") +
xlab("Time (hr)")
``````

``````ggplot(s,aes(time,eff,color=Percentile)) + geom_line(size=1) + ylab("Effect") +
xlab("Time (hr)") + guides(color=FALSE)
``````

Note that you can see the parameters that were simulated for the example

``````head(sim\$param)
``````
``````##   sim.id   V2  V3  TCl      eta.Cl    KA    Q Kin Kout EC50
## 1      1 40.2 297 18.6 -0.23883576 0.294 10.5   1    1  200
## 2      2 40.2 297 18.6  0.37204656 0.294 10.5   1    1  200
## 3      3 40.2 297 18.6 -0.60662015 0.294 10.5   1    1  200
## 4      4 40.2 297 18.6 -0.08291157 0.294 10.5   1    1  200
## 5      5 40.2 297 18.6 -0.18610920 0.294 10.5   1    1  200
## 6      6 40.2 297 18.6 -0.24490319 0.294 10.5   1    1  200
``````

You can also supply a data-frame of parameters to simulate instead of using an omega simulation. In this contrived example we will use the previously simulated data.

``````theta <- sim\$param;
sim  <- rxSolve(mod,theta,ev)
rxHtml(sim)
``````

Solved RxODE object
Parameters (\$params):
sim.id V2 V3 TCl eta.Cl KA Q Kin Kout EC50
1 40.2 297 18.6 -0.2388358 0.294 10.5 1 1 200
2 40.2 297 18.6 0.3720466 0.294 10.5 1 1 200
3 40.2 297 18.6 -0.6066202 0.294 10.5 1 1 200
4 40.2 297 18.6 -0.0829116 0.294 10.5 1 1 200
5 40.2 297 18.6 -0.1861092 0.294 10.5 1 1 200
6 40.2 297 18.6 -0.2449032 0.294 10.5 1 1 200
7 40.2 297 18.6 -0.3878451 0.294 10.5 1 1 200
8 40.2 297 18.6 -0.3590569 0.294 10.5 1 1 200
9 40.2 297 18.6 -0.1473807 0.294 10.5 1 1 200
10 40.2 297 18.6 -0.4905109 0.294 10.5 1 1 200
11 40.2 297 18.6 -0.0369530 0.294 10.5 1 1 200
12 40.2 297 18.6 -0.1401009 0.294 10.5 1 1 200
13 40.2 297 18.6 -0.1195365 0.294 10.5 1 1 200
14 40.2 297 18.6 0.4905827 0.294 10.5 1 1 200
15 40.2 297 18.6 -0.1131729 0.294 10.5 1 1 200
16 40.2 297 18.6 -0.0821916 0.294 10.5 1 1 200
17 40.2 297 18.6 0.0142137 0.294 10.5 1 1 200
18 40.2 297 18.6 -0.0392171 0.294 10.5 1 1 200
19 40.2 297 18.6 0.5021410 0.294 10.5 1 1 200
20 40.2 297 18.6 -0.1727723 0.294 10.5 1 1 200
21 40.2 297 18.6 -0.2967152 0.294 10.5 1 1 200
22 40.2 297 18.6 0.0620091 0.294 10.5 1 1 200
23 40.2 297 18.6 -0.5435836 0.294 10.5 1 1 200
24 40.2 297 18.6 0.5977396 0.294 10.5 1 1 200
25 40.2 297 18.6 -0.0018323 0.294 10.5 1 1 200
26 40.2 297 18.6 0.6313187 0.294 10.5 1 1 200
27 40.2 297 18.6 0.1402706 0.294 10.5 1 1 200
28 40.2 297 18.6 0.2899946 0.294 10.5 1 1 200
29 40.2 297 18.6 0.2129941 0.294 10.5 1 1 200
30 40.2 297 18.6 -0.0857317 0.294 10.5 1 1 200
31 40.2 297 18.6 -0.3140020 0.294 10.5 1 1 200
32 40.2 297 18.6 -0.1207519 0.294 10.5 1 1 200
33 40.2 297 18.6 -0.1179634 0.294 10.5 1 1 200
34 40.2 297 18.6 -0.2244224 0.294 10.5 1 1 200
35 40.2 297 18.6 0.4130835 0.294 10.5 1 1 200
36 40.2 297 18.6 0.1940724 0.294 10.5 1 1 200
37 40.2 297 18.6 0.0251928 0.294 10.5 1 1 200
38 40.2 297 18.6 0.5530819 0.294 10.5 1 1 200
39 40.2 297 18.6 0.1859351 0.294 10.5 1 1 200
40 40.2 297 18.6 -0.4337755 0.294 10.5 1 1 200
41 40.2 297 18.6 -0.0000088 0.294 10.5 1 1 200
42 40.2 297 18.6 -0.6034226 0.294 10.5 1 1 200
43 40.2 297 18.6 0.1519508 0.294 10.5 1 1 200
44 40.2 297 18.6 0.1749234 0.294 10.5 1 1 200
45 40.2 297 18.6 -0.3118630 0.294 10.5 1 1 200
46 40.2 297 18.6 -0.1306790 0.294 10.5 1 1 200
47 40.2 297 18.6 -0.5312595 0.294 10.5 1 1 200
48 40.2 297 18.6 0.4469370 0.294 10.5 1 1 200
49 40.2 297 18.6 -0.2630219 0.294 10.5 1 1 200
50 40.2 297 18.6 -0.3841502 0.294 10.5 1 1 200
51 40.2 297 18.6 -0.2753404 0.294 10.5 1 1 200
52 40.2 297 18.6 0.2000730 0.294 10.5 1 1 200
53 40.2 297 18.6 -0.1045640 0.294 10.5 1 1 200
54 40.2 297 18.6 -0.0569493 0.294 10.5 1 1 200
55 40.2 297 18.6 0.3799748 0.294 10.5 1 1 200
56 40.2 297 18.6 0.3228048 0.294 10.5 1 1 200
57 40.2 297 18.6 0.2941710 0.294 10.5 1 1 200
58 40.2 297 18.6 0.1695037 0.294 10.5 1 1 200
59 40.2 297 18.6 -0.1270347 0.294 10.5 1 1 200
60 40.2 297 18.6 -0.9867963 0.294 10.5 1 1 200
61 40.2 297 18.6 0.1261524 0.294 10.5 1 1 200
62 40.2 297 18.6 -0.0597588 0.294 10.5 1 1 200
63 40.2 297 18.6 -0.6197404 0.294 10.5 1 1 200
64 40.2 297 18.6 -0.0578333 0.294 10.5 1 1 200
65 40.2 297 18.6 -0.2424002 0.294 10.5 1 1 200
66 40.2 297 18.6 -0.4192675 0.294 10.5 1 1 200
67 40.2 297 18.6 0.1207482 0.294 10.5 1 1 200
68 40.2 297 18.6 -0.1301417 0.294 10.5 1 1 200
69 40.2 297 18.6 -0.3905998 0.294 10.5 1 1 200
70 40.2 297 18.6 0.1374133 0.294 10.5 1 1 200
71 40.2 297 18.6 -0.4687107 0.294 10.5 1 1 200
72 40.2 297 18.6 -0.6350925 0.294 10.5 1 1 200
73 40.2 297 18.6 -0.1451134 0.294 10.5 1 1 200
74 40.2 297 18.6 -0.4076332 0.294 10.5 1 1 200
75 40.2 297 18.6 -0.2901100 0.294 10.5 1 1 200
76 40.2 297 18.6 0.1754801 0.294 10.5 1 1 200
77 40.2 297 18.6 0.1047433 0.294 10.5 1 1 200
78 40.2 297 18.6 0.3609997 0.294 10.5 1 1 200
79 40.2 297 18.6 -0.1554073 0.294 10.5 1 1 200
80 40.2 297 18.6 -0.8185079 0.294 10.5 1 1 200
81 40.2 297 18.6 -0.5547246 0.294 10.5 1 1 200
82 40.2 297 18.6 -0.0485547 0.294 10.5 1 1 200
83 40.2 297 18.6 -0.0475025 0.294 10.5 1 1 200
84 40.2 297 18.6 -0.4736681 0.294 10.5 1 1 200
85 40.2 297 18.6 0.0953561 0.294 10.5 1 1 200
86 40.2 297 18.6 0.2578280 0.294 10.5 1 1 200
87 40.2 297 18.6 -0.6064463 0.294 10.5 1 1 200
88 40.2 297 18.6 -0.1607954 0.294 10.5 1 1 200
89 40.2 297 18.6 0.1186228 0.294 10.5 1 1 200
90 40.2 297 18.6 -0.0028003 0.294 10.5 1 1 200
91 40.2 297 18.6 0.3718016 0.294 10.5 1 1 200
92 40.2 297 18.6 0.3931947 0.294 10.5 1 1 200
93 40.2 297 18.6 -0.2489252 0.294 10.5 1 1 200
94 40.2 297 18.6 0.1313894 0.294 10.5 1 1 200
95 40.2 297 18.6 0.4374996 0.294 10.5 1 1 200
96 40.2 297 18.6 -0.2869060 0.294 10.5 1 1 200
97 40.2 297 18.6 -0.3970466 0.294 10.5 1 1 200
98 40.2 297 18.6 0.3146381 0.294 10.5 1 1 200
99 40.2 297 18.6 -0.4496222 0.294 10.5 1 1 200
100 40.2 297 18.6 0.2486597 0.294 10.5 1 1 200
Initial Conditions ( \$inits):
depot centr peri eff
0 0 0 1
First part of data (object):
sim.id time C2 C3 CL depot centr peri eff
1 0.0000000 248.75622 0.000000 14.64832 0 10000.000 0.000 1.000000
1 0.4848485 183.89373 3.646408 14.64832 0 7392.528 1082.983 1.222014
1 0.9696970 136.33885 6.283599 14.64832 0 5480.822 1866.229 1.355838
1 1.4545455 101.46938 8.181449 14.64832 0 4079.069 2429.890 1.415866
1 1.9393939 75.89756 9.537755 14.64832 0 3051.082 2832.713 1.421514
1 2.4242424 57.14041 10.497481 14.64832 0 2297.044 3117.752 1.392573

Even though multiple subjects were simulated, this is still a reactive data frame, meaning you can change things about the model on the fly.

For example, if the effect at time 0 should have been 100, you can fix this by:

``````sim\$eff0 <- 100
rxHtml(sim)
``````

Solved RxODE object
Parameters (\$params):
sim.id V2 V3 TCl eta.Cl KA Q Kin Kout EC50
1 40.2 297 18.6 -0.2388358 0.294 10.5 1 1 200
2 40.2 297 18.6 0.3720466 0.294 10.5 1 1 200
3 40.2 297 18.6 -0.6066202 0.294 10.5 1 1 200
4 40.2 297 18.6 -0.0829116 0.294 10.5 1 1 200
5 40.2 297 18.6 -0.1861092 0.294 10.5 1 1 200
6 40.2 297 18.6 -0.2449032 0.294 10.5 1 1 200
7 40.2 297 18.6 -0.3878451 0.294 10.5 1 1 200
8 40.2 297 18.6 -0.3590569 0.294 10.5 1 1 200
9 40.2 297 18.6 -0.1473807 0.294 10.5 1 1 200
10 40.2 297 18.6 -0.4905109 0.294 10.5 1 1 200
11 40.2 297 18.6 -0.0369530 0.294 10.5 1 1 200
12 40.2 297 18.6 -0.1401009 0.294 10.5 1 1 200
13 40.2 297 18.6 -0.1195365 0.294 10.5 1 1 200
14 40.2 297 18.6 0.4905827 0.294 10.5 1 1 200
15 40.2 297 18.6 -0.1131729 0.294 10.5 1 1 200
16 40.2 297 18.6 -0.0821916 0.294 10.5 1 1 200
17 40.2 297 18.6 0.0142137 0.294 10.5 1 1 200
18 40.2 297 18.6 -0.0392171 0.294 10.5 1 1 200
19 40.2 297 18.6 0.5021410 0.294 10.5 1 1 200
20 40.2 297 18.6 -0.1727723 0.294 10.5 1 1 200
21 40.2 297 18.6 -0.2967152 0.294 10.5 1 1 200
22 40.2 297 18.6 0.0620091 0.294 10.5 1 1 200
23 40.2 297 18.6 -0.5435836 0.294 10.5 1 1 200
24 40.2 297 18.6 0.5977396 0.294 10.5 1 1 200
25 40.2 297 18.6 -0.0018323 0.294 10.5 1 1 200
26 40.2 297 18.6 0.6313187 0.294 10.5 1 1 200
27 40.2 297 18.6 0.1402706 0.294 10.5 1 1 200
28 40.2 297 18.6 0.2899946 0.294 10.5 1 1 200
29 40.2 297 18.6 0.2129941 0.294 10.5 1 1 200
30 40.2 297 18.6 -0.0857317 0.294 10.5 1 1 200
31 40.2 297 18.6 -0.3140020 0.294 10.5 1 1 200
32 40.2 297 18.6 -0.1207519 0.294 10.5 1 1 200
33 40.2 297 18.6 -0.1179634 0.294 10.5 1 1 200
34 40.2 297 18.6 -0.2244224 0.294 10.5 1 1 200
35 40.2 297 18.6 0.4130835 0.294 10.5 1 1 200
36 40.2 297 18.6 0.1940724 0.294 10.5 1 1 200
37 40.2 297 18.6 0.0251928 0.294 10.5 1 1 200
38 40.2 297 18.6 0.5530819 0.294 10.5 1 1 200
39 40.2 297 18.6 0.1859351 0.294 10.5 1 1 200
40 40.2 297 18.6 -0.4337755 0.294 10.5 1 1 200
41 40.2 297 18.6 -0.0000088 0.294 10.5 1 1 200
42 40.2 297 18.6 -0.6034226 0.294 10.5 1 1 200
43 40.2 297 18.6 0.1519508 0.294 10.5 1 1 200
44 40.2 297 18.6 0.1749234 0.294 10.5 1 1 200
45 40.2 297 18.6 -0.3118630 0.294 10.5 1 1 200
46 40.2 297 18.6 -0.1306790 0.294 10.5 1 1 200
47 40.2 297 18.6 -0.5312595 0.294 10.5 1 1 200
48 40.2 297 18.6 0.4469370 0.294 10.5 1 1 200
49 40.2 297 18.6 -0.2630219 0.294 10.5 1 1 200
50 40.2 297 18.6 -0.3841502 0.294 10.5 1 1 200
51 40.2 297 18.6 -0.2753404 0.294 10.5 1 1 200
52 40.2 297 18.6 0.2000730 0.294 10.5 1 1 200
53 40.2 297 18.6 -0.1045640 0.294 10.5 1 1 200
54 40.2 297 18.6 -0.0569493 0.294 10.5 1 1 200
55 40.2 297 18.6 0.3799748 0.294 10.5 1 1 200
56 40.2 297 18.6 0.3228048 0.294 10.5 1 1 200
57 40.2 297 18.6 0.2941710 0.294 10.5 1 1 200
58 40.2 297 18.6 0.1695037 0.294 10.5 1 1 200
59 40.2 297 18.6 -0.1270347 0.294 10.5 1 1 200
60 40.2 297 18.6 -0.9867963 0.294 10.5 1 1 200
61 40.2 297 18.6 0.1261524 0.294 10.5 1 1 200
62 40.2 297 18.6 -0.0597588 0.294 10.5 1 1 200
63 40.2 297 18.6 -0.6197404 0.294 10.5 1 1 200
64 40.2 297 18.6 -0.0578333 0.294 10.5 1 1 200
65 40.2 297 18.6 -0.2424002 0.294 10.5 1 1 200
66 40.2 297 18.6 -0.4192675 0.294 10.5 1 1 200
67 40.2 297 18.6 0.1207482 0.294 10.5 1 1 200
68 40.2 297 18.6 -0.1301417 0.294 10.5 1 1 200
69 40.2 297 18.6 -0.3905998 0.294 10.5 1 1 200
70 40.2 297 18.6 0.1374133 0.294 10.5 1 1 200
71 40.2 297 18.6 -0.4687107 0.294 10.5 1 1 200
72 40.2 297 18.6 -0.6350925 0.294 10.5 1 1 200
73 40.2 297 18.6 -0.1451134 0.294 10.5 1 1 200
74 40.2 297 18.6 -0.4076332 0.294 10.5 1 1 200
75 40.2 297 18.6 -0.2901100 0.294 10.5 1 1 200
76 40.2 297 18.6 0.1754801 0.294 10.5 1 1 200
77 40.2 297 18.6 0.1047433 0.294 10.5 1 1 200
78 40.2 297 18.6 0.3609997 0.294 10.5 1 1 200
79 40.2 297 18.6 -0.1554073 0.294 10.5 1 1 200
80 40.2 297 18.6 -0.8185079 0.294 10.5 1 1 200
81 40.2 297 18.6 -0.5547246 0.294 10.5 1 1 200
82 40.2 297 18.6 -0.0485547 0.294 10.5 1 1 200
83 40.2 297 18.6 -0.0475025 0.294 10.5 1 1 200
84 40.2 297 18.6 -0.4736681 0.294 10.5 1 1 200
85 40.2 297 18.6 0.0953561 0.294 10.5 1 1 200
86 40.2 297 18.6 0.2578280 0.294 10.5 1 1 200
87 40.2 297 18.6 -0.6064463 0.294 10.5 1 1 200
88 40.2 297 18.6 -0.1607954 0.294 10.5 1 1 200
89 40.2 297 18.6 0.1186228 0.294 10.5 1 1 200
90 40.2 297 18.6 -0.0028003 0.294 10.5 1 1 200
91 40.2 297 18.6 0.3718016 0.294 10.5 1 1 200
92 40.2 297 18.6 0.3931947 0.294 10.5 1 1 200
93 40.2 297 18.6 -0.2489252 0.294 10.5 1 1 200
94 40.2 297 18.6 0.1313894 0.294 10.5 1 1 200
95 40.2 297 18.6 0.4374996 0.294 10.5 1 1 200
96 40.2 297 18.6 -0.2869060 0.294 10.5 1 1 200
97 40.2 297 18.6 -0.3970466 0.294 10.5 1 1 200
98 40.2 297 18.6 0.3146381 0.294 10.5 1 1 200
99 40.2 297 18.6 -0.4496222 0.294 10.5 1 1 200
100 40.2 297 18.6 0.2486597 0.294 10.5 1 1 200
Initial Conditions ( \$inits):
depot centr peri eff
0 0 0 100
First part of data (object):
sim.id time C2 C3 CL depot centr peri eff
1 0.0000000 248.75622 0.000000 14.64832 0 10000.000 0.000 100.00000
1 0.4848485 183.89373 3.646408 14.64832 0 7392.528 1082.983 79.54061
1 0.9696970 136.33885 6.283599 14.64832 0 5480.822 1866.229 61.10747
1 1.4545455 101.46938 8.181449 14.64832 0 4079.069 2429.890 45.44976
1 1.9393939 75.89755 9.537755 14.64832 0 3051.082 2832.713 32.86105
1 2.4242424 57.14040 10.497482 14.64832 0 2297.044 3117.752 23.22543

#### Simulation of unexplained variability

In addition to conveniently simulating between subject variability, you can also easily simulate unexplained variability.

``````mod <- RxODE({
eff(0) = 1
C2 = centr/V2;
C3 = peri/V3;
CL =  TCl*exp(eta.Cl) ## This is coded as a variable in the model
d/dt(depot) =-KA*depot;
d/dt(centr) = KA*depot - CL*C2 - Q*C2 + Q*C3;
d/dt(peri)  =                    Q*C2 - Q*C3;
d/dt(eff)  = Kin - Kout*(1-C2/(EC50+C2))*eff;
e = eff+eff.err
cp = centr*(1+cp.err)
})

theta <- c(KA=2.94E-01, TCl=1.86E+01, V2=4.02E+01,  # central
Q=1.05E+01, V3=2.97E+02,                # peripheral
Kin=1, Kout=1, EC50=200)                # effects

sigma <- diag(2)*0.1
dimnames(sigma) <- list(NULL, c("eff.err","cp.err"))

sim  <- rxSolve(mod, theta, ev, omega=omega, nSub=100, sigma=sigma)

p <- c(0.05, 0.5, 0.95);
s <-sim %>% group_by(time) %>%
do(data.frame(p=p, eff=quantile(.\$e, probs=p),
eff.n = length(.\$e), eff.avg = mean(.\$e),
centr=quantile(.\$cp, probs=p),
centr.n=length(.\$cp),centr.avg = mean(.\$cp))) %>%
mutate(Percentile=factor(sprintf("%d%%",p*100),levels=c("5%","50%","95%")))

ggplot(s,aes(time,centr,color=Percentile)) + geom_line(size=1) + coord_trans(y = "log10") + ylab("Central Concentration") +
xlab("Time (hr)")
``````

``````ggplot(s,aes(time,eff,color=Percentile)) + geom_line(size=1) + ylab("Effect") +
xlab("Time (hr)") + guides(color=FALSE)
``````

#### Simulation of Individuals

Sometimes you may want to match the dosing and observations of individuals in a clinical trial. To do this you will have to create a data.frame using the `RxODE` event specification as well as an `ID` column to indicate an individual. The RxODE event vignette talks more about how these datasets should be created.

If you have a NONMEM/Monlix dataset and the package `nlmixr`, you can convert the `NONMEM``dataset to a`RxODE```compatible dataset to use for simulation with the```nmDataConvert function.

Instead of using nlmixr for this simple example, I will combine two RxODE event tables.

``````ev1 <- eventTable(amount.units="mg", time.units="hours") %>%

ev2 <- eventTable(amount.units="mg", time.units="hours") %>%

dat <- rbind(data.frame(ID=1, ev1\$get.EventTable()),
data.frame(ID=2, ev2\$get.EventTable()))

## Note the number of subject is not needed since it is determined by the data
sim  <- rxSolve(mod, theta, dat, omega=omega, sigma=sigma)

sim %>% select(id, time, e, cp)
``````
``````##    id      time         e          cp
## 1   1  0.000000 1.1062517 4456.646344
## 2   1  5.333333 1.3497727  341.324470
## 3   1 10.666667 1.0046824  193.351840
## 4   1 16.000000 0.9320987   54.932311
## 5   1 21.333333 0.9184267   62.081461
## 6   1 26.666667 0.8216998   21.285549
## 7   1 32.000000 0.5856710   89.440303
## 8   1 37.333333 0.9387878   62.370854
## 9   1 42.666667 0.9791421   78.795529
## 10  1 48.000000 0.8765182   94.805161
## 11  2  0.000000 1.1109873 4643.837151
## 12  2  6.857143 0.9887559   75.655019
## 13  2 13.714286 1.1978570   23.609678
## 14  2 20.571429 1.4758890   54.211423
## 15  2 27.428571 1.5740115   30.579824
## 16  2 34.285714 0.4820187   20.341122
## 17  2 41.142857 0.9270456   10.241666
## 18  2 48.000000 0.9148543    5.287749
``````

#### Simulation of Clinical Trials

By either using a simple single event table, or data from a clinical trial as described above, a complete clinical trial simulation can be performed.

Typically in clinical trial simulations you want to account for the uncertainty in the fixed parameter estimates, and even the uncertainty in both your between subject variability as well as the unexplained variability.

`RxODE` allows you to account for these uncertainties by simulating multiple virtual “studies,” specified by the parameter `nStud`. In a single virtual study:

• A Population effect parameter is sampled from a multivariate normal distribution with mean given by the parameter estimates and the variance specified by the named matrix `thetaMat`.

• A between subject variability/covariance matrix is sampled from either a scaled inverse chi-squared distribution (for the univariate case) or a inverse Wishart that is parameterized to scale to the conjugate prior covariance term, as described by the wikipedia article. (This is not the same as the scaled inverse Wishart distribution ). In the case of the between subject variability, the variance/covariance matrix is given by the 'omega' matrix parameter and the degrees of freedom is the number of subjects in the simulation.

• Unexplained variability is also simulated from the scaled inverse chi squared distribution or inverse Wishart distribution with the variance/covariance matrix given by the 'sigma' matrix parameter and the degrees of freedom given by the number of observations being simulated.

The covariance/variance prior is simulated from `RxODE`s `cvPost` function.

An example of this simulation is below:

``````## Creating covariance matrix
tmp <- matrix(rnorm(8^2), 8, 8)
tMat <- tcrossprod(tmp, tmp) / (8 ^ 2)
dimnames(tMat) <- list(NULL, names(theta))

sim  <- rxSolve(mod, theta, ev, omega=omega, nSub=100, sigma=sigma, thetaMat=tMat, nStud=10,
dfSub=10, dfObs=100)

p <- c(0.05, 0.5, 0.95);
s <-sim %>% group_by(time) %>%
do(data.frame(p=p, eff=quantile(.\$e, probs=p),
eff.n = length(.\$e), eff.avg = mean(.\$e),
centr=quantile(.\$cp, probs=p),
centr.n=length(.\$cp),centr.avg = mean(.\$cp))) %>%
mutate(Percentile=factor(sprintf("%d%%",p*100),levels=c("5%","50%","95%")))

ggplot(s,aes(time,centr,color=Percentile)) + geom_line(size=1) + coord_trans(y = "log10") + ylab("Central Concentration") +
xlab("Time (hr)")
``````

``````ggplot(s,aes(time,eff,color=Percentile)) + geom_line(size=1) + ylab("Effect") +
xlab("Time (hr)") + guides(color=FALSE)
``````

If you wish you can see what `omega` and `sigma` was used for each virtual study by accessing them in the solved data object with `\$omega.list` and `\$sigma.list`:

``````head(sim\$omega.list)
``````
``````## [[1]]
##           [,1]
## [1,] 0.1147344
##
## [[2]]
##           [,1]
## [1,] 0.1200283
##
## [[3]]
##           [,1]
## [1,] 0.1298885
##
## [[4]]
##           [,1]
## [1,] 0.3125934
##
## [[5]]
##           [,1]
## [1,] 0.1273521
##
## [[6]]
##           [,1]
## [1,] 0.1207694
``````
``````head(sim\$sigma.list)
``````
``````## [[1]]
##              [,1]         [,2]
## [1,]  0.115055008 -0.001300522
## [2,] -0.001300522  0.095544021
##
## [[2]]
##             [,1]        [,2]
## [1,] 0.115403508 0.009912083
## [2,] 0.009912083 0.103826250
##
## [[3]]
##              [,1]         [,2]
## [1,]  0.104344964 -0.002497709
## [2,] -0.002497709  0.104757997
##
## [[4]]
##             [,1]        [,2]
## [1,] 0.091556183 0.005414474
## [2,] 0.005414474 0.087806016
##
## [[5]]
##             [,1]        [,2]
## [1,] 0.092347188 0.004894803
## [2,] 0.004894803 0.091091435
##
## [[6]]
##               [,1]          [,2]
## [1,]  0.0906104390 -0.0009193655
## [2,] -0.0009193655  0.1207593030
``````

You can also see the parameter realizations from the `\$params` data frame.

If you do not wish to sample from the prior distributions of either the `omega` or `sigma` matrices, you can turn off this feature by specifying the `simVariability = FALSE` option when solving:

``````sim  <- rxSolve(mod, theta, ev, omega=omega, nSub=100, sigma=sigma, thetaMat=tMat, nStud=10,
simVariability=FALSE);

p <- c(0.05, 0.5, 0.95);
s <-sim %>% group_by(time) %>%
do(data.frame(p=p, eff=quantile(.\$e, probs=p),
eff.n = length(.\$e), eff.avg = mean(.\$e),
centr=quantile(.\$cp, probs=p),
centr.n=length(.\$cp),centr.avg = mean(.\$cp))) %>%
mutate(Percentile=factor(sprintf("%d%%",p*100),levels=c("5%","50%","95%")))

ggplot(s,aes(time,centr,color=Percentile)) + geom_line(size=1) + coord_trans(y = "log10") + ylab("Central Concentration") +
xlab("Time (hr)")
``````

``````ggplot(s,aes(time,eff,color=Percentile)) + geom_line(size=1) + ylab("Effect") +
xlab("Time (hr)") + guides(color=FALSE)
``````

Note since realizations of `omega` and `sigma` were not simulated, `\$omega.list` and `\$sigma.list` both return `NULL`.

#### RxODE multi-threaded solving and simulation

RxODE now supports multi-threaded solving on OpenMP supported compilers, including linux and windows. Mac OSX can also be supported but takes additional setup. By default it uses all your available cores for solving as determined by `rxCores()`. This may be overkill depending on your system, at a certain point the speed of solving is limited by things other than computing power.

You can also speed up simulation by using the multi-cores to generate random deviates with `mvnfast`. This is controlled by the `nCoresRV` parameter. For example:

``````sim  <- rxSolve(mod, theta, ev, omega=omega, nSub=100, sigma=sigma, thetaMat=tMat, nStud=10,
nCoresRV=2);

p <- c(0.05, 0.5, 0.95);
s <-sim %>% group_by(time) %>%
do(data.frame(p=p, eff=quantile(.\$e, probs=p),
eff.n = length(.\$e), eff.avg = mean(.\$e),
centr=quantile(.\$cp, probs=p),
centr.n=length(.\$cp),centr.avg = mean(.\$cp))) %>%
mutate(Percentile=factor(sprintf("%d%%",p*100),levels=c("5%","50%","95%")))

ggplot(s,aes(time,centr,color=Percentile)) + geom_line(size=1) + coord_trans(y = "log10") + ylab("Central Concentration") +
xlab("Time (hr)")
``````

``````ggplot(s,aes(time,eff,color=Percentile)) + geom_line(size=1) + ylab("Effect") +
xlab("Time (hr)") + guides(color=FALSE)
``````

The default for this is `1` core since the result depends on the number of cores and the random seed you use in your simulation. However, you can always speed up this process with more cores if you are sure your collaborators have the same number of cores available to them and have OpenMP thread-capable compile.