Generating subjects and extracting subject-level information from a simulation

[john-harrold]

I have two related objectives. First take this example:

library(rxode2)

mod <- function() {
  ini({
    KA <- 2.94E-01
    TCl <- 1.86E+01
    # between subject variability
    eta.Cl ~ 0.4^2
    V2 <- 4.02E+01
    Q <- 1.05E+01
    V3 <- 2.97E+02
    Kin <- 1
    Kout <- 1
    EC50 <- 200
  })
  model({
    C2 <- centr/V2
    C3 <- peri/V3
    CL <-  TCl*exp(eta.Cl) ## This is coded as a variable in the model
    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
    eff(0) <- 1
  })
}


ev <- et(amountUnits="mg", timeUnits="hours") %>%
  et(amt=10000, cmt="centr")

sim  <- rxSolve(mod, ev, nSub=100) 

unique(sim$CL)

What is the best way to extract the individual subject parameters from this simulation? If I try something like this:

head(as.data.frame(sim))

It just gives me the time course of the pieces from the model section.

The second question is can I generate these subjects without running a simulation? Is there a lower-level function I can use to just generate a data frame of subjects?

[mattfidler]

This seems like a bug, if you have 2 etas it works and you can access it by $params

library(rxode2)
#> rxode2 2.1.0 using 8 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`

mod <- function() {
  ini({
    KA <- 2.94E-01
    TCl <- 1.86E+01
    # between subject variability
    eta.Cl ~ 0.4^2
    V2 <- 4.02E+01
    eta.V2 ~ 0.4^2
    Q <- 1.05E+01
    V3 <- 2.97E+02
    Kin <- 1
    Kout <- 1
    EC50 <- 200
  })
  model({
    C2 <- centr/(V2*exp(eta.V2))
    C3 <- peri/V3
    CL <-  TCl*exp(eta.Cl) ## This is coded as a variable in the model
    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
    eff(0) <- 1
  })
}


ev <- et(amountUnits="mg", timeUnits="hours") %>%
  et(amt=10000, cmt="centr")

sim  <- rxSolve(mod, ev, nSub=100, nStud=4)
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
print(sim)
#> ── Solved rxode2 object ──
#> ── Parameters ($params): ──
#> # A tibble: 400 × 11
#>    sim.id    KA   TCl    V2     Q    V3   Kin  Kout  EC50 eta.Cl  eta.V2
#>     <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>   <dbl>
#>  1      1 0.294  18.6  40.2  10.5   297     1     1   200 -0.230  0.0820
#>  2      2 0.294  18.6  40.2  10.5   297     1     1   200 -0.433  0.0675
#>  3      3 0.294  18.6  40.2  10.5   297     1     1   200 -0.669 -0.448 
#>  4      4 0.294  18.6  40.2  10.5   297     1     1   200  0.292  0.372 
#>  5      5 0.294  18.6  40.2  10.5   297     1     1   200 -0.493 -0.0845
#>  6      6 0.294  18.6  40.2  10.5   297     1     1   200  0.281  1.01  
#>  7      7 0.294  18.6  40.2  10.5   297     1     1   200  0.594  0.229 
#>  8      8 0.294  18.6  40.2  10.5   297     1     1   200 -0.349 -0.0942
#>  9      9 0.294  18.6  40.2  10.5   297     1     1   200  0.331 -0.246 
#> 10     10 0.294  18.6  40.2  10.5   297     1     1   200 -0.118 -0.433 
#> # ℹ 390 more rows
#> ── Initial Conditions ($inits): ──
#> depot centr  peri   eff 
#>     0     0     0     1 
#> 
#> Simulation without uncertainty in parameters, omega, or sigma matricies
#> 
#> ── First part of data (object): ──
#> # A tibble: 80,000 × 9
#>   sim.id  time    C2    C3    CL depot  centr  peri   eff
#>    <int>   [h] <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>
#> 1      1 0      229. 0      14.8     0 10000     0   1   
#> 2      1 0.121  214. 0.942  14.8     0  9326.  280.  1.06
#> 3      1 0.241  199. 1.82   14.8     0  8698.  539.  1.12
#> 4      1 0.362  186. 2.63   14.8     0  8114.  780.  1.17
#> 5      1 0.482  173. 3.38   14.8     0  7570. 1004.  1.21
#> 6      1 0.603  162. 4.08   14.8     0  7064. 1212.  1.25
#> # ℹ 79,994 more rows

