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.038   0.000   0.038
print(sum)
.        ID              time      RESPONSE           EBASE         
.  Min.   :  1.00   Min.   :0   Min.   :  11.90   Min.   :-2.12900  
.  1st Qu.: 25.75   1st Qu.:0   1st Qu.:  54.34   1st Qu.:-0.60985  
.  Median : 50.50   Median :0   Median : 115.32   Median : 0.14254  
.  Mean   : 50.50   Mean   :0   Mean   : 177.58   Mean   : 0.08727  
.  3rd Qu.: 75.25   3rd Qu.:0   3rd Qu.: 195.23   3rd Qu.: 0.66872  
.  Max.   :100.00   Max.   :0   Max.   :1448.53   Max.   : 2.67313  
.       BASE        
.  Min.   :  11.90  
.  1st Qu.:  54.34  
.  Median : 115.32  
.  Mean   : 177.58  
.  3rd Qu.: 195.23  
.  Max.   :1448.53

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.005   0.000   0.004
print(sum)
.        ID              time      RESPONSE          EBASE         
.  Min.   :  1.00   Min.   :0   Min.   : 80.73   Min.   :-0.21404  
.  1st Qu.: 25.75   1st Qu.:0   1st Qu.: 87.86   1st Qu.:-0.12940  
.  Median : 50.50   Median :0   Median : 96.07   Median :-0.04005  
.  Mean   : 50.50   Mean   :0   Mean   : 97.30   Mean   :-0.03396  
.  3rd Qu.: 75.25   3rd Qu.:0   3rd Qu.:105.14   3rd Qu.: 0.05008  
.  Max.   :100.00   Max.   :0   Max.   :119.90   Max.   : 0.18152  
.       BASE       
.  Min.   : 80.73  
.  1st Qu.: 87.86  
.  Median : 96.07  
.  Mean   : 97.30  
.  3rd Qu.:105.14  
.  Max.   :119.90

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)
Loading model from cache.

Simulate without simeps

system.time({
  out <- mod %>% ev(amt=100) %>% mrgsim(end=48)
  sum <- summary(out)
})
.    user  system elapsed 
.   0.008   0.000   0.008
print(sum)
.        ID         time            CENT              CP         
.  Min.   :1   Min.   : 0.00   Min.   :  0.00   Min.   :-14.460  
.  1st Qu.:1   1st Qu.:11.25   1st Qu.: 15.93   1st Qu.: -3.555  
.  Median :1   Median :23.50   Median : 29.38   Median :  2.544  
.  Mean   :1   Mean   :23.52   Mean   : 37.47   Mean   :  1.528  
.  3rd Qu.:1   3rd Qu.:35.75   3rd Qu.: 54.21   3rd Qu.:  6.951  
.  Max.   :1   Max.   :48.00   Max.   :100.00   Max.   : 15.472

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.034   0.000   0.034 
print(sum)
       ID         time            CENT              CP          
 Min.   :1   Min.   : 0.00   Min.   :  0.00   Min.   : 0.03585  
 1st Qu.:1   1st Qu.:11.25   1st Qu.: 15.93   1st Qu.: 3.71079  
 Median :1   Median :23.50   Median : 29.38   Median : 6.69697  
 Mean   :1   Mean   :23.52   Mean   : 37.47   Mean   : 7.02250  
 3rd Qu.:1   3rd Qu.:35.75   3rd Qu.: 54.21   3rd Qu.: 9.72255  
 Max.   :1   Max.   :48.00   Max.   :100.00   Max.   :18.61057  

Better … all concentrations are positive.