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.1.2 using 2 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
})
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’

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.1114888313  0.080564822  0.00342595  0.100977903 -0.0002587565
#> TCl   0.0805648220  0.358725754 -0.03900873  0.081456063  0.0127513557
#> V2    0.0034259499 -0.039008731  0.16186347  0.021455884  0.0901981860
#> Q     0.1009779031  0.081456063  0.02145588  0.118507891 -0.0074253634
#> V3   -0.0002587565  0.012751356  0.09019819 -0.007425363  0.0962043635
#> Kin  -0.0345435870 -0.019640912  0.04347588 -0.014156924  0.0631010458
#> Kout -0.0054267621  0.002847269  0.03538826  0.041898784 -0.0498996575
#> EC50 -0.0257430877 -0.036836442  0.02598046 -0.035344383  0.0194233764
#>               Kin         Kout        EC50
#> TKA  -0.034543587 -0.005426762 -0.02574309
#> TCl  -0.019640912  0.002847269 -0.03683644
#> V2    0.043475876  0.035388256  0.02598046
#> Q    -0.014156924  0.041898784 -0.03534438
#> V3    0.063101046 -0.049899658  0.01942338
#> Kin   0.154494363 -0.005810415  0.01795091
#> Kout -0.005810415  0.163452698  0.01746295
#> EC50  0.017950912  0.017462948  0.08392045

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)
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> unhandled error message: EE:[lsoda] 70000 steps taken before reaching tout
#>  @(lsoda.c:751
#> 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.Cl(inv==1)` `inv.Cl(inv==2)` `inv.Ka(inv==1)`
#>     <int> <fct>            <dbl>            <dbl>            <dbl>
#>  1      1 1              -0.0109          -0.0679            0.407
#>  2      1 2              -0.0109          -0.0679            0.407
#>  3      1 3              -0.0109          -0.0679            0.407
#>  4      1 4              -0.0109          -0.0679            0.407
#>  5      1 5              -0.0109          -0.0679            0.407
#>  6      1 6              -0.0109          -0.0679            0.407
#>  7      1 7              -0.0109          -0.0679            0.407
#>  8      1 8              -0.0109          -0.0679            0.407
#>  9      1 9              -0.0109          -0.0679            0.407
#> 10      1 10             -0.0109          -0.0679            0.407
#> # i 7,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: 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.0109  0.407  0.119 0.0298 0.0766 0.0414   0    0     
#> 2      1     1  0.1 -0.0109  0.407  0.210 0.0251 0.0766 0.0414  30.2  0.0514
#> 3      1     1  4   -0.0109  0.407  0.119 0.0298 0.0766 0.0414  53.8 11.0   
#> 4      1     1  4.1 -0.0109  0.407  0.210 0.0251 0.0766 0.0414  21.8 11.1   
#> 5      1     1  8   -0.0109  0.407  0.119 0.0298 0.0766 0.0414  11.6 11.7   
#> 6      1     1  8.1 -0.0109  0.407  0.210 0.0251 0.0766 0.0414  10.8 11.7   
#> # i 975,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.0108638    -0.06785478      0.4072602     -0.1817774
#> 2      1  2     -0.0108638    -0.06785478      0.4072602     -0.1817774
#> 3      1  3     -0.0108638    -0.06785478      0.4072602     -0.1817774
#> 4      1  4     -0.0108638    -0.06785478      0.4072602     -0.1817774
#> 5      1  5     -0.0108638    -0.06785478      0.4072602     -0.1817774
#> 6      1  6     -0.0108638    -0.06785478      0.4072602     -0.1817774
#>   eye.Cl(eye==1) eye.Cl(eye==2) eye.Ka(eye==1) eye.Ka(eye==2) iov.Cl(occ==1)
#> 1     0.11873491     0.20997720    0.029765309     0.02506394    0.076596752
#> 2     0.06691781    -0.16346135   -0.143332809     0.14465129    0.003167789
#> 3     0.06028580    -0.29573229   -0.087171501     0.13719106   -0.173910706
#> 4     0.13723694    -0.01996766    0.336740151     0.09382886   -0.068865861
#> 5    -0.07910533     0.28493577    0.018776938    -0.11316552    0.137662442
#> 6    -0.46238288     0.18217218    0.008058002    -0.01851742   -0.244161771
#>   iov.Cl(occ==2) iov.Ka(occ==1) iov.Ka(occ==2)       V2       V3      TCl
#> 1    -0.01810611    0.041428006   -0.103783885 39.87406 297.0258 19.14037
#> 2     0.05041521    0.284705060   -0.001486575 39.87406 297.0258 19.14037
#> 3    -0.07260717   -0.079264146   -0.037515597 39.87406 297.0258 19.14037
#> 4     0.06555440   -0.001528273    0.013843836 39.87406 297.0258 19.14037
#> 5    -0.01615194   -0.016123395    0.048488620 39.87406 297.0258 19.14037
#> 6     0.02486997    0.117116580    0.083084290 39.87406 297.0258 19.14037
#>        eta.Cl       TKA      eta.Ka        Q       Kin      Kout     EC50
#> 1 -0.45918070 0.8313329 -0.15619340 10.90492 0.8587868 0.5779411 199.9883
#> 2 -0.55695799 0.8313329  0.36025072 10.90492 0.8587868 0.5779411 199.9883
#> 3 -0.23351796 0.8313329  0.05065181 10.90492 0.8587868 0.5779411 199.9883
#> 4 -0.02409583 0.8313329 -0.34464102 10.90492 0.8587868 0.5779411 199.9883
#> 5 -0.06957227 0.8313329  0.25833714 10.90492 0.8587868 0.5779411 199.9883
#> 6  0.08344246 0.8313329 -0.29313409 10.90492 0.8587868 0.5779411 199.9883
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.1271798    -0.09660709    -0.08707027     0.09813521
#> 2      2  2      0.1271798    -0.09660709    -0.08707027     0.09813521
#> 3      2  3      0.1271798    -0.09660709    -0.08707027     0.09813521
#> 4      2  4      0.1271798    -0.09660709    -0.08707027     0.09813521
#> 5      2  5      0.1271798    -0.09660709    -0.08707027     0.09813521
#> 6      2  6      0.1271798    -0.09660709    -0.08707027     0.09813521
#>   eye.Cl(eye==1) eye.Cl(eye==2) eye.Ka(eye==1) eye.Ka(eye==2) iov.Cl(occ==1)
#> 1     -0.7612991     0.21460130    -0.14170136      0.2185514   -0.133352219
#> 2     -0.3509146     0.29965733    -0.02441195     -0.4906191   -0.004736529
#> 3      0.3515331    -0.05861322     0.06739830     -0.4111834    0.011317771
#> 4     -0.1919729    -0.29530910    -0.34831020     -0.1778014   -0.080863900
#> 5      0.3858965    -0.03144946     0.18695232     -0.3648079   -0.197356015
#> 6      0.2135519    -0.09750608     0.26602736      0.3183694    0.123762037
#>   iov.Cl(occ==2) iov.Ka(occ==1) iov.Ka(occ==2)       V2       V3      TCl
#> 1     0.03591827   -0.190188676    -0.04532262 39.76378 296.9413 18.79856
#> 2    -0.09693353    0.006120717     0.05658106 39.76378 296.9413 18.79856
#> 3    -0.12074594   -0.131954560    -0.11378586 39.76378 296.9413 18.79856
#> 4     0.15911676   -0.145297471     0.01248529 39.76378 296.9413 18.79856
#> 5     0.01885821   -0.026566485    -0.07063365 39.76378 296.9413 18.79856
#> 6     0.04808696    0.013804412     0.02895194 39.76378 296.9413 18.79856
#>          eta.Cl       TKA      eta.Ka        Q      Kin      Kout     EC50
#> 1 -0.4034608778 0.3669201  0.54424611 10.65578 1.818653 0.8634281 199.7095
#> 2 -0.0584168782 0.3669201  0.17676029 10.65578 1.818653 0.8634281 199.7095
#> 3  0.0002969481 0.3669201 -0.39414255 10.65578 1.818653 0.8634281 199.7095
#> 4  0.2375981756 0.3669201  0.04203328 10.65578 1.818653 0.8634281 199.7095
#> 5  0.1480406727 0.3669201 -0.56547145 10.65578 1.818653 0.8634281 199.7095
#> 6 -0.5661562486 0.3669201 -0.21037180 10.65578 1.818653 0.8634281 199.7095

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.