^{Created on 2023-12-10 with reprex v2.0.2}

There may be a lower-level function to simulate (I think there may be an example even) but I can't recall it off hand.

[mattfidler]

Hm. I tried with one eta.Cl and it worked ok for me now:

library(rxode2)
#> rxode2 2.1.0 using 8 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`

mod <- function() {
  ini({
    KA <- 2.94E-01
    TCl <- 1.86E+01
    # between subject variability
    eta.Cl ~ 0.4^2
    V2 <- 4.02E+01
    Q <- 1.05E+01
    V3 <- 2.97E+02
    Kin <- 1
    Kout <- 1
    EC50 <- 200
  })
  model({
    C2 <- centr/(V2)
    C3 <- peri/V3
    CL <-  TCl*exp(eta.Cl) ## This is coded as a variable in the model
    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
    eff(0) <- 1
  })
}


ev <- et(amountUnits="mg", timeUnits="hours") %>%
  et(amt=10000, cmt="centr")

sim  <- rxSolve(mod, ev, nSub=100)
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
print(sim)
#> ── Solved rxode2 object ──
#> ── Parameters ($params): ──
#> # A tibble: 100 × 10
#>    sim.id    KA   TCl    V2     Q    V3   Kin  Kout  EC50 eta.Cl
#>     <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
#>  1      1 0.294  18.6  40.2  10.5   297     1     1   200 -0.379
#>  2      2 0.294  18.6  40.2  10.5   297     1     1   200 -0.108
#>  3      3 0.294  18.6  40.2  10.5   297     1     1   200  0.263
#>  4      4 0.294  18.6  40.2  10.5   297     1     1   200 -0.823
#>  5      5 0.294  18.6  40.2  10.5   297     1     1   200  0.101
#>  6      6 0.294  18.6  40.2  10.5   297     1     1   200  0.238
#>  7      7 0.294  18.6  40.2  10.5   297     1     1   200  0.745
#>  8      8 0.294  18.6  40.2  10.5   297     1     1   200 -0.257
#>  9      9 0.294  18.6  40.2  10.5   297     1     1   200  0.272
#> 10     10 0.294  18.6  40.2  10.5   297     1     1   200 -0.522
#> # ℹ 90 more rows
#> ── Initial Conditions ($inits): ──
#> depot centr  peri   eff 
#>     0     0     0     1 
#> 
#> Simulation without uncertainty in parameters, omega, or sigma matricies
#> 
#> ── First part of data (object): ──
#> # A tibble: 20,000 × 9
#>   sim.id  time    C2    C3    CL depot  centr  peri   eff
#>    <int>   [h] <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>
#> 1      1 0      249.  0     12.7     0 10000     0   1   
#> 2      1 0.121  232.  1.02  12.7     0  9328.  304.  1.06
#> 3      1 0.241  216.  1.97  12.7     0  8702.  586.  1.12
#> 4      1 0.362  202.  2.85  12.7     0  8119.  847.  1.18
#> 5      1 0.482  188.  3.67  12.7     0  7576. 1090.  1.22
#> 6      1 0.603  176.  4.43  12.7     0  7071. 1316.  1.26
#> # ℹ 19,994 more rows

^{Created on 2023-12-10 with reprex v2.0.2}

Could you try with the version on r-universe?

[mattfidler]

This is the documentation on how to simulate with prior data, which is what would be what you simulated:

https://nlmixr2.github.io/rxode2/articles/rxode2-prior-data.html

This is the lower level function you were probably looking for:

https://nlmixr2.github.io/rxode2/reference/rxSimThetaOmega.html

Here is an example bringing the two together:

https://nlmixr2.github.io/rxode2/articles/rxode2-pipeline.html?q=rxSimThetaOmega#simulating-from-a-data-frame-of-parameters

[john-harrold]

The example above now works with the r-universe install. I looked at the documentation for rxSimThetaOmega and I'm not sure how to use it ;). It doesn't seem very straight forward. I'll use the simulation method for now.