Nesting levels in rxode2

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

set.seed(1014)

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

#> unhandled error message: EE:[lsoda] 70000 steps taken before reaching tout
#>  @(lsoda.c:752)
#> unhandled error message: EE:[lsoda] 70000 steps taken before reaching tout
#>  @(lsoda.c:752)
#> unhandled error message: EE:[lsoda] 70000 steps taken before reaching tout
#>  @(lsoda.c:752)
#> unhandled error message: EE:[lsoda] 70000 steps taken before reaching tout
#>  @(lsoda.c:752)
#> unhandled error message: EE:[lsoda] 70000 steps taken before reaching tout
#>  @(lsoda.c:752)
#> unhandled error message: EE:[lsoda] 70000 steps taken before reaching tout
#>  @(lsoda.c:752)
#> unhandled error message: EE:[lsoda] 70000 steps taken before reaching tout
#>  @(lsoda.c:752)
#> unhandled error message: EE:[lsoda] 70000 steps taken before reaching tout
#>  @(lsoda.c:752)
#> unhandled error message: EE:[lsoda] 70000 steps taken before reaching tout
#>  @(lsoda.c:752)
#> unhandled error message: EE:[lsoda] 70000 steps taken before reaching tout
#>  @(lsoda.c:752)
#> unhandled error message: EE:[lsoda] 70000 steps taken before reaching tout
#>  @(lsoda.c:752)
#> unhandled error message: EE:[lsoda] 70000 steps taken before reaching tout
#>  @(lsoda.c:752)
#> unhandled error message: EE:[lsoda] 70000 steps taken before reaching tout
#>  @(lsoda.c:752)

#> 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.153           0.0701
#>  2      1 2               -0.125           -0.153           0.0701
#>  3      1 3               -0.125           -0.153           0.0701
#>  4      1 4               -0.125           -0.153           0.0701
#>  5      1 5               -0.125           -0.153           0.0701
#>  6      1 6               -0.125           -0.153           0.0701
#>  7      1 7               -0.125           -0.153           0.0701
#>  8      1 8               -0.125           -0.153           0.0701
#>  9      1 9               -0.125           -0.153           0.0701
#> 10      1 10              -0.125           -0.153           0.0701
#> # 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>, 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: 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.0701 -0.177 -0.273 -0.0212 -0.0875  0     0      
#> 2      1     1  0.1 -0.125 0.0701  0.132  0.394 -0.0212 -0.0875 -0.197 0.00724
#> 3      1     1  4   -0.125 0.0701 -0.177 -0.273 -0.0212 -0.0875 61.6   5.76   
#> 4      1     1  4.1 -0.125 0.0701  0.132  0.394 -0.0212 -0.0875 51.0   5.88   
#> 5      1     1  8   -0.125 0.0701 -0.177 -0.273 -0.0212 -0.0875 15.9   8.18   
#> 6      1     1  8.1 -0.125 0.0701  0.132  0.394 -0.0212 -0.0875 19.8   8.20   
#> # 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.1529692     0.07014459     0.04335434
#> 2      1  2     -0.1254652     -0.1529692     0.07014459     0.04335434
#> 3      1  3     -0.1254652     -0.1529692     0.07014459     0.04335434
#> 4      1  4     -0.1254652     -0.1529692     0.07014459     0.04335434
#> 5      1  5     -0.1254652     -0.1529692     0.07014459     0.04335434
#> 6      1  6     -0.1254652     -0.1529692     0.07014459     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.131706240    -0.27302477     0.39404583    -0.02123662
#> 2   -0.009839091   -0.217949757     0.31190569     0.38603090     0.04989188
#> 3   -0.124445741    0.477197022    -0.12809443    -0.13919935    -0.21847189
#> 4   -0.193074978    0.001424707    -0.18425474     0.08945107    -0.11929660
#> 5    0.156461380    0.518130536    -0.01610728    -0.18191604     0.01041120
#> 6    0.193470464   -0.129242249     0.06840428     0.23350313     0.11968596
#>   iov.Cl(occ==2) iov.Ka(occ==1) iov.Ka(occ==2)       V2       V3      TCl
#> 1   -0.046015932   -0.087454725    -0.07362632 40.73185 296.8681 18.50917
#> 2    0.008798674   -0.123778069     0.07342570 40.73185 296.8681 18.50917
#> 3   -0.102271286    0.005686142     0.15971079 40.73185 296.8681 18.50917
#> 4    0.068260157   -0.124760703     0.07412590 40.73185 296.8681 18.50917
#> 5   -0.114937885    0.085510837     0.03127980 40.73185 296.8681 18.50917
#> 6   -0.010831963    0.084431763     0.09175564 40.73185 296.8681 18.50917
#>        eta.Cl       TKA      eta.Ka        Q      Kin      Kout     EC50
#> 1  0.22676346 0.2354877 -0.06421990 10.20557 0.650145 0.6719848 200.3501
#> 2  0.01720316 0.2354877 -0.15204055 10.20557 0.650145 0.6719848 200.3501
#> 3 -0.12666922 0.2354877 -0.04664825 10.20557 0.650145 0.6719848 200.3501
#> 4 -0.11380876 0.2354877 -0.36471000 10.20557 0.650145 0.6719848 200.3501
#> 5 -0.29184818 0.2354877 -0.24201355 10.20557 0.650145 0.6719848 200.3501
#> 6  0.02707820 0.2354877 -0.25173810 10.20557 0.650145 0.6719848 200.3501

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

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.