Nesting levels (IOV and beyond)

2026-10-27

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

library(rxode2)

#> rxode2 5.1.8.9000 using 2 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`

library(dplyr)

#> 
#> 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 <- function() {
  ini({
    TKA <- 0.294
    TCl <- 18.6
    V2 <- 40.2
    Q <- 10.5
    V3 <- 297
    Kin <- 1
    Kout <- 1
    EC50 <- 200
  })
  model({
    ## 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
  })
}
mod <- mod()

Uncertainty in Model parameters

set.seed(1014)

theta <- mod$theta

## 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.173571236 -0.1003204607  0.038185010 -0.004108928 -0.095032973
#> TCl  -0.100320461  0.2195710868 -0.043849095  0.013295549 -0.007477895
#> V2    0.038185010 -0.0438490948  0.129784612 -0.017270432 -0.038004762
#> Q    -0.004108928  0.0132955493 -0.017270432  0.022145634  0.020376451
#> V3   -0.095032973 -0.0074778948 -0.038004762  0.020376451  0.165568340
#> Kin  -0.040867119 -0.0492597458 -0.003056722 -0.033468634 -0.003021883
#> Kout  0.035469225  0.0275087955  0.033725901  0.027668205  0.005497301
#> EC50  0.026158042  0.0009434711  0.039426946 -0.036283167 -0.093134292
#>               Kin         Kout          EC50
#> TKA  -0.040867119  0.035469225  0.0261580416
#> TCl  -0.049259746  0.027508796  0.0009434711
#> V2   -0.003056722  0.033725901  0.0394269465
#> Q    -0.033468634  0.027668205 -0.0362831667
#> V3   -0.003021883  0.005497301 -0.0931342922
#> Kin   0.226735493 -0.083447793  0.0659884544
#> Kout -0.083447793  0.117195570 -0.0291684598
#> EC50  0.065988454 -0.029168460  0.0928611407

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, ev,
             thetaMat=tMat, omega=omega,
             sigma=sigma, sigmaDf=400,
             nStud=200)

#> [IDID=-1, excess work done on this call (perhaps wrong jt).]: 11 warning(s) for subject(s): 19 (sim 4), 1 (sim 43), 12 (sim 43), 19 (sim 122), 5 (sim 160), ... (6 more)
#> [intdy -- t = 182.256 +/- 39.8635 (88/240.1) illegal]: 164 warning(s) for subject(s): 18 (sim 43), 13 (sim 111), 12 (sim 186)
#> [lsoda -- internal t + h = t (h too small for machine precision)]: 756852 warning(s) for subject(s): 19 (sim 4), 4 (sim 31), 6 (sim 31), 7 (sim 31), 9 (sim 31), ... (74 more)

