metrumresearchgroup / mrgsolve

Simulate from ODE-based population PK/PD and QSP models in R
https://mrgsolve.org
GNU General Public License v2.0
131 stars 36 forks source link

Considering simulating OMEGA regardless of etasrc #1138

Closed kylebaron closed 9 months ago

kylebaron commented 11 months ago

The issue is that EPS isn't reproducible when simulating the same problem with etasrc set to "OMEGA" and "data.all", for example. There isn't a bug here, it's just the random number generator isn't getting tapped in the same way to make EPS reproducible regardless of what we do with ETA.

Example of EPS that is not reproducible when switchin etasrc

library(mrgsolve)
#> 
#> Attaching package: 'mrgsolve'
#> The following object is masked from 'package:stats':
#> 
#>     filter

data <- expand.ev(
  amt = 100, 
  ETA1 = 0.1, ETA2 = 0.2, ETA3 = 0.3
)

mod <- modlib("popex", capture = "EPS(1)")
#> Building popex ...
#> done.
mod <- smat(mod, matrix(1))

set.seed(112)
out1 <- mrgsim(mod, data)
set.seed(112)
out2 <- mrgsim(mod, data, etasrc = "data.all")

out1
#> Model:  popex 
#> Dim:    482 x 10 
#> Time:   0 to 240 
#> ID:     1 
#>     ID time    GUT  CENT    CL     V   ECL  IPRED      DV   EPS_1
#> 1:   1  0.0   0.00  0.00 1.458 22.52 0.377 0.0000  0.0000 -0.9965
#> 2:   1  0.0 100.00  0.00 1.458 22.52 0.377 0.0000  0.0000 -0.2696
#> 3:   1  0.5  81.24 18.45 1.458 22.52 0.377 0.8193  0.2427 -1.2165
#> 4:   1  1.0  66.00 32.85 1.458 22.52 0.377 1.4588  0.7215 -0.7040
#> 5:   1  1.5  53.61 43.99 1.458 22.52 0.377 1.9531  4.1797  0.7608
#> 6:   1  2.0  43.55 52.48 1.458 22.52 0.377 2.3301  2.9877  0.2486
#> 7:   1  2.5  35.38 58.84 1.458 22.52 0.377 2.6128 19.0116  1.9846
#> 8:   1  3.0  28.74 63.50 1.458 22.52 0.377 2.8194 17.3175  1.8152
out2
#> Model:  popex 
#> Dim:    482 x 10 
#> Time:   0 to 240 
#> ID:     1 
#>     ID time    GUT  CENT    CL     V ECL  IPRED     DV   EPS_1
#> 1:   1  0.0   0.00  0.00 1.105 29.31 0.1 0.0000 0.0000 -0.2012
#> 2:   1  0.0 100.00  0.00 1.105 29.31 0.1 0.0000 0.0000  0.6884
#> 3:   1  0.5  71.36 28.36 1.105 29.31 0.1 0.9674 0.7448 -0.2615
#> 4:   1  1.0  50.92 48.07 1.105 29.31 0.1 1.6397 1.1106 -0.3896
#> 5:   1  1.5  36.33 61.61 1.105 29.31 0.1 2.1017 0.7759 -0.9965
#> 6:   1  2.0  25.93 70.76 1.105 29.31 0.1 2.4140 1.8435 -0.2696
#> 7:   1  2.5  18.50 76.79 1.105 29.31 0.1 2.6197 0.7762 -1.2165
#> 8:   1  3.0  13.20 80.61 1.105 29.31 0.1 2.7498 1.3600 -0.7040

Created on 2024-02-03 with reprex v2.0.2

Now it's reproducible

library(mrgsolve)
#> 
#> Attaching package: 'mrgsolve'
#> The following object is masked from 'package:stats':
#> 
#>     filter

data <- expand.ev(
  amt = 100, 
  ETA1 = 0.1, ETA2 = 0.2, ETA3 = 0.3
)

mod <- modlib("popex", capture = "EPS(1)")
#> Building popex ...
#> done.
mod <- smat(mod, matrix(1))

set.seed(112)
out1 <- mrgsim(mod, data)
set.seed(112)
out2 <- mrgsim(mod, data, etasrc = "data.all")

out1
#> Model:  popex 
#> Dim:    482 x 10 
#> Time:   0 to 240 
#> ID:     1 
#>     ID time    GUT  CENT    CL     V   ECL  IPRED      DV   EPS_1
#> 1:   1  0.0   0.00  0.00 1.458 22.52 0.377 0.0000  0.0000 -0.9965
#> 2:   1  0.0 100.00  0.00 1.458 22.52 0.377 0.0000  0.0000 -0.2696
#> 3:   1  0.5  81.24 18.45 1.458 22.52 0.377 0.8193  0.2427 -1.2165
#> 4:   1  1.0  66.00 32.85 1.458 22.52 0.377 1.4588  0.7215 -0.7040
#> 5:   1  1.5  53.61 43.99 1.458 22.52 0.377 1.9531  4.1797  0.7608
#> 6:   1  2.0  43.55 52.48 1.458 22.52 0.377 2.3301  2.9877  0.2486
#> 7:   1  2.5  35.38 58.84 1.458 22.52 0.377 2.6128 19.0116  1.9846
#> 8:   1  3.0  28.74 63.50 1.458 22.52 0.377 2.8194 17.3175  1.8152
out2
#> Model:  popex 
#> Dim:    482 x 10 
#> Time:   0 to 240 
#> ID:     1 
#>     ID time    GUT  CENT    CL     V ECL  IPRED      DV   EPS_1
#> 1:   1  0.0   0.00  0.00 1.105 29.31 0.1 0.0000  0.0000 -0.9965
#> 2:   1  0.0 100.00  0.00 1.105 29.31 0.1 0.0000  0.0000 -0.2696
#> 3:   1  0.5  71.36 28.36 1.105 29.31 0.1 0.9674  0.2866 -1.2165
#> 4:   1  1.0  50.92 48.07 1.105 29.31 0.1 1.6397  0.8110 -0.7040
#> 5:   1  1.5  36.33 61.61 1.105 29.31 0.1 2.1017  4.4977  0.7608
#> 6:   1  2.0  25.93 70.76 1.105 29.31 0.1 2.4140  3.0952  0.2486
#> 7:   1  2.5  18.50 76.79 1.105 29.31 0.1 2.6197 19.0622  1.9846
#> 8:   1  3.0  13.20 80.61 1.105 29.31 0.1 2.7498 16.8896  1.8152

Created on 2024-02-03 with reprex v2.0.2