More than one level of nesting is possible in RxODE; In this example we will be using the following uncertainties and sources of variability:
| Level | Variable | Matrix specified | Integrated Matrix |
|---|---|---|---|
| Model uncertainty | NA | thetaMat |
thetaMat |
| Investigator |
inv.Cl, inv.Ka
|
omega |
theta |
| Subject |
eta.Cl, eta.Ka
|
omega |
omega |
| Eye |
eye.Cl, eye.Ka
|
omega |
omega |
| Occasion |
iov.Cl, occ.Ka
|
omega |
omega |
| Unexplained Concentration | prop.sd |
sigma |
sigma |
| Unexplained Effect | add.sd |
sigma |
sigma |
This event table contains nesting variables: - inv: investigator id - id: subject id - eye: eye id (left or right) - occ: occasion
#> RxODE 1.0.2 using 1 threads (see ?getRxThreads)
#> no cache: create with `rxCreateCache()`#>
#> Attaching package: 'dplyr'#> The following objects are masked from 'package:stats':
#>
#> filter, lag#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
et(amountUnits="mg", timeUnits="hours") %>%
et(amt=10000, addl=9,ii=12,cmt="depot") %>%
et(time=120, amt=2000, addl=4, ii=14, cmt="depot") %>%
et(seq(0, 240, by=4)) %>% # Assumes sampling when there is no dosing information
et(seq(0, 240, by=4) + 0.1) %>% ## adds 0.1 for separate eye
et(id=1:20) %>%
## Add an occasion per dose
mutate(occ=cumsum(!is.na(amt))) %>%
mutate(occ=ifelse(occ == 0, 1, occ)) %>%
mutate(occ=2- occ %% 2) %>%
mutate(eye=ifelse(round(time) == time, 1, 2)) %>%
mutate(inv=ifelse(id < 10, 1, 2)) %>% as_tibble ->
evThis creates the RxODE model with multi-level nesting. Note the variables inv.Cl, inv.Ka, eta.Cl etc; You only need one variable for each level of nesting.
mod <- RxODE({
## Clearance with individuals
eff(0) = 1
C2 = centr/V2*(1+prop.sd);
C3 = peri/V3;
CL = TCl*exp(eta.Cl + eye.Cl + iov.Cl + inv.Cl)
KA = TKA * exp(eta.Ka + eye.Ka + iov.Cl + inv.Ka)
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;
ef0 = eff + add.sd
})#> qs v0.23.5.
theta <- c("TKA"=0.294, "TCl"=18.6, "V2"=40.2,
"Q"=10.5, "V3"=297, "Kin"=1, "Kout"=1, "EC50"=200)
## Creating covariance matrix
tmp <- matrix(rnorm(8^2), 8, 8)
tMat <- tcrossprod(tmp, tmp) / (8 ^ 2)
dimnames(tMat) <- list(names(theta), names(theta))
tMat#> TKA TCl V2 Q V3
#> TKA 0.107924293 -0.007655364 -0.002449253 -0.085719657 -0.04775912
#> TCl -0.007655364 0.198122861 0.021732655 0.007892284 -0.04167297
#> V2 -0.002449253 0.021732655 0.139499928 -0.003829790 0.05130865
#> Q -0.085719657 0.007892284 -0.003829790 0.159093401 0.01927711
#> V3 -0.047759123 -0.041672967 0.051308645 0.019277114 0.09019371
#> Kin 0.008093638 0.046149166 0.040853353 -0.036598728 0.03697735
#> Kout -0.023869002 0.100602141 0.005184793 0.033359694 -0.02314619
#> EC50 0.039314747 -0.008432118 -0.010809544 0.003862315 -0.01862489
#> Kin Kout EC50
#> TKA 0.008093638 -0.023869002 0.039314747
#> TCl 0.046149166 0.100602141 -0.008432118
#> V2 0.040853353 0.005184793 -0.010809544
#> Q -0.036598728 0.033359694 0.003862315
#> V3 0.036977346 -0.023146186 -0.018624885
#> Kin 0.139912288 0.014241132 -0.014681951
#> Kout 0.014241132 0.089777221 0.003002065
#> EC50 -0.014681951 0.003002065 0.035658813To specify multiple levels of nesting, you can specify it as a nested lotri matrix; When using this approach you use the condition operator | to specify what variable nesting occurs on; For the Bayesian simulation we need to specify how much information we have for each parameter; For RxODE this is the nu parameter.
In this case: - id, nu=100 or the model came from 100 subjects - eye, nu=200 or the model came from 200 eyes - occ, nu=200 or the model came from 200 occasions - inv, nu=10 or the model came from 10 investigators
To specify this in lotri you can use | var(nu=X), or:
omega <- lotri(lotri(eta.Cl ~ 0.1,
eta.Ka ~ 0.1) | id(nu=100),
lotri(eye.Cl ~ 0.05,
eye.Ka ~ 0.05) | eye(nu=200),
lotri(iov.Cl ~ 0.01,
iov.Ka ~ 0.01) | occ(nu=200),
lotri(inv.Cl ~ 0.02,
inv.Ka ~ 0.02) | inv(nu=10))
omega#> $id
#> eta.Cl eta.Ka
#> eta.Cl 0.1 0.0
#> eta.Ka 0.0 0.1
#>
#> $eye
#> eye.Cl eye.Ka
#> eye.Cl 0.05 0.00
#> eye.Ka 0.00 0.05
#>
#> $occ
#> iov.Cl iov.Ka
#> iov.Cl 0.01 0.00
#> iov.Ka 0.00 0.01
#>
#> $inv
#> inv.Cl inv.Ka
#> inv.Cl 0.02 0.00
#> inv.Ka 0.00 0.02
#>
#> Properties: nuThe last piece of variability to specify is the unexplained variability
sigma <- lotri(prop.sd ~ .25,
add.sd~ 0.125)
s <- rxSolve(mod, theta, ev,
thetaMat=tMat, omega=omega,
sigma=sigma, sigmaDf=400,
nStud=400)#> intdy -- t = 72.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 76 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 76.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 80 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 80.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 84 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 84.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 88 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 88.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 92 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 92.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 96 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 96.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 100 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 100.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 104 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 104.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 108 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 108.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 112 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 112.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 116 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 116.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 120 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 120.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 124 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 124.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 128 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 128.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 132 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 132.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 134 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 136 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 136.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 140 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 140.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 144 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 144.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 148 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 148.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 152 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 152.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 156 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 156.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 160 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 160.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 162 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 164 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 164.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 168 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 168.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 172 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 172.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 176 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 176.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 180 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 180.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 184 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 184.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 188 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 188.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 192 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 192.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 196 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 196.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 200 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 200.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 204 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 204.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 208 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 208.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 212 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 212.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 216 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 216.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 220 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 220.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 224 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 224.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 228 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 228.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 232 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 232.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 236 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 236.1 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 240 illegal. t not in interval tcur - _C(hu) to tcur
#> intdy -- t = 240.1 illegal. t not in interval tcur - _C(hu) to tcur
#> unhandled error message: EE:[lsoda] trouble from intdy, itask = 1, tout = 240.1
#> @(lsoda.c:689
#> unhandled error message: EE:[lsoda] 70000 steps taken before reaching tout
#> @(lsoda.c:748#> Warning: some ID(s) could not solve the ODEs correctly; These values are
#> replaced with 'NA'
print(s)#> ______________________________ Solved RxODE object _____________________________
#> -- Parameters ($params): -------------------------------------------------------
#> # A tibble: 8,000 x 24
#> sim.id id `inv.Cl(inv==1)` `inv.Cl(inv==2)` `inv.Ka(inv==1)`
#> <int> <fct> <dbl> <dbl> <dbl>
#> 1 1 1 0.276 0.0471 -0.0740
#> 2 1 2 0.276 0.0471 -0.0740
#> 3 1 3 0.276 0.0471 -0.0740
#> 4 1 4 0.276 0.0471 -0.0740
#> 5 1 5 0.276 0.0471 -0.0740
#> 6 1 6 0.276 0.0471 -0.0740
#> 7 1 7 0.276 0.0471 -0.0740
#> 8 1 8 0.276 0.0471 -0.0740
#> 9 1 9 0.276 0.0471 -0.0740
#> 10 1 10 0.276 0.0471 -0.0740
#> # ... with 7,990 more rows, and 19 more variables: `inv.Ka(inv==2)` <dbl>,
#> # `eye.Cl(eye==1)` <dbl>, `eye.Cl(eye==2)` <dbl>, `eye.Ka(eye==1)` <dbl>,
#> # `eye.Ka(eye==2)` <dbl>, `iov.Cl(occ==1)` <dbl>, `iov.Cl(occ==2)` <dbl>,
#> # `iov.Ka(occ==1)` <dbl>, `iov.Ka(occ==2)` <dbl>, V2 <dbl>, V3 <dbl>,
#> # TCl <dbl>, eta.Cl <dbl>, TKA <dbl>, eta.Ka <dbl>, Q <dbl>, Kin <dbl>,
#> # Kout <dbl>, EC50 <dbl>
#> -- Initial Conditions ($inits): ------------------------------------------------
#> depot centr peri eff
#> 0 0 0 1
#>
#> Simulation with uncertainty in:
#> * parameters (s$thetaMat for changes)
#> * omega matrix (s$omegaList)
#>
#> -- First part of data (object): ------------------------------------------------
#> # A tibble: 976,000 x 18
#> sim.id id time inv.Cl inv.Ka eye.Cl eye.Ka iov.Cl iov.Ka C2
#> <int> <int> [h] <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 1 0.0 0.276 -0.0740 -0.0638 -0.179 0.0785 -0.0910 0
#> 2 1 1 0.1 0.276 -0.0740 0.151 -0.152 0.0785 -0.0910 -6.53
#> 3 1 1 4.0 0.276 -0.0740 -0.0638 -0.179 0.0785 -0.0910 -113.
#> 4 1 1 4.1 0.276 -0.0740 0.151 -0.152 0.0785 -0.0910 -229.
#> 5 1 1 8.0 0.276 -0.0740 -0.0638 -0.179 0.0785 -0.0910 -329.
#> 6 1 1 8.1 0.276 -0.0740 0.151 -0.152 0.0785 -0.0910 17.9
#> # ... with 975,994 more rows, and 8 more variables: C3 <dbl>, CL <dbl>,
#> # KA <dbl>, ef0 <dbl>, depot <dbl>, centr <dbl>, peri <dbl>, eff <dbl>
#> ________________________________________________________________________________There are multiple investigators in a study; Each investigator has a number of individuals enrolled at their site. RxODE automatically determines the number of investigators and then will simulate an effect for each investigator. With the output, inv.Cl(inv==1) will be the inv.Cl for investigator 1, inv.Cl(inv==2) will be the inv.Cl for investigator 2, etc.
inv.Cl(inv==1), inv.Cl(inv==2), etc will be simulated for each study and then combined to form the between investigator variability. In equation form these represent the following:
inv.Cl = (inv == 1) * `inv.Cl(inv==1)` + (inv == 2) * `inv.Cl(inv==2)`If you look at the simulated parameters you can see inv.Cl(inv==1) and inv.Cl(inv==2) are in the s$params; They are the same for each study:
#> sim.id id inv.Cl(inv==1) inv.Cl(inv==2) inv.Ka(inv==1) inv.Ka(inv==2)
#> 1 1 1 0.2755711 0.04714543 -0.0740288 -0.008040342
#> 2 1 2 0.2755711 0.04714543 -0.0740288 -0.008040342
#> 3 1 3 0.2755711 0.04714543 -0.0740288 -0.008040342
#> 4 1 4 0.2755711 0.04714543 -0.0740288 -0.008040342
#> 5 1 5 0.2755711 0.04714543 -0.0740288 -0.008040342
#> 6 1 6 0.2755711 0.04714543 -0.0740288 -0.008040342
#> eye.Cl(eye==1) eye.Cl(eye==2) eye.Ka(eye==1) eye.Ka(eye==2) iov.Cl(occ==1)
#> 1 -0.06377666 0.15105685 -0.1786692 -0.1520267 0.07846531
#> 2 0.08693985 -0.15573672 0.1557618 -0.1099414 0.10478330
#> 3 0.21274973 -0.15195973 -0.2579102 0.1296449 0.11643389
#> 4 0.23157918 0.07613762 -0.1328889 -0.1628823 0.07722827
#> 5 -0.24177606 -0.21347582 -0.0706487 0.1202125 -0.14104513
#> 6 -0.20997741 -0.08491303 -0.1173625 0.1611406 0.02696981
#> iov.Cl(occ==2) iov.Ka(occ==1) iov.Ka(occ==2) V2 V3 TCl
#> 1 0.08962458 -0.09101293 -0.058152279 40.5109 297.5977 18.85686
#> 2 -0.05666680 0.13200667 -0.091116589 40.5109 297.5977 18.85686
#> 3 0.11274088 0.06402374 0.005631993 40.5109 297.5977 18.85686
#> 4 -0.08106882 0.15253238 0.050789569 40.5109 297.5977 18.85686
#> 5 -0.04479171 -0.09903988 -0.072261824 40.5109 297.5977 18.85686
#> 6 0.09823473 0.15181825 -0.042597511 40.5109 297.5977 18.85686
#> eta.Cl TKA eta.Ka Q Kin Kout EC50
#> 1 0.15896813 -0.3099298 0.092231031 10.52823 1.106277 1.032151 199.6466
#> 2 -0.26522434 -0.3099298 0.251118936 10.52823 1.106277 1.032151 199.6466
#> 3 -0.39932187 -0.3099298 0.184285841 10.52823 1.106277 1.032151 199.6466
#> 4 0.01730799 -0.3099298 -0.012962317 10.52823 1.106277 1.032151 199.6466
#> 5 0.09462363 -0.3099298 -0.053433989 10.52823 1.106277 1.032151 199.6466
#> 6 -0.10657556 -0.3099298 0.005568882 10.52823 1.106277 1.032151 199.6466#> sim.id id inv.Cl(inv==1) inv.Cl(inv==2) inv.Ka(inv==1) inv.Ka(inv==2)
#> 1 2 1 0.03731579 -0.2119719 0.1565425 0.1257644
#> 2 2 2 0.03731579 -0.2119719 0.1565425 0.1257644
#> 3 2 3 0.03731579 -0.2119719 0.1565425 0.1257644
#> 4 2 4 0.03731579 -0.2119719 0.1565425 0.1257644
#> 5 2 5 0.03731579 -0.2119719 0.1565425 0.1257644
#> 6 2 6 0.03731579 -0.2119719 0.1565425 0.1257644
#> eye.Cl(eye==1) eye.Cl(eye==2) eye.Ka(eye==1) eye.Ka(eye==2) iov.Cl(occ==1)
#> 1 0.03603852 -0.02281387 -0.11280622 -0.09882465 -0.132766001
#> 2 0.12439895 0.11196998 0.03405942 -0.18839087 0.043255148
#> 3 0.08829691 -0.04782738 0.19514921 0.05632306 0.005924353
#> 4 -0.20881694 0.16064290 0.57907698 0.04341776 0.252851591
#> 5 0.25403448 -0.32433222 -0.13117535 0.14646807 -0.054967817
#> 6 0.07558715 0.04246528 -0.08851438 0.30159674 -0.044378421
#> iov.Cl(occ==2) iov.Ka(occ==1) iov.Ka(occ==2) V2 V3 TCl
#> 1 -0.067357439 -0.003930598 0.13455989 40.41579 296.6274 19.85728
#> 2 -0.019515686 -0.166474175 -0.10432113 40.41579 296.6274 19.85728
#> 3 -0.071215057 -0.044724440 -0.01847172 40.41579 296.6274 19.85728
#> 4 -0.005191752 0.102101931 0.03393396 40.41579 296.6274 19.85728
#> 5 -0.069187495 -0.030763626 -0.03931005 40.41579 296.6274 19.85728
#> 6 0.094235252 0.009555823 -0.05426202 40.41579 296.6274 19.85728
#> eta.Cl TKA eta.Ka Q Kin Kout EC50
#> 1 0.02154255 0.4509381 0.20742231 10.06732 1.063482 1.573717 199.9159
#> 2 -0.22965387 0.4509381 0.53253312 10.06732 1.063482 1.573717 199.9159
#> 3 0.21979206 0.4509381 0.73271925 10.06732 1.063482 1.573717 199.9159
#> 4 -0.43465238 0.4509381 0.14643090 10.06732 1.063482 1.573717 199.9159
#> 5 -0.01564301 0.4509381 -0.06340707 10.06732 1.063482 1.573717 199.9159
#> 6 0.04505743 0.4509381 0.58289795 10.06732 1.063482 1.573717 199.9159For between eye variability and between occasion variability each individual simulates a number of variables that become the between eye and between occasion variability; In the case of the eye:
eye.Cl = (eye == 1) * `eye.Cl(eye==1)` + (eye == 2) * `eye.Cl(eye==2)`So when you look the simulation each of these variables (ie eye.Cl(eye==1), eye.Cl(eye==2), etc) they change for each individual and when combined make the between eye variability or the between occasion variability that can be seen in some pharamcometric models.