Chapter 11 Simulation
11.1 Single Subject solving
Originally, RxODE was only created to solve ODEs for one individual. That is a single system without any changes in individual parameters.
Of course this is still supported, the classic examples are found in RxODE intro.
This article discusses the differences between multiple subject and single subject solving. There are three differences:
- Single solving does not solve each ID in parallel
- Single solving lacks the
idcolumn in parameters($params) as well as in the actual dataset. - Single solving allows parameter exploration easier because each parameter can be modified. With multiple subject solves, you have to make sure to update each individual parameter.
The first obvious difference is in speed; With multiple subjects you can run each subject ID in parallel. For more information and examples of the speed gains with multiple subject solving see the Speeding up RxODE vignette.
The next difference is the amount of information output in the final data.
Taking the 2 compartment indirect response model originally in the tutorial:
library(RxODE)
mod1 <-RxODE({
KA=2.94E-01
CL=1.86E+01
V2=4.02E+01
Q=1.05E+01
V3=2.97E+02
Kin=1
Kout=1
EC50=200
C2 = centr/V2
C3 = peri/V3
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
eff(0) = 1
})
et <- et(amount.units='mg', time.units='hours') %>%
et(dose=10000, addl=9, ii=12) %>%
et(amt=20000, nbr.doses=5, start.time=120, dosing.interval=24) %>%
et(0:240) # samplingNow a simple solve
x <- rxSolve(mod1, et)
x#> ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ Solved RxODE object ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂
#> ── Parameters (x$params): ──────────────────────────────────
#> KA CL V2 Q V3 Kin Kout EC50
#> 0.294 18.600 40.200 10.500 297.000 1.000 1.000 200.000
#> ── Initial Conditions (x$inits): ───────────────────────────
#> depot centr peri eff
#> 0 0 0 1
#> ── First part of data (object): ────────────────────────────
#> # A tibble: 241 x 7
#> time C2 C3 depot centr peri eff
#> [h] <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 0 0 0 10000 0 0 1
#> 2 1 44.4 0.920 7453. 1784. 273. 1.08
#> 3 2 54.9 2.67 5554. 2206. 794. 1.18
#> 4 3 51.9 4.46 4140. 2087. 1324. 1.23
#> 5 4 44.5 5.98 3085. 1789. 1776. 1.23
#> 6 5 36.5 7.18 2299. 1467. 2132. 1.21
#> # … with 235 more rows
#> ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂print(x)#> ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ Solved RxODE object ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂
#> ── Parameters ($params): ───────────────────────────────────
#> KA CL V2 Q V3 Kin Kout EC50
#> 0.294 18.600 40.200 10.500 297.000 1.000 1.000 200.000
#> ── Initial Conditions ($inits): ────────────────────────────
#> depot centr peri eff
#> 0 0 0 1
#> ── First part of data (object): ────────────────────────────
#> # A tibble: 241 x 7
#> time C2 C3 depot centr peri eff
#> [h] <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 0 0 0 10000 0 0 1
#> 2 1 44.4 0.920 7453. 1784. 273. 1.08
#> 3 2 54.9 2.67 5554. 2206. 794. 1.18
#> 4 3 51.9 4.46 4140. 2087. 1324. 1.23
#> 5 4 44.5 5.98 3085. 1789. 1776. 1.23
#> 6 5 36.5 7.18 2299. 1467. 2132. 1.21
#> # … with 235 more rows
#> ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂plot(x, C2, eff)
To better see the differences between the single solve, you can solve for 2 individuals
x2 <- rxSolve(mod1, et %>% et(id=1:2), params=data.frame(CL=c(18.6, 7.6)))#> Warning: 'ID' missing in 'parameters' dataset
#> individual parameters are assumed to have the same order as the event datasetprint(x2)#> ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ Solved RxODE object ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂
#> ── Parameters ($params): ───────────────────────────────────
#> # A tibble: 2 x 9
#> id KA CL V2 Q V3 Kin Kout EC50
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0.294 18.6 40.2 10.5 297 1 1 200
#> 2 2 0.294 7.6 40.2 10.5 297 1 1 200
#> ── Initial Conditions ($inits): ────────────────────────────
#> depot centr peri eff
#> 0 0 0 1
#> ── First part of data (object): ────────────────────────────
#> # A tibble: 482 x 8
#> id time C2 C3 depot centr peri eff
#> <int> [h] <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0 0 0 10000 0 0 1
#> 2 1 1 44.4 0.920 7453. 1784. 273. 1.08
#> 3 1 2 54.9 2.67 5554. 2206. 794. 1.18
#> 4 1 3 51.9 4.46 4140. 2087. 1324. 1.23
#> 5 1 4 44.5 5.98 3085. 1789. 1776. 1.23
#> 6 1 5 36.5 7.18 2299. 1467. 2132. 1.21
#> # … with 476 more rows
#> ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂plot(x2, C2, eff)
By observing the two solves, you can see:
- A multiple subject solve contains the
idcolumn both in the data frame and then data frame of parameters for each subject.
The last feature that is not as obvious, modifying the individual parameters. For single subject data, you can modify the RxODE data frame changing initial conditions and parameter values as if they were part of the data frame, as described in the RxODE Data Frames.
For multiple subject solving, this feature still works, but requires care when supplying each individual’s parameter value, otherwise you may change the solve and drop parameter for key individuals.
11.1.1 Summary of Single solve vs Multiple subject solving
| Feature | Single Subject Solve | Multiple Subject Solve |
|---|---|---|
| Parallel | None | Each Subject |
| $params | data.frame with one parameter value | data.frame with one parameter per subject (w/ID column) |
| solved data | Can modify individual parameters with $ syntax | Have to modify all the parameters to update solved object |
11.2 Population Simulations with RxODE
11.2.1 Simulation of Variability with RxODE
In pharmacometrics the nonlinear-mixed effect modeling software (like nlmixr) characterizes the between-subject variability. With this between subject variability you can simulate new subjects.
Assuming that you have a 2-compartment, indirect response model, you can set create an RxODE model describing this system below:
11.2.1.1 Setting up the RxODE model
library(RxODE)
set.seed(32)
mod <- RxODE({
eff(0) = 1
C2 = centr/V2*(1+prop.err);
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;
})11.2.1.2 Adding the parameter estimates
The next step is to get the parameters into R so that you can start the simulation:
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, prop.err=0) # effectsIn this case, I use lotri to specify the omega since it uses similar
lower-triangular matrix specification as nlmixr (also similar to
NONMEM):
### the column names of the omega matrix need to match the parameters specified by RxODE
omega <- lotri(eta.Cl ~ 0.4^2)
omega#> eta.Cl
#> eta.Cl 0.1611.2.1.3 Simulating
The next step to simulate is to create the dosing regimen for overall simulation:
ev <- et(amount.units="mg", time.units="hours") %>%
et(amt=10000, cmt="centr")If you wish, you can also add sampling times (though now RxODE can fill these in for you):
ev <- ev %>% et(0,48, length.out=100)Note the et takes similar arguments as seq when adding sampling
times. There are more methods to adding sampling times and events to
make complex dosing regimens (See the event
vignette). This includes ways to add variability
to the both the sampling and dosing
times).
Once this is complete you can simulate using the rxSolve routine:
sim <- rxSolve(mod,theta,ev,omega=omega,nSub=100)To quickly look and customize your simulation you use the default
plot routine. Since this is an RxODE object, it will create a
ggplot2 object that you can modify as you wish. The extra parameter
to the plot tells RxODE/R what piece of information you are
interested in plotting. In this case, we are interested in looking at
the derived parameter C2:
11.2.1.4 Checking the simulation with plot
library(ggplot2)
### The plots from RxODE are ggplots so they can be modified with
### standard ggplot commands.
plot(sim, C2, log="y") +
ylab("Central Compartment") 
Of course this additional parameter could also be a state value, like eff:
### They also takes many of the standard plot arguments; See ?plot
plot(sim, eff, ylab="Effect")
Or you could even look at the two side-by-side:
plot(sim, C2, eff)
Or stack them with patchwork
library(patchwork)
plot(sim, C2, log="y") / plot(sim, eff)
11.2.1.5 Processing the data to create summary plots
Usually in pharmacometric simulations it is not enough to simply simulate the system. We have to do something easier to digest, like look at the central and extreme tendencies of the simulation.
Since the RxODE solve object is a type of data
frame
It is now straightforward to perform calculations and generate plots with the simulated data. You can
Below, the 5th, 50th, and 95th percentiles of the simulated data are plotted.
confint(sim, "C2", level=0.95) %>%
plot(ylab="Central Concentration", log="y")#> ! in order to put confidence bands around the intervals, you need at least 2500 simulations#> summarizing data...done
confint(sim, "eff", level=0.95) %>%
plot(ylab="Effect")#> ! in order to put confidence bands around the intervals, you need at least 2500 simulations#> summarizing data...done
Note that you can see the parameters that were simulated for the example
head(sim$param)#> sim.id V2 prop.err V3 TCl eta.Cl KA Q Kin Kout EC50
#> 1 1 40.2 0 297 18.6 0.2368417 0.294 10.5 1 1 200
#> 2 2 40.2 0 297 18.6 0.5454099 0.294 10.5 1 1 200
#> 3 3 40.2 0 297 18.6 0.1828379 0.294 10.5 1 1 200
#> 4 4 40.2 0 297 18.6 -0.2237885 0.294 10.5 1 1 200
#> 5 5 40.2 0 297 18.6 0.4640872 0.294 10.5 1 1 200
#> 6 6 40.2 0 297 18.6 -0.2748536 0.294 10.5 1 1 20011.2.1.6 Simulation of unexplained variability (sigma)
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 <- lotri(eff.err ~ 0.1, cp.err ~ 0.1)
sim <- rxSolve(mod, theta, ev, omega=omega, nSub=100, sigma=sigma)
s <- confint(sim, c("eff", "centr"));#> ! in order to put confidence bands around the intervals, you need at least 2500 simulations#> summarizing data...doneplot(s)
11.2.1.7 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.
library(dplyr)
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 [h] 1.0563542 11329.59098
#> 2 1 5.333333 [h] 1.4003578 376.07820
#> 3 1 10.666667 [h] 0.0510544 117.09167
#> 4 1 16.000000 [h] 1.4589483 141.30089
#> 5 1 21.333333 [h] 1.1416624 84.85403
#> 6 1 26.666667 [h] 1.2504412 95.93320
#> 7 1 32.000000 [h] 0.9425509 144.84771
#> 8 1 37.333333 [h] 1.5173332 148.73731
#> 9 1 42.666667 [h] 1.2391798 60.69626
#> 10 1 48.000000 [h] 1.3173971 53.60546
#> 11 2 0.000000 [h] 0.7351683 5471.03043
#> 12 2 6.857143 [h] 0.7138482 109.19130
#> 13 2 13.714286 [h] 1.2041123 137.81498
#> 14 2 20.571429 [h] 1.1766657 81.08167
#> 15 2 27.428571 [h] 1.7274978 57.74205
#> 16 2 34.285714 [h] 0.4546936 60.74535
#> 17 2 41.142857 [h] 0.7159257 44.59950
#> 18 2 48.000000 [h] 1.3206859 42.0386011.3 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. Each
of these studies samples a realization of fixed effect parameters and
covariance matrices for the between subject variability(omega) and
unexplained variabilities (sigma). Depending on the information you
have from the models, there are a few strategies for simulating a
realization of the omega and sigma matrices.
The first strategy occurs when either there is not any standard errors for standard deviations (or related parameters), or there is a modeled correlation in the model you are simulating from. In that case the suggested strategy is to use the inverse Wishart (parameterized to scale to the conjugate prior)/scaled inverse chi distribution. this approach uses a single parameter to inform the variability of the covariance matrix sampled (the degrees of freedom).
The second strategy occurs if you have standard errors on the
variance/standard deviation with no modeled correlations in the
covariance matrix. In this approach you perform separate simulations
for the standard deviations and the correlation matrix. First you
simulate the variance/standard deviation components in the thetaMat
multivariate normal simulation. After simulation and transformation
to standard deviations, a correlation matrix is simulated using the
degrees of freedom of your covariance matrix. Combining the simulated
standard deviation with the simulated correlation matrix will give a
simulated covariance matrix. For smaller dimension covariance matrices
(dimension < 10x10) it is recommended you use the lkj distribution
to simulate the correlation matrix. For higher dimension covariance
matrices it is suggested you use the inverse wishart distribution
(transformed to a correlation matrix) for the simulations.
The covariance/variance prior is simulated from RxODEs cvPost()
function.
11.3.1 Simulation from inverse Wishart correlations
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)
s <-sim %>% confint(c("centr", "eff"))#> summarizing data...doneplot(s)
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.5728809
#>
#> [[2]]
#> [,1]
#> [1,] 0.3465021
#>
#> [[3]]
#> [,1]
#> [1,] 0.1386869
#>
#> [[4]]
#> [,1]
#> [1,] 0.1570577
#>
#> [[5]]
#> [,1]
#> [1,] 0.1677731
#>
#> [[6]]
#> [,1]
#> [1,] 0.3184372head(sim$sigma.list)#> [[1]]
#> [,1] [,2]
#> [1,] 0.093539238 0.007270049
#> [2,] 0.007270049 0.098648424
#>
#> [[2]]
#> [,1] [,2]
#> [1,] 0.109020277 -0.004127612
#> [2,] -0.004127612 0.087054268
#>
#> [[3]]
#> [,1] [,2]
#> [1,] 0.10606530 0.01457913
#> [2,] 0.01457913 0.10189653
#>
#> [[4]]
#> [,1] [,2]
#> [1,] 0.1025867133 -0.0007429996
#> [2,] -0.0007429996 0.0962922149
#>
#> [[5]]
#> [,1] [,2]
#> [1,] 0.098080929 -0.006730568
#> [2,] -0.006730568 0.112366768
#>
#> [[6]]
#> [,1] [,2]
#> [1,] 0.1123437 0.0188019
#> [2,] 0.0188019 0.1021367You can also see the parameter realizations from the $params data frame.
11.3.2 Simulate using variance/standard deviation standard errors
Lets assume we wish to simulate from the nonmem run included in xpose
First we setup the model:
rx1 <- RxODE({
cl <- tcl*(1+crcl.cl*(CLCR-65)) * exp(eta.v)
v <- tv * WT * exp(eta.v)
ka <- tka * exp(eta.ka)
ipred <- linCmt()
obs <- ipred * (1 + prop.sd) + add.sd
})Next we input the estimated parameters:
theta <- c(tcl=2.63E+01, tv=1.35E+00, tka=4.20E+00, tlag=2.08E-01,
prop.sd=2.05E-01, add.sd=1.06E-02, crcl.cl=7.17E-03,
## Note that since we are using the separation strategy the ETA variances are here too
eta.cl=7.30E-02, eta.v=3.80E-02, eta.ka=1.91E+00)And also their covariances; To me, the easiest way to create a named
covariance matrix is to use lotri():
thetaMat <- lotri(
tcl + tv + tka + tlag + prop.sd + add.sd + crcl.cl + eta.cl + eta.v + eta.ka ~
c(7.95E-01,
2.05E-02, 1.92E-03,
7.22E-02, -8.30E-03, 6.55E-01,
-3.45E-03, -6.42E-05, 3.22E-03, 2.47E-04,
8.71E-04, 2.53E-04, -4.71E-03, -5.79E-05, 5.04E-04,
6.30E-04, -3.17E-06, -6.52E-04, -1.53E-05, -3.14E-05, 1.34E-05,
-3.30E-04, 5.46E-06, -3.15E-04, 2.46E-06, 3.15E-06, -1.58E-06, 2.88E-06,
-1.29E-03, -7.97E-05, 1.68E-03, -2.75E-05, -8.26E-05, 1.13E-05, -1.66E-06, 1.58E-04,
-1.23E-03, -1.27E-05, -1.33E-03, -1.47E-05, -1.03E-04, 1.02E-05, 1.67E-06, 6.68E-05, 1.56E-04,
7.69E-02, -7.23E-03, 3.74E-01, 1.79E-03, -2.85E-03, 1.18E-05, -2.54E-04, 1.61E-03, -9.03E-04, 3.12E-01))
evw <- et(amount.units="mg", time.units="hours") %>%
et(amt=100) %>%
## For this problem we will simulate with sampling windows
et(list(c(0, 0.5),
c(0.5, 1),
c(1, 3),
c(3, 6),
c(6, 12))) %>%
et(id=1:1000)
### From the run we know that:
### total number of observations is: 476
### Total number of individuals: 74
sim <- rxSolve(rx1, theta, evw, nSub=100, nStud=10,
thetaMat=thetaMat,
## Match boundaries of problem
thetaLower=0,
sigma=c("prop.sd", "add.sd"), ## Sigmas are standard deviations
sigmaXform="identity", # default sigma xform="identity"
omega=c("eta.cl", "eta.v", "eta.ka"), ## etas are variances
omegaXform="variance", # default omega xform="variance"
iCov=data.frame(WT=rnorm(1000, 70, 15), CLCR=rnorm(1000, 65, 25)),
dfSub=74, dfObs=476);
print(sim)#> ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ Solved RxODE object ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂
#> ── Parameters ($params): ───────────────────────────────────
#> # A tibble: 10,000 x 10
#> sim.id id tcl crcl.cl CLCR eta.v tv WT tka eta.ka
#> <int> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 1 27.0 1.04 54.0 0.907 2.00 71.8 5.69 -0.153
#> 2 1 2 27.0 1.04 19.7 -0.225 2.00 80.2 5.69 0.249
#> 3 1 3 27.0 1.04 45.7 1.66 2.00 66.3 5.69 0.236
#> 4 1 4 27.0 1.04 73.9 0.556 2.00 69.4 5.69 -0.156
#> 5 1 5 27.0 1.04 91.4 0.296 2.00 45.5 5.69 -0.331
#> 6 1 6 27.0 1.04 94.9 -0.680 2.00 35.8 5.69 0.372
#> 7 1 7 27.0 1.04 13.6 -0.327 2.00 95.9 5.69 -0.0760
#> 8 1 8 27.0 1.04 66.2 0.589 2.00 57.3 5.69 0.688
#> 9 1 9 27.0 1.04 71.7 -0.611 2.00 41.0 5.69 0.212
#> 10 1 10 27.0 1.04 76.6 0.00250 2.00 66.5 5.69 0.243
#> # … with 9,990 more rows
#> ── Initial Conditions ($inits): ────────────────────────────
#> named numeric(0)
#>
#> Simulation with uncertainty in:
#> ● parameters (sim$thetaMat for changes)
#> ● omega matrix (sim$omegaList)
#> ● sigma matrix (sim$sigmaList)
#>
#> ── First part of data (object): ────────────────────────────
#> # A tibble: 50,000 x 8
#> sim.id id time cl v ka ipred obs
#> <int> <int> [h] <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 1 0.20072222 -696. 356. 4.88 NA NA
#> 2 1 1 0.79938985 -696. 356. 4.88 NA NA
#> 3 1 1 2.50526151 -696. 356. 4.88 NA NA
#> 4 1 1 3.38595486 -696. 356. 4.88 NA NA
#> 5 1 1 9.28579107 -696. 356. 4.88 NA NA
#> 6 1 2 0.04197341 -992. 128. 7.30 NA NA
#> # … with 49,994 more rows
#> ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂### Notice that the simulation time-points change for the individual
### If you want the same sampling time-points you can do that as well:
evw <- et(amount.units="mg", time.units="hours") %>%
et(amt=100) %>%
et(0, 24, length.out=50) %>%
et(id=1:100)
sim <- rxSolve(rx1, theta, evw, nSub=100, nStud=10,
thetaMat=thetaMat,
## Match boundaries of problem
thetaLower=0,
sigma=c("prop.sd", "add.sd"), ## Sigmas are standard deviations
sigmaXform="identity", # default sigma xform="identity"
omega=c("eta.cl", "eta.v", "eta.ka"), ## etas are variances
omegaXform="variance", # default omega xform="variance"
iCov=data.frame(WT=rnorm(1000, 70, 15), CLCR=rnorm(1000, 65, 25)),
dfSub=74, dfObs=476)
s <-sim %>% confint(c("ipred"))#> summarizing data...doneplot(s)
11.3.3 Simulate without uncertainty in omega or sigma parameters
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:
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 <- lotri(eff.err ~ 0.1, cp.err ~ 0.1)
sim <- rxSolve(mod, theta, ev, omega=omega, nSub=100, sigma=sigma,
thetaMat=tMat, nStud=10,
simVariability=FALSE)
s <-sim %>% confint(c("centr", "eff"))#> summarizing data...doneplot(s)
Note since realizations of omega and sigma were not simulated,
$omega.list and $sigma.list both return NULL.
11.3.3.0.1 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
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 (either mvnfast::rmvn() or
mvnfast::rmvt()). 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);
s <-sim %>% confint(c("eff", "centr"))#> summarizing data...doneThe default for this is 1 core since the result depends on the
number of cores and the random seed you use in your simulation as well
as the work-load each thread is sharing/architecture. 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.
11.4 Using prior data for solving
RxODE can use a single subject or
multiple subjects with a single event table to
solve ODEs. Additionally, RxODE can use an arbitrary data frame with
individualized events. For example when using nlmixr, you could use
the RxODE/vignettes/theo_sd data frame
library(RxODE)
### Load data from nlmixr
d <- qs::qread("RxODE/vignettes/theo_sd.qs")
### Create RxODE model
theo <- RxODE({
tka ~ 0.45 # Log Ka
tcl ~ 1 # Log Cl
tv ~ 3.45 # Log V
eta.ka ~ 0.6
eta.cl ~ 0.3
eta.v ~ 0.1
ka <- exp(tka + eta.ka)
cl <- exp(tcl + eta.cl)
v <- exp(tv + eta.v)
d/dt(depot) = -ka * depot
d/dt(center) = ka * depot - cl / v * center
cp = center / v
})
### Create parameter dataset
library(dplyr)
parsDf <- tribble(
~ eta.ka, ~ eta.cl, ~ eta.v,
0.105, -0.487, -0.080,
0.221, 0.144, 0.021,
0.368, 0.031, 0.058,
-0.277, -0.015, -0.007,
-0.046, -0.155, -0.142,
-0.382, 0.367, 0.203,
-0.791, 0.160, 0.047,
-0.181, 0.168, 0.096,
1.420, 0.042, 0.012,
-0.738, -0.391, -0.170,
0.790, 0.281, 0.146,
-0.527, -0.126, -0.198) %>%
mutate(tka = 0.451, tcl = 1.017, tv = 3.449)
### Now solve the dataset
solveData <- rxSolve(theo, parsDf, d)
plot(solveData, cp)
print(solveData)#> ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ Solved RxODE object ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂
#> ── Parameters ($params): ───────────────────────────────────
#> # A tibble: 12 x 1
#> id
#> <fct>
#> 1 1
#> 2 2
#> 3 3
#> 4 4
#> 5 5
#> 6 6
#> 7 7
#> 8 8
#> 9 9
#> 10 10
#> 11 11
#> 12 12
#> ── Initial Conditions ($inits): ────────────────────────────
#> depot center
#> 0 0
#> ── First part of data (object): ────────────────────────────
#> # A tibble: 132 x 8
#> id time ka cl v cp depot center
#> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0 2.86 3.67 34.8 0 320. 0
#> 2 1 0.25 2.86 3.67 34.8 4.62 157. 161.
#> 3 1 0.570 2.86 3.67 34.8 7.12 62.8 248.
#> 4 1 1.12 2.86 3.67 34.8 8.09 13.0 282.
#> 5 1 2.02 2.86 3.67 34.8 7.68 0.996 267.
#> 6 1 3.82 2.86 3.67 34.8 6.38 0.00581 222.
#> # … with 126 more rows
#> ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂### Of course the fasest way to solve if you don't care about the RxODE extra parameters is
solveData <- rxSolve(theo, parsDf, d, returnType="data.frame")
### solved data
dplyr::as.tbl(solveData)#> # A tibble: 132 x 8
#> id time ka cl v cp depot center
#> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0 2.86 3.67 34.8 0 3.20e+2 0
#> 2 1 0.25 2.86 3.67 34.8 4.62 1.57e+2 161.
#> 3 1 0.570 2.86 3.67 34.8 7.12 6.28e+1 248.
#> 4 1 1.12 2.86 3.67 34.8 8.09 1.30e+1 282.
#> 5 1 2.02 2.86 3.67 34.8 7.68 9.96e-1 267.
#> 6 1 3.82 2.86 3.67 34.8 6.38 5.81e-3 222.
#> 7 1 5.1 2.86 3.67 34.8 5.58 1.50e-4 194.
#> 8 1 7.03 2.86 3.67 34.8 4.55 6.02e-7 158.
#> 9 1 9.05 2.86 3.67 34.8 3.68 1.77e-9 128.
#> 10 1 12.1 2.86 3.67 34.8 2.66 9.43e-9 92.6
#> # … with 122 more rowsdata.table::data.table(solveData)#> id time ka cl v cp depot center
#> 1: 1 0.00 2.857651 3.669297 34.81332 0.0000000 3.199920e+02 0.00000
#> 2: 1 0.25 2.857651 3.669297 34.81332 4.6240421 1.566295e+02 160.97825
#> 3: 1 0.57 2.857651 3.669297 34.81332 7.1151647 6.276731e+01 247.70249
#> 4: 1 1.12 2.857651 3.669297 34.81332 8.0922106 1.303613e+01 281.71670
#> 5: 1 2.02 2.857651 3.669297 34.81332 7.6837844 9.958446e-01 267.49803
#> ---
#> 128: 12 5.07 2.857651 3.669297 34.81332 5.6044213 1.636210e-04 195.10850
#> 129: 12 7.07 2.857651 3.669297 34.81332 4.5392337 5.385697e-07 158.02579
#> 130: 12 9.03 2.857651 3.669297 34.81332 3.6920276 1.882087e-09 128.53173
#> 131: 12 12.05 2.857651 3.669297 34.81332 2.6855080 8.461424e-09 93.49144
#> 132: 12 24.15 2.857651 3.669297 34.81332 0.7501667 -4.775222e-10 26.11579