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 2.1.2.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 <- 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.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=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'

print(s)

#> -- 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.149           0.0666          -0.0608
#>  2      1 2                0.149           0.0666          -0.0608
#>  3      1 3                0.149           0.0666          -0.0608
#>  4      1 4                0.149           0.0666          -0.0608
#>  5      1 5                0.149           0.0666          -0.0608
#>  6      1 6                0.149           0.0666          -0.0608
#>  7      1 7                0.149           0.0666          -0.0608
#>  8      1 8                0.149           0.0666          -0.0608
#>  9      1 9                0.149           0.0666          -0.0608
#> 10      1 10               0.149           0.0666          -0.0608
#> # 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.149 -0.0608  0.0137 -0.105 -0.0627 -0.0588  0    0     
#> 2      1     1  0.1  0.149 -0.0608 -0.0758 -0.157 -0.0627 -0.0588  5.35 0.0108
#> 3      1     1  4    0.149 -0.0608  0.0137 -0.105 -0.0627 -0.0588 46.4  5.61  
#> 4      1     1  4.1  0.149 -0.0608 -0.0758 -0.157 -0.0627 -0.0588 45.7  5.75  
#> 5      1     1  8    0.149 -0.0608  0.0137 -0.105 -0.0627 -0.0588 23.7  9.29  
#> 6      1     1  8.1  0.149 -0.0608 -0.0758 -0.157 -0.0627 -0.0588 39.5  9.34  
#> # 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.1492479      0.0665656    -0.06079813    0.004733472
#> 2      1  2      0.1492479      0.0665656    -0.06079813    0.004733472
#> 3      1  3      0.1492479      0.0665656    -0.06079813    0.004733472
#> 4      1  4      0.1492479      0.0665656    -0.06079813    0.004733472
#> 5      1  5      0.1492479      0.0665656    -0.06079813    0.004733472
#> 6      1  6      0.1492479      0.0665656    -0.06079813    0.004733472
#>   eye.Cl(eye==1) eye.Cl(eye==2) eye.Ka(eye==1) eye.Ka(eye==2) iov.Cl(occ==1)
#> 1    0.013722304    -0.07575309    -0.10537117   -0.157208354    -0.06269727
#> 2   -0.310461174    -0.32965956    -0.32748735   -0.063160091     0.14493514
#> 3    0.004236311    -0.08037476     0.09903067    0.007048495    -0.02303274
#> 4    0.185836927    -0.20310126     0.28466188   -0.201217795    -0.04494725
#> 5   -0.057520782     0.27098390     0.17417693    0.206112899    -0.06608225
#> 6   -0.338946095     0.18769487     0.03454896   -0.029638991     0.02065775
#>   iov.Cl(occ==2) iov.Ka(occ==1) iov.Ka(occ==2)       V2       V3      TCl
#> 1    -0.11214693   -0.058819964    0.094094101 40.13563 297.4763 17.89084
#> 2     0.09591295   -0.053492080   -0.088904917 40.13563 297.4763 17.89084
#> 3    -0.23254972    0.091325853    0.007085379 40.13563 297.4763 17.89084
#> 4     0.17245303    0.003427374    0.015688993 40.13563 297.4763 17.89084
#> 5    -0.04083866    0.063690605   -0.076675254 40.13563 297.4763 17.89084
#> 6     0.09025181    0.157316987   -0.075983567 40.13563 297.4763 17.89084
#>        eta.Cl       TKA      eta.Ka        Q      Kin     Kout     EC50
#> 1 -0.12215654 0.3084588  0.03304006 10.55699 1.535729 1.129002 199.8009
#> 2 -0.06450068 0.3084588  0.29264958 10.55699 1.535729 1.129002 199.8009
#> 3 -0.01511460 0.3084588 -0.49100154 10.55699 1.535729 1.129002 199.8009
#> 4 -0.10174580 0.3084588  0.04016224 10.55699 1.535729 1.129002 199.8009
#> 5  0.43976675 0.3084588 -0.07847621 10.55699 1.535729 1.129002 199.8009
#> 6 -0.06100883 0.3084588 -0.48437792 10.55699 1.535729 1.129002 199.8009

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.03797007   -0.003500523     0.09175848   -0.001407258
#> 2      2  2    -0.03797007   -0.003500523     0.09175848   -0.001407258
#> 3      2  3    -0.03797007   -0.003500523     0.09175848   -0.001407258
#> 4      2  4    -0.03797007   -0.003500523     0.09175848   -0.001407258
#> 5      2  5    -0.03797007   -0.003500523     0.09175848   -0.001407258
#> 6      2  6    -0.03797007   -0.003500523     0.09175848   -0.001407258
#>   eye.Cl(eye==1) eye.Cl(eye==2) eye.Ka(eye==1) eye.Ka(eye==2) iov.Cl(occ==1)
#> 1    -0.14051326    -0.29144557      0.4051760    -0.03631610    -0.18082050
#> 2    -0.11565333     0.33346437      0.3163433    -0.03394341     0.10221965
#> 3     0.05882822     0.08843513     -0.4566997    -0.20076567     0.10053381
#> 4    -0.04654595    -0.02775519      0.2206998     0.27772712     0.06816234
#> 5     0.17343983    -0.03893844      0.1984094    -0.07858389     0.08071230
#> 6    -0.05639225     0.09009346      0.3991894     0.34541555    -0.09126549
#>   iov.Cl(occ==2) iov.Ka(occ==1) iov.Ka(occ==2)       V2       V3      TCl
#> 1    -0.03998668   -0.085438687    -0.08477932 40.07637 297.2877 19.59749
#> 2    -0.04647827   -0.142005266    -0.04867860 40.07637 297.2877 19.59749
#> 3     0.02113789   -0.049955478     0.04501613 40.07637 297.2877 19.59749
#> 4    -0.09208714   -0.025512110    -0.13304296 40.07637 297.2877 19.59749
#> 5     0.15775875   -0.039322706     0.06220148 40.07637 297.2877 19.59749
#> 6     0.06104810    0.007667834    -0.10824797 40.07637 297.2877 19.59749
#>        eta.Cl       TKA     eta.Ka        Q       Kin     Kout     EC50
#> 1 -0.16234260 0.2005236 -0.3951311 10.56206 0.4232414 1.716212 199.8744
#> 2 -0.03181361 0.2005236  0.2448144 10.56206 0.4232414 1.716212 199.8744
#> 3 -0.03800938 0.2005236 -0.4449192 10.56206 0.4232414 1.716212 199.8744
#> 4  0.51756138 0.2005236 -0.1933814 10.56206 0.4232414 1.716212 199.8744
#> 5 -0.35070400 0.2005236 -0.2593210 10.56206 0.4232414 1.716212 199.8744
#> 6  0.35870375 0.2005236 -0.3610563 10.56206 0.4232414 1.716212 199.8744

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.

2024-04-16