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.12 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
})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
#> TKA 0.143263775 0.032638666 0.002765035 0.091854052 0.022512898
#> TCl 0.032638666 0.127002940 -0.014574890 0.114719467 0.075347713
#> V2 0.002765035 -0.014574890 0.096211535 -0.058107773 0.008753197
#> Q 0.091854052 0.114719467 -0.058107773 0.180588829 0.072189358
#> V3 0.022512898 0.075347713 0.008753197 0.072189358 0.220977991
#> Kin -0.050665239 -0.006714258 0.060492544 -0.039674513 0.138031230
#> Kout 0.009982166 0.061994646 -0.021442701 0.053535835 -0.007252253
#> EC50 0.018161679 -0.035483756 -0.052638269 -0.004006376 -0.029344993
#> Kin Kout EC50
#> TKA -0.050665239 0.009982166 0.018161679
#> TCl -0.006714258 0.061994646 -0.035483756
#> V2 0.060492544 -0.021442701 -0.052638269
#> Q -0.039674513 0.053535835 -0.004006376
#> V3 0.138031230 -0.007252253 -0.029344993
#> Kin 0.174879222 -0.053948772 -0.035743773
#> Kout -0.053948772 0.085935965 -0.023561776
#> EC50 -0.035743773 -0.023561776 0.095329782Nesting 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)#> 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.C~1 inv.C~2 inv.K~3 inv.K~4 eye.C~5 eye.C~6 eye.Ka~7 eye.K~8
#> <int> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 1 0.183 -0.494 0.0475 0.0396 0.261 0.0620 -0.162 0.395
#> 2 1 2 0.183 -0.494 0.0475 0.0396 0.0542 0.0688 0.230 -0.0537
#> 3 1 3 0.183 -0.494 0.0475 0.0396 -0.0977 0.501 -0.149 -0.175
#> 4 1 4 0.183 -0.494 0.0475 0.0396 -0.193 -0.134 0.00757 0.126
#> 5 1 5 0.183 -0.494 0.0475 0.0396 0.0261 -0.196 -0.160 0.126
#> 6 1 6 0.183 -0.494 0.0475 0.0396 -0.543 -0.0827 0.0768 0.0935
#> 7 1 7 0.183 -0.494 0.0475 0.0396 -0.215 -0.0186 -0.0442 0.492
#> 8 1 8 0.183 -0.494 0.0475 0.0396 0.0103 0.197 0.249 0.273
#> 9 1 9 0.183 -0.494 0.0475 0.0396 -0.312 -0.419 0.0205 0.0238
#> 10 1 10 0.183 -0.494 0.0475 0.0396 -0.0504 -0.421 -0.0986 0.176
#> # ... with 7,990 more rows, 14 more variables: `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>, and abbreviated variable names
#> # 1: `inv.Cl(inv==1)`, 2: `inv.Cl(inv==2)`, 3: `inv.Ka(inv==1)`,
#> # 4: `inv.Ka(inv==2)`, 5: `eye.Cl(eye==1)`, 6: `eye.Cl(eye==2)`,
#> # 7: `eye.Ka(eye==1)`, 8: `eye.Ka(eye==2)`
#> -- 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.183 0.0475 0.261 -0.162 0.0698 -0.0543 0 0
#> 2 1 1 0.1 0.183 0.0475 0.0620 0.395 0.0698 -0.0543 7.94 0.0105
#> 3 1 1 4 0.183 0.0475 0.261 -0.162 0.0698 -0.0543 8.36 7.29
#> 4 1 1 4.1 0.183 0.0475 0.0620 0.395 0.0698 -0.0543 69.0 7.43
#> 5 1 1 8 0.183 0.0475 0.261 -0.162 0.0698 -0.0543 16.0 9.89
#> 6 1 1 8.1 0.183 0.0475 0.0620 0.395 0.0698 -0.0543 16.7 9.90
#> # ... with 975,994 more rows, and 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.1828788 -0.49416 0.04748327 0.03958731
#> 2 1 2 0.1828788 -0.49416 0.04748327 0.03958731
#> 3 1 3 0.1828788 -0.49416 0.04748327 0.03958731
#> 4 1 4 0.1828788 -0.49416 0.04748327 0.03958731
#> 5 1 5 0.1828788 -0.49416 0.04748327 0.03958731
#> 6 1 6 0.1828788 -0.49416 0.04748327 0.03958731
#> eye.Cl(eye==1) eye.Cl(eye==2) eye.Ka(eye==1) eye.Ka(eye==2) iov.Cl(occ==1)
#> 1 0.26141829 0.06202272 -0.161521508 0.39454650 0.069799681
#> 2 0.05416414 0.06882105 0.230048974 -0.05368092 -0.003263889
#> 3 -0.09772280 0.50069918 -0.148895693 -0.17451566 0.024280279
#> 4 -0.19277842 -0.13433412 0.007565801 0.12647856 -0.027506157
#> 5 0.02613360 -0.19578933 -0.160099675 0.12571014 0.006267818
#> 6 -0.54338345 -0.08267844 0.076824293 0.09347925 0.046881622
#> iov.Cl(occ==2) iov.Ka(occ==1) iov.Ka(occ==2) V2 V3 TCl
#> 1 -0.009559067 -0.05428404 -0.15376368 40.25594 296.6028 18.25838
#> 2 -0.202023239 0.07924445 -0.25508621 40.25594 296.6028 18.25838
#> 3 0.164674628 0.03314495 0.03869573 40.25594 296.6028 18.25838
#> 4 -0.082827559 0.11986668 -0.03384211 40.25594 296.6028 18.25838
#> 5 -0.005249714 0.07687374 0.06823757 40.25594 296.6028 18.25838
#> 6 -0.104984851 -0.24464808 -0.08148293 40.25594 296.6028 18.25838
#> eta.Cl TKA eta.Ka Q Kin Kout EC50
#> 1 -0.29692069 0.3551825 -0.30186930 10.3325 0.803143 1.073291 200.0174
#> 2 0.09406344 0.3551825 0.15140058 10.3325 0.803143 1.073291 200.0174
#> 3 0.08693360 0.3551825 -0.34674767 10.3325 0.803143 1.073291 200.0174
#> 4 0.05726492 0.3551825 0.50194619 10.3325 0.803143 1.073291 200.0174
#> 5 0.24499407 0.3551825 -0.06303613 10.3325 0.803143 1.073291 200.0174
#> 6 -0.32278229 0.3551825 -0.56998453 10.3325 0.803143 1.073291 200.0174#> sim.id id inv.Cl(inv==1) inv.Cl(inv==2) inv.Ka(inv==1) inv.Ka(inv==2)
#> 1 2 1 -0.0626712 0.04608162 0.2914937 0.1598596
#> 2 2 2 -0.0626712 0.04608162 0.2914937 0.1598596
#> 3 2 3 -0.0626712 0.04608162 0.2914937 0.1598596
#> 4 2 4 -0.0626712 0.04608162 0.2914937 0.1598596
#> 5 2 5 -0.0626712 0.04608162 0.2914937 0.1598596
#> 6 2 6 -0.0626712 0.04608162 0.2914937 0.1598596
#> eye.Cl(eye==1) eye.Cl(eye==2) eye.Ka(eye==1) eye.Ka(eye==2) iov.Cl(occ==1)
#> 1 -0.47066149 -0.07791196 -0.01788343 0.4801737 -0.05451005
#> 2 0.02699322 0.03402449 -0.36996275 0.1539896 0.03384445
#> 3 -0.21239472 -0.08143373 -0.20084445 -0.1065339 -0.08395143
#> 4 0.05363381 0.26774906 0.15334080 -0.1657511 0.04232139
#> 5 0.18275346 -0.08879176 0.19565336 0.3499795 -0.01870560
#> 6 -0.27624928 -0.05645154 0.08980599 0.3289861 -0.14739417
#> iov.Cl(occ==2) iov.Ka(occ==1) iov.Ka(occ==2) V2 V3 TCl
#> 1 0.07876265 -0.03815685 -0.06341663 40.01196 296.9316 18.50566
#> 2 -0.09701967 0.09904382 0.12186930 40.01196 296.9316 18.50566
#> 3 -0.11566919 -0.05425941 0.06129329 40.01196 296.9316 18.50566
#> 4 -0.20088668 -0.06478261 -0.05139012 40.01196 296.9316 18.50566
#> 5 0.08216436 0.13550242 -0.11181202 40.01196 296.9316 18.50566
#> 6 0.07764325 -0.08838865 -0.13538536 40.01196 296.9316 18.50566
#> eta.Cl TKA eta.Ka Q Kin Kout EC50
#> 1 -0.2304925 0.5525777 0.193522870 10.47074 0.753052 0.8732414 200.4423
#> 2 -0.1029807 0.5525777 0.239879773 10.47074 0.753052 0.8732414 200.4423
#> 3 0.2189875 0.5525777 0.004535902 10.47074 0.753052 0.8732414 200.4423
#> 4 -0.1789055 0.5525777 -0.120521524 10.47074 0.753052 0.8732414 200.4423
#> 5 0.4077604 0.5525777 -0.163146034 10.47074 0.753052 0.8732414 200.4423
#> 6 -0.3168620 0.5525777 0.176973472 10.47074 0.753052 0.8732414 200.4423For 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.