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:
#> RxODE 1.1.0 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
})
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.0729555976 0.029714842 -0.009723169 -0.064199455 0.01069553
#> TCl 0.0297148418 0.167692128 0.037094746 0.009623620 0.08034273
#> V2 -0.0097231694 0.037094746 0.070742277 0.087151404 0.07359617
#> Q -0.0641994547 0.009623620 0.087151404 0.200429705 0.07271673
#> V3 0.0106955306 0.080342726 0.073596173 0.072716732 0.14556855
#> Kin 0.0070440826 0.034146017 0.021796701 -0.060232686 0.06325513
#> Kout 0.0008130782 -0.007198921 -0.004865182 0.006195073 -0.02146333
#> EC50 0.0158844256 -0.087004181 0.009428023 -0.034654007 -0.01422844
#> Kin Kout EC50
#> TKA 0.007044083 0.0008130782 0.015884426
#> TCl 0.034146017 -0.0071989210 -0.087004181
#> V2 0.021796701 -0.0048651824 0.009428023
#> Q -0.060232686 0.0061950729 -0.034654007
#> V3 0.063255125 -0.0214633256 -0.014228439
#> Kin 0.189831227 -0.0698227583 0.081319255
#> Kout -0.069822758 0.0472685846 -0.039719447
#> EC50 0.081319255 -0.0397194466 0.174344448To 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)#> 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.420 -0.213 0.120
#> 2 1 2 -0.420 -0.213 0.120
#> 3 1 3 -0.420 -0.213 0.120
#> 4 1 4 -0.420 -0.213 0.120
#> 5 1 5 -0.420 -0.213 0.120
#> 6 1 6 -0.420 -0.213 0.120
#> 7 1 7 -0.420 -0.213 0.120
#> 8 1 8 -0.420 -0.213 0.120
#> 9 1 9 -0.420 -0.213 0.120
#> 10 1 10 -0.420 -0.213 0.120
#> # ... 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 C3
#> <int> <int> [h] <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 1 0 -0.420 0.120 0.391 -0.0835 -0.197 0.0166 0 0
#> 2 1 1 0.1 -0.420 0.120 -0.491 -0.0844 -0.197 0.0166 4.48 0.00651
#> 3 1 1 4 -0.420 0.120 0.391 -0.0835 -0.197 0.0166 62.8 5.18
#> 4 1 1 4.1 -0.420 0.120 -0.491 -0.0844 -0.197 0.0166 84.9 5.35
#> 5 1 1 8 -0.420 0.120 0.391 -0.0835 -0.197 0.0166 -11.1 11.1
#> 6 1 1 8.1 -0.420 0.120 -0.491 -0.0844 -0.197 0.0166 59.7 11.2
#> # ... with 975,994 more rows, and 7 more variables: 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.4203958 -0.2130085 0.1204075 -0.006785869
#> 2 1 2 -0.4203958 -0.2130085 0.1204075 -0.006785869
#> 3 1 3 -0.4203958 -0.2130085 0.1204075 -0.006785869
#> 4 1 4 -0.4203958 -0.2130085 0.1204075 -0.006785869
#> 5 1 5 -0.4203958 -0.2130085 0.1204075 -0.006785869
#> 6 1 6 -0.4203958 -0.2130085 0.1204075 -0.006785869
#> eye.Cl(eye==1) eye.Cl(eye==2) eye.Ka(eye==1) eye.Ka(eye==2) iov.Cl(occ==1)
#> 1 0.3912619 -0.49111561 -0.083490021 -0.08442088 -0.19740065
#> 2 0.2959090 0.26700321 -0.159890637 -0.29441094 0.04844301
#> 3 -0.3777342 0.50029603 0.072327594 0.10659733 -0.01474296
#> 4 -0.1336540 -0.25191280 0.320628003 -0.01745451 0.04597216
#> 5 0.1686718 -0.32068438 0.009463607 0.11629676 -0.24519393
#> 6 -0.2957091 -0.06857382 -0.444767143 -0.13273211 -0.14083012
#> iov.Cl(occ==2) iov.Ka(occ==1) iov.Ka(occ==2) V2 V3 TCl
#> 1 -0.077611118 0.01657936 0.04427318 40.33104 297.0264 18.49128
#> 2 -0.130268048 -0.08077711 -0.06269579 40.33104 297.0264 18.49128
#> 3 -0.089252982 0.13271402 -0.16917194 40.33104 297.0264 18.49128
#> 4 -0.038870212 -0.03054084 0.18321811 40.33104 297.0264 18.49128
#> 5 -0.003213317 -0.04190360 -0.03084874 40.33104 297.0264 18.49128
#> 6 0.086041699 0.13654626 -0.04577458 40.33104 297.0264 18.49128
#> eta.Cl TKA eta.Ka Q Kin Kout EC50
#> 1 -0.15304459 0.122184 0.32696580 11.09409 0.1773626 1.445233 199.5471
#> 2 0.30366032 0.122184 -0.14854962 11.09409 0.1773626 1.445233 199.5471
#> 3 0.08392156 0.122184 0.05032808 11.09409 0.1773626 1.445233 199.5471
#> 4 -0.12655485 0.122184 0.39933195 11.09409 0.1773626 1.445233 199.5471
#> 5 -0.05820359 0.122184 -0.36646997 11.09409 0.1773626 1.445233 199.5471
#> 6 0.06330689 0.122184 -0.33706060 11.09409 0.1773626 1.445233 199.5471#> sim.id id inv.Cl(inv==1) inv.Cl(inv==2) inv.Ka(inv==1) inv.Ka(inv==2)
#> 1 2 1 0.1145848 0.09457461 -0.05362635 0.07696309
#> 2 2 2 0.1145848 0.09457461 -0.05362635 0.07696309
#> 3 2 3 0.1145848 0.09457461 -0.05362635 0.07696309
#> 4 2 4 0.1145848 0.09457461 -0.05362635 0.07696309
#> 5 2 5 0.1145848 0.09457461 -0.05362635 0.07696309
#> 6 2 6 0.1145848 0.09457461 -0.05362635 0.07696309
#> eye.Cl(eye==1) eye.Cl(eye==2) eye.Ka(eye==1) eye.Ka(eye==2) iov.Cl(occ==1)
#> 1 -0.11866351 -0.18202287 -0.04725720 0.01369497 0.02384389
#> 2 -0.17858595 -0.10887058 -0.05337936 -0.17527610 0.02023782
#> 3 -0.11622757 -0.10336398 0.15710881 0.12572783 -0.17043609
#> 4 0.11627163 -0.33468907 0.08544801 0.17281907 -0.21272115
#> 5 -0.16203306 -0.05647415 0.24757833 0.26241044 -0.02580006
#> 6 0.07297872 -0.05111628 -0.12391057 0.06350245 0.09464866
#> iov.Cl(occ==2) iov.Ka(occ==1) iov.Ka(occ==2) V2 V3 TCl
#> 1 0.0080529792 0.031210414 -0.103996745 40.36877 297.4956 18.63833
#> 2 -0.1149723077 0.041828770 0.145628662 40.36877 297.4956 18.63833
#> 3 0.2085954354 -0.139951130 -0.176427055 40.36877 297.4956 18.63833
#> 4 0.0300091296 -0.132232256 0.013298749 40.36877 297.4956 18.63833
#> 5 -0.0001739385 0.086994738 0.123551590 40.36877 297.4956 18.63833
#> 6 0.1097545737 0.006786617 -0.007701577 40.36877 297.4956 18.63833
#> eta.Cl TKA eta.Ka Q Kin Kout EC50
#> 1 -0.19738241 0.009073191 -0.025243233 10.8825 1.521469 0.6886986 199.9741
#> 2 0.49397764 0.009073191 0.113357588 10.8825 1.521469 0.6886986 199.9741
#> 3 -0.03006838 0.009073191 -0.161780483 10.8825 1.521469 0.6886986 199.9741
#> 4 -0.28412610 0.009073191 0.480492561 10.8825 1.521469 0.6886986 199.9741
#> 5 0.36163390 0.009073191 0.008356588 10.8825 1.521469 0.6886986 199.9741
#> 6 -0.36472429 0.009073191 0.001865549 10.8825 1.521469 0.6886986 199.9741For 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.