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.1 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.052299027 -0.0020034683 -0.013430812 0.0042390076 0.0413906822
#> TCl -0.002003468 0.0218331047 0.002766279 0.0182346413 -0.0332541792
#> V2 -0.013430812 0.0027662795 0.078878019 0.0428913911 -0.0064044043
#> Q 0.004239008 0.0182346413 0.042891391 0.0688584959 0.0005038676
#> V3 0.041390682 -0.0332541792 -0.006404404 0.0005038676 0.1168521908
#> Kin -0.028991201 0.0249455744 0.019150688 0.0346070136 -0.0700758592
#> Kout 0.003501976 -0.0009969396 -0.047376005 0.0044577667 0.0140943439
#> EC50 -0.012011849 -0.0314580393 -0.079212914 -0.0642195971 0.0287538400
#> Kin Kout EC50
#> TKA -0.028991201 0.0035019757 -0.01201185
#> TCl 0.024945574 -0.0009969396 -0.03145804
#> V2 0.019150688 -0.0473760047 -0.07921291
#> Q 0.034607014 0.0044577667 -0.06421960
#> V3 -0.070075859 0.0140943439 0.02875384
#> Kin 0.106305332 0.0091986293 -0.04788276
#> Kout 0.009198629 0.0845471861 0.06074777
#> EC50 -0.047882765 0.0607477694 0.16205720To 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:749#> 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.00275 -0.0420 -0.151
#> 2 1 2 -0.00275 -0.0420 -0.151
#> 3 1 3 -0.00275 -0.0420 -0.151
#> 4 1 4 -0.00275 -0.0420 -0.151
#> 5 1 5 -0.00275 -0.0420 -0.151
#> 6 1 6 -0.00275 -0.0420 -0.151
#> 7 1 7 -0.00275 -0.0420 -0.151
#> 8 1 8 -0.00275 -0.0420 -0.151
#> 9 1 9 -0.00275 -0.0420 -0.151
#> 10 1 10 -0.00275 -0.0420 -0.151
#> # ... 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.00275 -0.151 0.0266 0.313 0.134 -0.0798 0 0
#> 2 1 1 0.1 -0.00275 -0.151 0.457 0.0464 0.134 -0.0798 0.664 0.000843
#> 3 1 1 4 -0.00275 -0.151 0.0266 0.313 0.134 -0.0798 3.37 0.342
#> 4 1 1 4.1 -0.00275 -0.151 0.457 0.0464 0.134 -0.0798 4.58 0.352
#> 5 1 1 8 -0.00275 -0.151 0.0266 0.313 0.134 -0.0798 5.06 0.717
#> 6 1 1 8.1 -0.00275 -0.151 0.457 0.0464 0.134 -0.0798 4.41 0.725
#> # ... 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.002747036 -0.04199234 -0.1511898 0.09129996
#> 2 1 2 -0.002747036 -0.04199234 -0.1511898 0.09129996
#> 3 1 3 -0.002747036 -0.04199234 -0.1511898 0.09129996
#> 4 1 4 -0.002747036 -0.04199234 -0.1511898 0.09129996
#> 5 1 5 -0.002747036 -0.04199234 -0.1511898 0.09129996
#> 6 1 6 -0.002747036 -0.04199234 -0.1511898 0.09129996
#> eye.Cl(eye==1) eye.Cl(eye==2) eye.Ka(eye==1) eye.Ka(eye==2) iov.Cl(occ==1)
#> 1 0.02663260 0.45652885 0.31266166 0.04643933 0.133931213
#> 2 -0.00645066 -0.14663791 -0.05418058 0.46416699 -0.187409314
#> 3 -0.19104833 0.28809622 -0.21654429 0.23997922 0.028992966
#> 4 -0.43362901 0.18609466 -0.28564274 -0.25316655 0.043293203
#> 5 0.20127512 0.14003835 0.06926179 0.02192560 0.006925008
#> 6 0.44923971 0.06974016 -0.16320977 -0.11844843 0.108659686
#> iov.Cl(occ==2) iov.Ka(occ==1) iov.Ka(occ==2) V2 V3 TCl
#> 1 0.056467287 -0.079792617 -0.02431260 40.29231 296.7368 18.7175
#> 2 -0.027264704 -0.056838177 -0.11527060 40.29231 296.7368 18.7175
#> 3 0.001101861 -0.002305014 -0.02606005 40.29231 296.7368 18.7175
#> 4 -0.040806968 0.052022440 -0.08489036 40.29231 296.7368 18.7175
#> 5 -0.167999980 -0.075263205 0.04186347 40.29231 296.7368 18.7175
#> 6 0.046880336 0.009432103 -0.06332094 40.29231 296.7368 18.7175
#> eta.Cl TKA eta.Ka Q Kin Kout EC50
#> 1 0.05249875 0.02854144 -0.66296533 10.50871 1.19586 0.7642861 199.6414
#> 2 0.69046590 0.02854144 0.34608629 10.50871 1.19586 0.7642861 199.6414
#> 3 -0.16526149 0.02854144 0.07753312 10.50871 1.19586 0.7642861 199.6414
#> 4 0.17141880 0.02854144 0.07579197 10.50871 1.19586 0.7642861 199.6414
#> 5 0.04265318 0.02854144 0.09043223 10.50871 1.19586 0.7642861 199.6414
#> 6 0.04398132 0.02854144 0.41274350 10.50871 1.19586 0.7642861 199.6414#> sim.id id inv.Cl(inv==1) inv.Cl(inv==2) inv.Ka(inv==1) inv.Ka(inv==2)
#> 1 2 1 0.2142626 0.005337268 0.2257368 -0.06238657
#> 2 2 2 0.2142626 0.005337268 0.2257368 -0.06238657
#> 3 2 3 0.2142626 0.005337268 0.2257368 -0.06238657
#> 4 2 4 0.2142626 0.005337268 0.2257368 -0.06238657
#> 5 2 5 0.2142626 0.005337268 0.2257368 -0.06238657
#> 6 2 6 0.2142626 0.005337268 0.2257368 -0.06238657
#> eye.Cl(eye==1) eye.Cl(eye==2) eye.Ka(eye==1) eye.Ka(eye==2) iov.Cl(occ==1)
#> 1 -0.28637459 -0.32602143 -0.18435100 -0.06883430 -0.11106234
#> 2 0.34495061 0.01464476 0.16817239 -0.04083867 -0.01601272
#> 3 0.24400619 0.19866998 -0.04173507 -0.29629002 -0.13689397
#> 4 0.02743539 0.05776011 -0.24078054 0.04593929 0.02620149
#> 5 0.19090509 -0.16212950 0.49114346 0.10336253 -0.13556974
#> 6 -0.15327550 -0.09805926 0.11269492 -0.01609979 0.02763238
#> iov.Cl(occ==2) iov.Ka(occ==1) iov.Ka(occ==2) V2 V3 TCl
#> 1 0.13902745 0.003094429 0.03594814 40.16754 297.1583 18.46715
#> 2 -0.02666789 0.007220282 -0.06790929 40.16754 297.1583 18.46715
#> 3 0.11351119 -0.086691401 -0.06195628 40.16754 297.1583 18.46715
#> 4 0.01901208 -0.112879021 0.18574793 40.16754 297.1583 18.46715
#> 5 -0.03040496 -0.006146741 0.03085491 40.16754 297.1583 18.46715
#> 6 -0.12782711 0.043930861 -0.08498676 40.16754 297.1583 18.46715
#> eta.Cl TKA eta.Ka Q Kin Kout EC50
#> 1 -0.40789384 0.3212433 -0.07189073 10.3295 0.9591919 0.9804559 200.2417
#> 2 0.11494752 0.3212433 0.30166435 10.3295 0.9591919 0.9804559 200.2417
#> 3 -0.07089084 0.3212433 -0.31349280 10.3295 0.9591919 0.9804559 200.2417
#> 4 0.25526236 0.3212433 -0.67862826 10.3295 0.9591919 0.9804559 200.2417
#> 5 -0.42579910 0.3212433 -0.18178978 10.3295 0.9591919 0.9804559 200.2417
#> 6 -0.24397398 0.3212433 0.28921473 10.3295 0.9591919 0.9804559 200.2417For 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.