Nesting in RxODE

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

Event table

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 ->
  ev

RxODE model

This 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.

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.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.035658813

Nesting 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 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: nu

Unexplained 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)
#> 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'
#> ______________________________ 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:

print(head(s$params))
#>   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
print(head(s$params %>% filter(sim.id == 2)))
#>   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.9159

For 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.