Nesting in rxode2
More than one level of nesting is possible in rxode2; 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 |
Event table
This event table contains nesting variables:
- inv: investigator id
- id: subject id
- eye: eye id (left or right)
- occ: occasion
#> rxode2 2.0.13.9000 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 ->
evrxode2 model
This creates the rxode2 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 <- rxode2({
## 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
})#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’Uncertainty in Model parameters
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 Kin
#> TKA 0.037214680 -0.01716076 -0.02124791 0.038978341 -0.01701485 -0.005130489
#> TCl -0.017160760 0.17260107 0.11885824 0.063639092 0.13671322 0.139484010
#> V2 -0.021247912 0.11885824 0.22814815 0.077206142 0.01832142 0.144048561
#> Q 0.038978341 0.06363909 0.07720614 0.141219272 0.01832059 0.111170832
#> V3 -0.017014852 0.13671322 0.01832142 0.018320587 0.21487246 0.072365954
#> Kin -0.005130489 0.13948401 0.14404856 0.111170832 0.07236595 0.317277559
#> Kout -0.006672981 0.05699381 0.06594624 0.001473466 0.02600052 0.068334568
#> EC50 -0.002670124 -0.06618161 -0.10807009 -0.070944904 0.01006890 -0.084067640
#> Kout EC50
#> TKA -0.006672981 -0.002670124
#> TCl 0.056993813 -0.066181613
#> V2 0.065946243 -0.108070094
#> Q 0.001473466 -0.070944904
#> V3 0.026000521 0.010068895
#> Kin 0.068334568 -0.084067640
#> Kout 0.053182657 -0.024078476
#> EC50 -0.024078476 0.110886887Nesting Variability
To 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 rxode2 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: nuUnexplained variability
The last piece of variability to specify is the unexplained variability
sigma <- lotri(prop.sd ~ .25,
add.sd~ 0.125)Solving the problem
s <- rxSolve(mod, theta, ev,
thetaMat=tMat, omega=omega,
sigma=sigma, sigmaDf=400,
nStud=400)#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’#> unhandled error message: EE:[lsoda] 70000 steps taken before reaching tout
#> @(lsoda.c:750#> Warning: some ID(s) could not solve the ODEs correctly; These values are
#> replaced with 'NA'
print(s)#> -- Solved rxode2 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.184 0.0339 -0.102
#> 2 1 2 0.184 0.0339 -0.102
#> 3 1 3 0.184 0.0339 -0.102
#> 4 1 4 0.184 0.0339 -0.102
#> 5 1 5 0.184 0.0339 -0.102
#> 6 1 6 0.184 0.0339 -0.102
#> 7 1 7 0.184 0.0339 -0.102
#> 8 1 8 0.184 0.0339 -0.102
#> 9 1 9 0.184 0.0339 -0.102
#> 10 1 10 0.184 0.0339 -0.102
#> # i 7,990 more rows
#> # i 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 ($thetaMat for changes)
#> * omega matrix ($omegaList)
#>
#> -- First part of data (object): --
#> # A tibble: 976,000 x 21
#> 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.184 -0.102 0.102 0.203 0.0284 0.241 0 0
#> 2 1 1 0.1 0.184 -0.102 -0.0211 -0.230 0.0284 0.241 8.62 0.0143
#> 3 1 1 4 0.184 -0.102 0.102 0.203 0.0284 0.241 43.1 4.25
#> 4 1 1 4.1 0.184 -0.102 -0.0211 -0.230 0.0284 0.241 34.7 4.34
#> 5 1 1 8 0.184 -0.102 0.102 0.203 0.0284 0.241 21.5 6.60
#> 6 1 1 8.1 0.184 -0.102 -0.0211 -0.230 0.0284 0.241 12.4 6.63
#> # i 975,994 more rows
#> # i 10 more variables: CL <dbl>, KA <dbl>, ef0 <dbl>, depot <dbl>, centr <dbl>,
#> # peri <dbl>, eff <dbl>, occ <fct>, eye <fct>, inv <fct>There are multiple investigators in a study; Each investigator has a
number of individuals enrolled at their site. rxode2
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.1840224 0.03394184 -0.1016289 0.2547654
#> 2 1 2 0.1840224 0.03394184 -0.1016289 0.2547654
#> 3 1 3 0.1840224 0.03394184 -0.1016289 0.2547654
#> 4 1 4 0.1840224 0.03394184 -0.1016289 0.2547654
#> 5 1 5 0.1840224 0.03394184 -0.1016289 0.2547654
#> 6 1 6 0.1840224 0.03394184 -0.1016289 0.2547654
#> eye.Cl(eye==1) eye.Cl(eye==2) eye.Ka(eye==1) eye.Ka(eye==2) iov.Cl(occ==1)
#> 1 0.10194347 -0.021115812 0.20337983 -0.2302690 0.02835565
#> 2 0.01130175 -0.021694714 -0.15154228 0.1358393 -0.06714322
#> 3 0.13707474 0.003585042 0.37430930 -0.3717166 -0.08521717
#> 4 -0.31073977 0.041739806 0.07300206 0.2223848 -0.06683794
#> 5 0.08384790 0.166605861 -0.13491960 -0.5548553 -0.09729366
#> 6 -0.01973425 0.024876177 -0.63509849 -0.3476314 0.06885542
#> iov.Cl(occ==2) iov.Ka(occ==1) iov.Ka(occ==2) V2 V3 TCl
#> 1 -0.005235509 0.24110229 0.04864245 39.77126 296.848 18.33307
#> 2 0.066536064 -0.02414637 -0.01071090 39.77126 296.848 18.33307
#> 3 -0.069816325 0.09914816 -0.02590684 39.77126 296.848 18.33307
#> 4 0.029893498 0.12992391 0.03215196 39.77126 296.848 18.33307
#> 5 -0.236222049 -0.01972638 0.08233437 39.77126 296.848 18.33307
#> 6 -0.153375815 0.02119405 -0.08195238 39.77126 296.848 18.33307
#> eta.Cl TKA eta.Ka Q Kin Kout EC50
#> 1 0.16708104 0.2729526 0.09068134 10.35831 0.3258594 0.6571626 200.0506
#> 2 -0.44146336 0.2729526 -0.29424356 10.35831 0.3258594 0.6571626 200.0506
#> 3 -0.06417486 0.2729526 -0.65670740 10.35831 0.3258594 0.6571626 200.0506
#> 4 0.18001278 0.2729526 0.46182787 10.35831 0.3258594 0.6571626 200.0506
#> 5 -0.14624908 0.2729526 0.20874497 10.35831 0.3258594 0.6571626 200.0506
#> 6 -0.43753621 0.2729526 -0.18666467 10.35831 0.3258594 0.6571626 200.0506#> sim.id id inv.Cl(inv==1) inv.Cl(inv==2) inv.Ka(inv==1) inv.Ka(inv==2)
#> 1 2 1 0.0805846 -0.01008132 -0.2526698 -0.04032135
#> 2 2 2 0.0805846 -0.01008132 -0.2526698 -0.04032135
#> 3 2 3 0.0805846 -0.01008132 -0.2526698 -0.04032135
#> 4 2 4 0.0805846 -0.01008132 -0.2526698 -0.04032135
#> 5 2 5 0.0805846 -0.01008132 -0.2526698 -0.04032135
#> 6 2 6 0.0805846 -0.01008132 -0.2526698 -0.04032135
#> eye.Cl(eye==1) eye.Cl(eye==2) eye.Ka(eye==1) eye.Ka(eye==2) iov.Cl(occ==1)
#> 1 -0.33748468 0.1629759 -0.1381456 -0.302243666 -0.00483812
#> 2 0.31142002 -0.2269723 0.1955361 -0.040894579 -0.24099209
#> 3 0.12619236 0.3489916 -0.3590797 0.088988973 0.02381477
#> 4 0.02725676 -0.1552299 0.0496031 -0.002903132 -0.06909340
#> 5 0.04460059 0.2733171 -0.2682419 -0.077056515 0.20675893
#> 6 0.05392161 0.1719959 0.3147468 0.216212260 -0.03419832
#> iov.Cl(occ==2) iov.Ka(occ==1) iov.Ka(occ==2) V2 V3 TCl
#> 1 0.234360899 0.048689271 -0.079185363 40.41734 297.2159 18.3585
#> 2 0.041281347 -0.007535906 0.038300729 40.41734 297.2159 18.3585
#> 3 0.109920971 -0.055068923 0.004376570 40.41734 297.2159 18.3585
#> 4 -0.007221718 -0.053405622 -0.092849121 40.41734 297.2159 18.3585
#> 5 0.148717727 -0.038365864 0.035105208 40.41734 297.2159 18.3585
#> 6 0.154494032 -0.145763596 0.000685484 40.41734 297.2159 18.3585
#> eta.Cl TKA eta.Ka Q Kin Kout EC50
#> 1 0.194063184 0.3235566 -0.08273807 10.09662 0.358013 1.112678 200.3166
#> 2 -0.007519027 0.3235566 -0.61521523 10.09662 0.358013 1.112678 200.3166
#> 3 0.254727124 0.3235566 0.15476468 10.09662 0.358013 1.112678 200.3166
#> 4 -0.105072455 0.3235566 0.06888678 10.09662 0.358013 1.112678 200.3166
#> 5 0.504153787 0.3235566 -0.12321080 10.09662 0.358013 1.112678 200.3166
#> 6 -0.215040153 0.3235566 -0.73844046 10.09662 0.358013 1.112678 200.3166For 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.