Skip to contents

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

rxode2 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         Kin
#> TKA   0.05379885  0.06084629 -0.044521758  0.012344708 -0.03631156  0.02449761
#> TCl   0.06084629  0.13030654 -0.130086996  0.013697570 -0.08071136 -0.01437930
#> V2   -0.04452176 -0.13008700  0.228490181  0.003044394  0.10291479  0.09704946
#> Q     0.01234471  0.01369757  0.003044394  0.055574557  0.01532365  0.03490612
#> V3   -0.03631156 -0.08071136  0.102914794  0.015323650  0.18823700  0.04485492
#> Kin   0.02449761 -0.01437930  0.097049457  0.034906120  0.04485492  0.11753936
#> Kout  0.03627966  0.06801639 -0.011657275  0.030516359 -0.06937874  0.01670058
#> EC50  0.02261036  0.04865038  0.005901667 -0.007439675 -0.02141328  0.04440321
#>             Kout         EC50
#> TKA   0.03627966  0.022610365
#> TCl   0.06801639  0.048650378
#> V2   -0.01165728  0.005901667
#> Q     0.03051636 -0.007439675
#> V3   -0.06937874 -0.021413281
#> Kin   0.01670058  0.044403211
#> Kout  0.17323877  0.011774338
#> EC50  0.01177434  0.094086997

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 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: 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)
#> 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'
#> -- 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.Cl~6 eye.K~7 eye.K~8
#>     <int> <fct>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>    <dbl>   <dbl>   <dbl>
#>  1      1 1      -0.171   0.283  -0.145 -0.0652 -0.144  -5.57e-2  0.153   0.430 
#>  2      1 2      -0.171   0.283  -0.145 -0.0652 -0.0267 -3.73e-1  0.438  -0.186 
#>  3      1 3      -0.171   0.283  -0.145 -0.0652  0.0563  2.89e-1 -0.151  -0.663 
#>  4      1 4      -0.171   0.283  -0.145 -0.0652  0.104   2.40e-1  0.0730  0.306 
#>  5      1 5      -0.171   0.283  -0.145 -0.0652 -0.574   8.93e-2  0.302   0.116 
#>  6      1 6      -0.171   0.283  -0.145 -0.0652  0.496   1.16e-1  0.0170 -0.381 
#>  7      1 7      -0.171   0.283  -0.145 -0.0652 -0.119  -2.38e-2  0.202  -0.0384
#>  8      1 8      -0.171   0.283  -0.145 -0.0652  0.304  -6.15e-4 -0.412   0.0316
#>  9      1 9      -0.171   0.283  -0.145 -0.0652 -0.322  -1.45e-1  0.276   0.0296
#> 10      1 10     -0.171   0.283  -0.145 -0.0652 -0.117  -2.03e-1 -0.327  -0.326 
#> # ... 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.171 -0.145 -0.144   0.153 -0.0370  0.278  0    0     
#> 2      1     1  0.1 -0.171 -0.145 -0.0557  0.430 -0.0370  0.278  3.39 0.0103
#> 3      1     1  4   -0.171 -0.145 -0.144   0.153 -0.0370  0.278 14.0  6.31  
#> 4      1     1  4.1 -0.171 -0.145 -0.0557  0.430 -0.0370  0.278 53.3  6.45  
#> 5      1     1  8   -0.171 -0.145 -0.144   0.153 -0.0370  0.278  1.80 9.47  
#> 6      1     1  8.1 -0.171 -0.145 -0.0557  0.430 -0.0370  0.278 14.9  9.50  
#> # ... 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:

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.1714654       0.282754     -0.1453517    -0.06517664
#> 2      1  2     -0.1714654       0.282754     -0.1453517    -0.06517664
#> 3      1  3     -0.1714654       0.282754     -0.1453517    -0.06517664
#> 4      1  4     -0.1714654       0.282754     -0.1453517    -0.06517664
#> 5      1  5     -0.1714654       0.282754     -0.1453517    -0.06517664
#> 6      1  6     -0.1714654       0.282754     -0.1453517    -0.06517664
#>   eye.Cl(eye==1) eye.Cl(eye==2) eye.Ka(eye==1) eye.Ka(eye==2) iov.Cl(occ==1)
#> 1    -0.14407698    -0.05566547     0.15264863      0.4302003    -0.03698016
#> 2    -0.02669572    -0.37309346     0.43839048     -0.1859381     0.07114705
#> 3     0.05632522     0.28903521    -0.15138282     -0.6625008     0.11211841
#> 4     0.10360968     0.24049667     0.07299190      0.3057251     0.12593732
#> 5    -0.57373099     0.08925899     0.30170989      0.1157607    -0.15709749
#> 6     0.49596789     0.11573179     0.01703636     -0.3811751     0.15971601
#>   iov.Cl(occ==2) iov.Ka(occ==1) iov.Ka(occ==2)       V2       V3      TCl
#> 1   -0.031395303     0.27755739     0.03961799 40.42595 297.0748 18.56616
#> 2   -0.082328119     0.13402419    -0.08399350 40.42595 297.0748 18.56616
#> 3   -0.001643445    -0.15837094    -0.01178689 40.42595 297.0748 18.56616
#> 4    0.054699698    -0.11279275     0.02410732 40.42595 297.0748 18.56616
#> 5   -0.020009758    -0.07987667     0.04242764 40.42595 297.0748 18.56616
#> 6    0.070971564     0.03901590    -0.04885075 40.42595 297.0748 18.56616
#>       eta.Cl       TKA     eta.Ka        Q      Kin     Kout    EC50
#> 1  0.2601271 0.3545425 -0.3758864 10.83319 1.296538 1.167748 199.978
#> 2 -0.3106944 0.3545425  0.2341073 10.83319 1.296538 1.167748 199.978
#> 3  0.2673731 0.3545425  0.5431282 10.83319 1.296538 1.167748 199.978
#> 4 -0.1175962 0.3545425 -0.3174643 10.83319 1.296538 1.167748 199.978
#> 5  0.2212713 0.3545425 -0.1173740 10.83319 1.296538 1.167748 199.978
#> 6 -0.0658035 0.3545425 -0.1737912 10.83319 1.296538 1.167748 199.978
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.06475492       0.186367     -0.1985893     -0.1511591
#> 2      2  2    -0.06475492       0.186367     -0.1985893     -0.1511591
#> 3      2  3    -0.06475492       0.186367     -0.1985893     -0.1511591
#> 4      2  4    -0.06475492       0.186367     -0.1985893     -0.1511591
#> 5      2  5    -0.06475492       0.186367     -0.1985893     -0.1511591
#> 6      2  6    -0.06475492       0.186367     -0.1985893     -0.1511591
#>   eye.Cl(eye==1) eye.Cl(eye==2) eye.Ka(eye==1) eye.Ka(eye==2) iov.Cl(occ==1)
#> 1    -0.37494372   0.2262433340    -0.27165305    0.065939808    0.008334125
#> 2     0.12104118  -0.5359501736     0.03259746   -0.587026496    0.068165001
#> 3     0.27610522  -0.2726351879    -0.24959602   -0.315299894    0.176437106
#> 4     0.26919772   0.2653386724    -0.05311869   -0.009057287    0.115410179
#> 5     0.03465298  -0.3612670385    -0.11797701   -0.284249056   -0.025605141
#> 6    -0.07812729  -0.0001582368    -0.01444489    0.152965614    0.252023096
#>   iov.Cl(occ==2) iov.Ka(occ==1) iov.Ka(occ==2)       V2       V3      TCl
#> 1    -0.10920478     0.12472419     0.07098692 40.60343 297.3651 18.10283
#> 2    -0.07432985    -0.11156728     0.02067388 40.60343 297.3651 18.10283
#> 3    -0.08753706    -0.02391582     0.11695665 40.60343 297.3651 18.10283
#> 4     0.03296760    -0.02576062     0.06976869 40.60343 297.3651 18.10283
#> 5     0.20362759     0.04315469    -0.08877331 40.60343 297.3651 18.10283
#> 6     0.06666725     0.07848039    -0.01381926 40.60343 297.3651 18.10283
#>        eta.Cl       TKA      eta.Ka        Q      Kin      Kout     EC50
#> 1  0.08085881 0.1119641  0.20564867 10.37129 1.022916 0.3254918 199.8773
#> 2 -0.02076663 0.1119641  0.04347798 10.37129 1.022916 0.3254918 199.8773
#> 3 -0.19128082 0.1119641  0.54529555 10.37129 1.022916 0.3254918 199.8773
#> 4 -0.61442709 0.1119641 -0.28786109 10.37129 1.022916 0.3254918 199.8773
#> 5  0.47745374 0.1119641  0.45806454 10.37129 1.022916 0.3254918 199.8773
#> 6  0.24892885 0.1119641  0.02002771 10.37129 1.022916 0.3254918 199.8773

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.