Reproducible results with set.seed

example

This post demonstrates that random simulation outputs can be made reproducible by using set.seed().

Author

Kyle Baron

Published

01/01/2017

library(mrgsolve)
library(dplyr)

1 Demo model

code <- '
$PKMODEL cmt="CENT"
$PARAM CL = 1, V = 20

$OMEGA 1 1

$TABLE
capture ETA1 = ETA(1);
capture ETA2 = ETA(2);
'
mod <- mcode("test", code)

2 Single dose

set.seed(9911)
single <- mod %>% ev(amt=100) %>% mrgsim(nid=100, end=24)

3 Multi dose

set.seed(9911)
multi <- mod %>% ev(amt=100, ii=24, addl=9) %>% mrgsim(nid=100, end=240)

4 Compare

single
. Model:  test 
. Dim:    2600 x 5 
. Time:   0 to 24 
. ID:     100 
.     ID time   CENT ETA1   ETA2
. 1:   1    0   0.00 1.75 -1.404
. 2:   1    0 100.00 1.75 -1.404
. 3:   1    1  95.12 1.75 -1.404
. 4:   1    2  90.48 1.75 -1.404
. 5:   1    3  86.07 1.75 -1.404
. 6:   1    4  81.87 1.75 -1.404
. 7:   1    5  77.88 1.75 -1.404
. 8:   1    6  74.08 1.75 -1.404
multi
. Model:  test 
. Dim:    24200 x 5 
. Time:   0 to 240 
. ID:     100 
.     ID time   CENT ETA1   ETA2
. 1:   1    0   0.00 1.75 -1.404
. 2:   1    0 100.00 1.75 -1.404
. 3:   1    1  95.12 1.75 -1.404
. 4:   1    2  90.48 1.75 -1.404
. 5:   1    3  86.07 1.75 -1.404
. 6:   1    4  81.87 1.75 -1.404
. 7:   1    5  77.88 1.75 -1.404
. 8:   1    6  74.08 1.75 -1.404
filter(single, CENT==0)
. # A tibble: 100 × 5
.       ID  time  CENT     ETA1   ETA2
.    <dbl> <dbl> <dbl>    <dbl>  <dbl>
.  1     1     0     0  1.75    -1.40 
.  2     2     0     0  0.00190 -0.417
.  3     3     0     0  1.26     1.88 
.  4     4     0     0  0.554    0.360
.  5     5     0     0  0.863   -2.79 
.  6     6     0     0 -0.913   -0.223
.  7     7     0     0 -1.12     0.548
.  8     8     0     0 -0.108   -1.19 
.  9     9     0     0 -0.405    0.660
. 10    10     0     0 -1.24    -1.06 
. # … with 90 more rows
filter(multi, CENT==0)
. # A tibble: 100 × 5
.       ID  time  CENT     ETA1   ETA2
.    <dbl> <dbl> <dbl>    <dbl>  <dbl>
.  1     1     0     0  1.75    -1.40 
.  2     2     0     0  0.00190 -0.417
.  3     3     0     0  1.26     1.88 
.  4     4     0     0  0.554    0.360
.  5     5     0     0  0.863   -2.79 
.  6     6     0     0 -0.913   -0.223
.  7     7     0     0 -1.12     0.548
.  8     8     0     0 -0.108   -1.19 
.  9     9     0     0 -0.405    0.660
. 10    10     0     0 -1.24    -1.06 
. # … with 90 more rows

4.1 Identical in the first day …

identical(filter(single, time < 24),filter(multi,time < 24))
. [1] TRUE

4.2 … but not after that

filter(single, time > 24)
. # A tibble: 0 × 5
. # … with 5 variables: ID <dbl>, time <dbl>, CENT <dbl>, ETA1 <dbl>, ETA2 <dbl>
filter(multi, time > 24)
. # A tibble: 21,600 × 5
.       ID  time  CENT  ETA1  ETA2
.    <dbl> <dbl> <dbl> <dbl> <dbl>
.  1     1    25 124.   1.75 -1.40
.  2     1    26 118.   1.75 -1.40
.  3     1    27 112.   1.75 -1.40
.  4     1    28 107.   1.75 -1.40
.  5     1    29 101.   1.75 -1.40
.  6     1    30  96.4  1.75 -1.40
.  7     1    31  91.7  1.75 -1.40
.  8     1    32  87.2  1.75 -1.40
.  9     1    33  83.0  1.75 -1.40
. 10     1    34  78.9  1.75 -1.40
. # … with 21,590 more rows