1 `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)

1.0.1 Simulate without `simeta`

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

.    user  system elapsed 
.   0.004   0.000   0.004

print(sum)

.        ID              time      RESPONSE          EBASE         
.  Min.   :  1.00   Min.   :0   Min.   : 13.45   Min.   :-2.00624  
.  1st Qu.: 25.75   1st Qu.:0   1st Qu.: 56.49   1st Qu.:-0.57110  
.  Median : 50.50   Median :0   Median :108.66   Median : 0.08301  
.  Mean   : 50.50   Mean   :0   Mean   :158.53   Mean   : 0.09402  
.  3rd Qu.: 75.25   3rd Qu.:0   3rd Qu.:198.08   3rd Qu.: 0.68328  
.  Max.   :100.00   Max.   :0   Max.   :739.21   Max.   : 2.00042  
.       BASE       
.  Min.   : 13.45  
.  1st Qu.: 56.49  
.  Median :108.66  
.  Mean   :158.53  
.  3rd Qu.:198.08  
.  Max.   :739.21

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

1.0.2 Simulate with `simeta`

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

.    user  system elapsed 
.   0.004   0.000   0.003

print(sum)

.        ID              time      RESPONSE          EBASE         
.  Min.   :  1.00   Min.   :0   Min.   : 80.77   Min.   :-0.21353  
.  1st Qu.: 25.75   1st Qu.:0   1st Qu.: 89.06   1st Qu.:-0.11582  
.  Median : 50.50   Median :0   Median : 96.87   Median :-0.03180  
.  Mean   : 50.50   Mean   :0   Mean   : 98.37   Mean   :-0.02274  
.  3rd Qu.: 75.25   3rd Qu.:0   3rd Qu.:108.52   3rd Qu.: 0.08178  
.  Max.   :100.00   Max.   :0   Max.   :119.58   Max.   : 0.17878  
.       BASE       
.  Min.   : 80.77  
.  1st Qu.: 89.06  
.  Median : 96.87  
.  Mean   : 98.37  
.  3rd Qu.:108.52  
.  Max.   :119.58

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

1.1 `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.

1.1.1 Simulate without `simeps`

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

.    user  system elapsed 
.   0.004   0.000   0.003

print(sum)

.        ID         time            CENT              CP         
.  Min.   :1   Min.   : 0.00   Min.   :  0.00   Min.   :-7.4847  
.  1st Qu.:1   1st Qu.:11.25   1st Qu.: 15.93   1st Qu.:-0.1713  
.  Median :1   Median :23.50   Median : 29.38   Median : 3.1417  
.  Mean   :1   Mean   :23.52   Mean   : 37.47   Mean   : 3.2181  
.  3rd Qu.:1   3rd Qu.:35.75   3rd Qu.: 54.21   3rd Qu.: 7.3496  
.  Max.   :1   Max.   :48.00   Max.   :100.00   Max.   :13.4523

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

1.1.2 Simulate with `simeps`

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

.    user  system elapsed 
.   0.004   0.000   0.003

print(sum)

.        ID         time            CENT              CP          
.  Min.   :1   Min.   : 0.00   Min.   :  0.00   Min.   : 0.03374  
.  1st Qu.:1   1st Qu.:11.25   1st Qu.: 15.93   1st Qu.: 2.50312  
.  Median :1   Median :23.50   Median : 29.38   Median : 5.67478  
.  Mean   :1   Mean   :23.52   Mean   : 37.47   Mean   : 6.85330  
.  3rd Qu.:1   3rd Qu.:35.75   3rd Qu.: 54.21   3rd Qu.:10.21087  
.  Max.   :1   Max.   :48.00   Max.   :100.00   Max.   :20.74302

Better … all concentrations are positive.

Resimulate random effect variates on demand

1 simeta example

1.0.1 Simulate without simeta

1.0.2 Simulate with simeta

1.1 simeps example