<!–html_preserve–><!–/html_preserve–>
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") %>%
add.dosing(dose=10000, nbr.doses=1, dosing.to=2) %>%
add.sampling(seq(0,48,length.out=100));
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): | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Initial Conditions ( $inits): | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
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): | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Initial Conditions ( $inits): | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
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)
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
RxODEcompatible 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") %>%
add.dosing(dose=10000, nbr.doses=1, dosing.to=2) %>%
add.sampling(seq(0,48,length.out=10));
ev2 <- eventTable(amount.units="mg", time.units="hours") %>%
add.dosing(dose=5000, nbr.doses=1, dosing.to=2) %>%
add.sampling(seq(0,48,length.out=8));
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
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 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.