Resimulate random effect variates on demand

library(mrgsolve)
library(dplyr)
options(mrgsolve.soloc="build")

`simeta` example

In this example, we want to simulate a patient-specific baseline response that is between 80 and 120.
In the code, we start a loop that calls simeta with no arguments until the baseline is between the specified bounds
For this example, we only calculate BASE when NEWIND <=1 … or whenever we are working on the first record of an individual. This ensures that we don’t re-simulate BASE at every simulation record.
We have also implemented a counter (i) that ensures we only try to resimulate 100 times. If a value for BASE cannot be generated in 100 tries, we give up.
You probably won’t need to implement FLAG for your problem. I am only using FLAG here so we can selectively call the simeta code to demonstrate how it is working.

code <- '
$PARAM TVBASE = 100, FLAG = 0

$CMT RESPONSE

$MAIN 

if(NEWIND <=1) {

  capture BASE = TVBASE*exp(EBASE);

  int i = 0;

  if(FLAG > 0) {
    while((BASE < 80) || (BASE > 120)) {
      if(++i > 100) {
        mrg::report("There was a problem simulating BASE.");
        break;
      }
      simeta();
      BASE = TVBASE*exp(EBASE);
    }
  }
  RESPONSE_0 = BASE;
}


$OMEGA @labels EBASE
1

$CAPTURE EBASE
'

mod <- mcode_cache("simeta",  code)

Simulate without `simeta`

system.time({
  out <- mod %>% mrgsim(nid=100, end=-1)
  sum <- summary(out)
})

.    user  system elapsed 
.   0.006   0.000   0.006

print(sum)

.        ID              time      RESPONSE            EBASE         
.  Min.   :  1.00   Min.   :0   Min.   :   9.529   Min.   :-2.35082  
.  1st Qu.: 25.75   1st Qu.:0   1st Qu.:  53.307   1st Qu.:-0.62910  
.  Median : 50.50   Median :0   Median : 101.460   Median : 0.01447  
.  Mean   : 50.50   Mean   :0   Mean   : 168.258   Mean   : 0.07127  
.  3rd Qu.: 75.25   3rd Qu.:0   3rd Qu.: 236.023   3rd Qu.: 0.85873  
.  Max.   :100.00   Max.   :0   Max.   :1006.399   Max.   : 2.30896  
.       BASE         
.  Min.   :   9.529  
.  1st Qu.:  53.307  
.  Median : 101.460  
.  Mean   : 168.258  
.  3rd Qu.: 236.023  
.  Max.   :1006.399

When we simulate with FLAG=0, the simeta code isn’t called and we BASE values all over the place.

Simulate with `simeta`

system.time({
  out <- mod %>% mrgsim(nid=100, end=-1, param=list(FLAG=1))
  sum <- summary(out)
})

.    user  system elapsed 
.   0.001   0.000   0.001

print(sum)

.        ID              time      RESPONSE          EBASE         
.  Min.   :  1.00   Min.   :0   Min.   : 80.78   Min.   :-0.21349  
.  1st Qu.: 25.75   1st Qu.:0   1st Qu.: 87.69   1st Qu.:-0.13137  
.  Median : 50.50   Median :0   Median : 96.81   Median :-0.03245  
.  Mean   : 50.50   Mean   :0   Mean   : 97.81   Mean   :-0.02905  
.  3rd Qu.: 75.25   3rd Qu.:0   3rd Qu.:107.45   3rd Qu.: 0.07187  
.  Max.   :100.00   Max.   :0   Max.   :119.59   Max.   : 0.17888  
.       BASE       
.  Min.   : 80.78  
.  1st Qu.: 87.69  
.  Median : 96.81  
.  Mean   : 97.81  
.  3rd Qu.:107.45  
.  Max.   :119.59

When we simulate with FLAG=1, the simeta code is called and we BASE values within the specified bounds.

`simeps` example

In this example, we want to re-simulate the residual error variate to make sure we have a concentration that is positive.
We set up a loop that looks like the simeta example, but we work in $TABLE this time because we are calculating CP.

code <- '
$PARAM CL = 1, V = 20, FLAG=0

$SIGMA 50

$PKMODEL cmt="CENT"

$TABLE
capture CP = CENT/V + EPS(1);

int i = 0;
while(CP < 0 && FLAG > 0) {
  if(++i > 100) {
    mrg::report("Problem simulating positive CP");
    break;
  }
  simeps();
  CP = CENT/V + EPS(1);
}

'

mod <- mcode_cache("simeps", code)

Building simeps ... done.

Simulate without `simeps`

system.time({
  out <- mod %>% ev(amt=100) %>% mrgsim(end=48)
  sum <- summary(out)
})

.    user  system elapsed 
.   0.003   0.000   0.003

print(sum)

.        ID         time            CENT              CP          
.  Min.   :1   Min.   : 0.00   Min.   :  0.00   Min.   :-15.4846  
.  1st Qu.:1   1st Qu.:11.25   1st Qu.: 15.93   1st Qu.: -1.4798  
.  Median :1   Median :23.50   Median : 29.38   Median :  0.9079  
.  Mean   :1   Mean   :23.52   Mean   : 37.47   Mean   :  2.0886  
.  3rd Qu.:1   3rd Qu.:35.75   3rd Qu.: 54.21   3rd Qu.:  6.3176  
.  Max.   :1   Max.   :48.00   Max.   :100.00   Max.   : 15.8507

Negative concentrations are simulated when we don’t call the simeps loop.

Simulate with `simeps`

system.time({
  out <- mod %>% ev(amt=100) %>% mrgsim(end=48, param=list(FLAG=1))
  sum <- summary(out)
})

   user  system elapsed 
  0.013   0.000   0.013

print(sum)

       ID         time            CENT              CP         
 Min.   :1   Min.   : 0.00   Min.   :  0.00   Min.   : 0.1761  
 1st Qu.:1   1st Qu.:11.25   1st Qu.: 15.93   1st Qu.: 1.9950  
 Median :1   Median :23.50   Median : 29.38   Median : 4.8258  
 Mean   :1   Mean   :23.52   Mean   : 37.47   Mean   : 5.6284  
 3rd Qu.:1   3rd Qu.:35.75   3rd Qu.: 54.21   3rd Qu.: 8.4775  
 Max.   :1   Max.   :48.00   Max.   :100.00   Max.   :16.6980

Better … all concentrations are positive.

simeta example

Simulate without simeta

Simulate with simeta

simeps example

Simulate without simeps