#> 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: 4,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.125           0.0701           -0.153
#>  2      1 2               -0.125           0.0701           -0.153
#>  3      1 3               -0.125           0.0701           -0.153
#>  4      1 4               -0.125           0.0701           -0.153
#>  5      1 5               -0.125           0.0701           -0.153
#>  6      1 6               -0.125           0.0701           -0.153
#>  7      1 7               -0.125           0.0701           -0.153
#>  8      1 8               -0.125           0.0701           -0.153
#>  9      1 9               -0.125           0.0701           -0.153
#> 10      1 10              -0.125           0.0701           -0.153
#> # i 3,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>, TKA <dbl>, TCl <dbl>, V2 <dbl>, Q <dbl>, V3 <dbl>,
#> #   Kin <dbl>, Kout <dbl>, EC50 <dbl>, eta.Cl <dbl>, eta.Ka <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: 488,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.125 -0.153 -0.177  0.132 -0.0212 -0.0460  0     0      
#> 2      1     1  0.1 -0.125 -0.153 -0.273  0.394 -0.0212 -0.0460 -0.236 0.00867
#> 3      1     1  4   -0.125 -0.153 -0.177  0.132 -0.0212 -0.0460 76.4   6.01   
#> 4      1     1  4.1 -0.125 -0.153 -0.273  0.394 -0.0212 -0.0460 63.9   6.16   
#> 5      1     1  8   -0.125 -0.153 -0.177  0.132 -0.0212 -0.0460 27.3   9.80   
#> 6      1     1  8.1 -0.125 -0.153 -0.273  0.394 -0.0212 -0.0460 34.1   9.85   
#> # i 487,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:

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.1254652     0.07014459     -0.1529692     0.04335434
#> 2      1  2     -0.1254652     0.07014459     -0.1529692     0.04335434
#> 3      1  3     -0.1254652     0.07014459     -0.1529692     0.04335434
#> 4      1  4     -0.1254652     0.07014459     -0.1529692     0.04335434
#> 5      1  5     -0.1254652     0.07014459     -0.1529692     0.04335434
#> 6      1  6     -0.1254652     0.07014459     -0.1529692     0.04335434
#>   eye.Cl(eye==1) eye.Cl(eye==2) eye.Ka(eye==1) eye.Ka(eye==2) iov.Cl(occ==1)
#> 1   -0.177088310    -0.27302477    0.131706240     0.39404583    -0.02123662
#> 2   -0.009839091     0.31190569   -0.217949757     0.38603090     0.04989188
#> 3   -0.124445741    -0.12809443    0.477197022    -0.13919935    -0.21847189
#> 4   -0.193074978    -0.18425474    0.001424707     0.08945107    -0.11929660
#> 5    0.156461380    -0.01610728    0.518130536    -0.18191604     0.01041120
#> 6    0.193470464     0.06840428   -0.129242249     0.23350313     0.11968596
#>   iov.Cl(occ==2) iov.Ka(occ==1) iov.Ka(occ==2)       TKA      TCl       V2
#> 1   -0.087454725   -0.046015932    -0.07362632 0.2354877 18.50917 40.73185
#> 2   -0.123778069    0.008798674     0.07342570 0.2354877 18.50917 40.73185
#> 3    0.005686142   -0.102271286     0.15971079 0.2354877 18.50917 40.73185
#> 4   -0.124760703    0.068260157     0.07412590 0.2354877 18.50917 40.73185
#> 5    0.085510837   -0.114937885     0.03127980 0.2354877 18.50917 40.73185
#> 6    0.084431763   -0.010831963     0.09175564 0.2354877 18.50917 40.73185
#>          Q       V3      Kin      Kout     EC50      eta.Cl      eta.Ka
#> 1 10.20557 296.8681 0.650145 0.6719848 200.3501  0.22676346 -0.06421990
#> 2 10.20557 296.8681 0.650145 0.6719848 200.3501  0.01720316 -0.15204055
#> 3 10.20557 296.8681 0.650145 0.6719848 200.3501 -0.12666922 -0.04664825
#> 4 10.20557 296.8681 0.650145 0.6719848 200.3501 -0.11380876 -0.36471000
#> 5 10.20557 296.8681 0.650145 0.6719848 200.3501 -0.29184818 -0.24201355
#> 6 10.20557 296.8681 0.650145 0.6719848 200.3501  0.02707820 -0.25173810

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.05438592      0.1504495      0.1738687     0.03488022
#> 2      2  2    -0.05438592      0.1504495      0.1738687     0.03488022
#> 3      2  3    -0.05438592      0.1504495      0.1738687     0.03488022
#> 4      2  4    -0.05438592      0.1504495      0.1738687     0.03488022
#> 5      2  5    -0.05438592      0.1504495      0.1738687     0.03488022
#> 6      2  6    -0.05438592      0.1504495      0.1738687     0.03488022
#>   eye.Cl(eye==1) eye.Cl(eye==2) eye.Ka(eye==1) eye.Ka(eye==2) iov.Cl(occ==1)
#> 1     0.42019765    -0.22563508     0.09755050   -0.449748606   -0.118757683
#> 2     0.13120948    -0.08428522    -0.05143635    0.197863209   -0.067081410
#> 3    -0.08106788    -0.19105418    -0.07627969    0.007062104    0.163568691
#> 4    -0.09836919    -0.17712473     0.58557232    0.086677902    0.018922047
#> 5     0.28692438     0.33372000    -0.07112331    0.072399824    0.006332733
#> 6     0.28255187    -0.23750571     0.03061533    0.525744685   -0.029028139
#>   iov.Cl(occ==2) iov.Ka(occ==1) iov.Ka(occ==2)       TKA     TCl       V2
#> 1   -0.074548652    0.145933281    0.032641505 0.5686675 18.2352 40.65042
#> 2    0.046241153   -0.161335163    0.060465829 0.5686675 18.2352 40.65042
#> 3   -0.121391365    0.006560036   -0.062182571 0.5686675 18.2352 40.65042
#> 4   -0.078563754    0.003162297    0.007408265 0.5686675 18.2352 40.65042
#> 5   -0.001768959    0.048703940    0.055628684 0.5686675 18.2352 40.65042
#> 6    0.197722775   -0.011293997   -0.112993469 0.5686675 18.2352 40.65042
#>          Q       V3      Kin      Kout     EC50      eta.Cl     eta.Ka
#> 1 10.33968 296.7639 1.612831 0.9919165 200.3462  0.68261056  0.5209415
#> 2 10.33968 296.7639 1.612831 0.9919165 200.3462 -0.19479212 -0.1309852
#> 3 10.33968 296.7639 1.612831 0.9919165 200.3462 -0.25227932 -0.3416660
#> 4 10.33968 296.7639 1.612831 0.9919165 200.3462  0.09409293  0.3663583
#> 5 10.33968 296.7639 1.612831 0.9919165 200.3462 -0.53525522 -0.2208024
#> 6 10.33968 296.7639 1.612831 0.9919165 200.3462  0.61207190  0.3556095

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.

Next steps

• Simulating populations
• Parameter uncertainty