[F] Add SIMETA functionality

[kylebaron]

In response to @Eliford question.

This isn't it ...but sketching out how we could use $INCLUDE type capability to do this.

1. Link back to RcppArmadillo and mrgsolve
2. Added MVNORM plugin ... similar function that mrgsolve uses, but for now just draws one vector from multivariate normal
3. Need to put a check in $INCLUDE so that, when Rcpp and RcppArmadillo are included, we only include RcppArmadillo
4. And probably need a way to pass the OMEGA matrix to $MAIN.

library(dplyr)
library(magrittr)
library(rbenchmark)
library(mrgsolve)

code <- '
$GLOBAL
arma::mat foo = arma::mat(2,2);
arma::mat eta = arma::mat(1,2);

$INCLUDE MVNORM

$MAIN
if(NEWIND <= 1) {
  foo(0,0) = omegacl;
  foo(0,1) = cov;
  foo(1,0) = cov;
  foo(1,1) = omegav;
  eta(0,0) = 2;

  while(fabs(eta(0,0)) > 1.86) {
    eta = MVNORM(foo);
  }
}

$PARAM TVCL = 1, TVV = 20,
omegacl = 1, cov = 0.03, omegav = 1

$MAIN
double CL = TVCL*exp(eta(0,1));

$TABLE
table(eta1) = eta(0,0);
table(eta2) = eta(0,1);

'

mod <- mcode("test", code)

## Warning: Could not find a $INIT or $CMT block

out <- mod %>% mrgsim(nid=10000, end=-1) 

summary(out$eta1)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -1.856000 -0.620200 -0.013330 -0.001329  0.635800  1.859000

summary(out$eta2)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -3.759000 -0.668600  0.008806  0.009800  0.676600  3.810000

code <- '
$INCLUDE DIST

$PARAM TVCL = 1, TVV = 20, omegacl = 1

$CMT CENT

$MAIN

if(NEWIND <=1) {
  double ETA1 = 2;
  while(fabs(ETA1) > 1.86) {
    ETA1 = DIST::rnorm(0,sqrt(omegacl));  
  }
} 

double CL = TVCL*exp(ETA1);

double V =  TVV;

$PKMODEL ncmt=1

$CAPTURE ETA1
'

mod <- mcode("test2", code)

out <- mod %>% mrgsim(nid=10000, end=-1) 

summary(out$ETA1)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -1.85900 -0.60910  0.01784  0.01003  0.63450  1.86000

[kylebaron]

@Eliford I just checked this into development version. It's working on some basic tests right now ... but still pretty new to me too: no documentation and I don't know if or where it could break. Also the interface could change .... but I'll be sure to keep you in the loop if you want to kick this around. This is purely exploratory right now ... but if you want to help test, please take crack.

The important points are:

1. $INCLUDE SIMETA this will link your model to RcppArmadillo. It will take longer to compile
2. Right now the statement SIMETA; simulates a new set of ETAs. This might change at some point ... not sure yet. But this is the simplest I could get it.
3. Nothing automatic in SIMETA that will prevent you from going into an unending loop (for example, if the condition to stop the while loop was to get an ETA > 1000000 or something like that. Not sure if this will be a problem or not. Until then, I put in a some safety code that will only allow 100 draws for each ID.

library(dplyr)
library(magrittr)
library(rbenchmark)
library(mrgsolve)

code <- '

$INCLUDE SIMETA

$PARAM BASE = 20

$CMT GUT CENT

$MAIN

CENT_0 = BASE;

double x = 0;
double y = 0;
double z = 0;
unsigned int i = 0;


if(NEWIND <= 1) {
  z = 10;

  // The ETA before SIMETA
  x = ETA(1);
  SIMETA;
  // The ETA after SIMETA
  y = ETA(1);

  // z should be between -1.5 and 1.5
  z = 10;
  i = 0;
  while(fabs(z) > 1.5 && i < 100) {
    SIMETA;
    z = ETA(2); 
    i++;
  }
}


$OMEGA
1 1 1

$TABLE
table(x) = x;
table(y) = y;
table(z) = z;
table(zz) = ETA(3);
'

mod <- mcode("simeta", code,warn=FALSE)

system.time(out <- mod %>% req() %>% mrgsim(nid=50000, end=-1))

##    user  system elapsed 
##   0.438   0.007   0.445

out

## Model:  simeta.cpp 
## Dim:    50000 x 6 
## Time:   0 to 0 
## ID:     50000 
##      ID time         x        y       z       zz
## [1,]  1    0  0.977651 -2.30265  0.4010 -0.87622
## [2,]  2    0 -0.665674  0.17995  0.5649  1.58915
## [3,]  3    0  1.488937  1.06421  1.0868 -2.13114
## [4,]  4    0  0.599754 -0.01641 -0.1809 -0.28113
## [5,]  5    0 -0.915229 -0.36190  0.1221 -0.49882
## [6,]  6    0 -1.310028  2.82191  1.2454 -1.08793
## [7,]  7    0 -0.949348  1.39963  1.2207 -0.02824
## [8,]  8    0  0.008718  1.74547 -0.4339 -1.15992

out %>% as.tbl %>% distinct(ID) %>% summary

##        ID             time         x                   y            
##  Min.   :    1   Min.   :0   Min.   :-4.586760   Min.   :-4.181208  
##  1st Qu.:12501   1st Qu.:0   1st Qu.:-0.672077   1st Qu.:-0.676097  
##  Median :25000   Median :0   Median : 0.003194   Median :-0.002884  
##  Mean   :25000   Mean   :0   Mean   : 0.003424   Mean   :-0.002405  
##  3rd Qu.:37500   3rd Qu.:0   3rd Qu.: 0.677402   3rd Qu.: 0.672442  
##  Max.   :50000   Max.   :0   Max.   : 3.927683   Max.   : 4.245212  
##        z                   zz           
##  Min.   :-1.499868   Min.   :-3.988251  
##  1st Qu.:-0.570521   1st Qu.:-0.680557  
##  Median : 0.005256   Median :-0.001042  
##  Mean   : 0.002864   Mean   :-0.004165  
##  3rd Qu.: 0.575868   3rd Qu.: 0.667387  
##  Max.   : 1.499615   Max.   : 3.994862

[Eliford]

Thank you very much for working on this. I will use the development version now. I will post any issue I encounter.

[kylebaron]

@Eliford My apologies if you're already using the simeta stuff. There was a change in the interface in the current development version.

You can see an example here:
https://github.com/metrumresearchgroup/mrgsolve/wiki/plugin

I've been refining this over the past week and feel better about the long-term prospects for this sort of setup. Hopefully there will be more advanced features available from the same / consistent interface now. And I'm more confident now that the interface will look something like this going forward.

If you're already using, let me know if you need help with the new interface. I wrote the vignette to help explain some stuff. Eventually there will be better documentation in the user_guide around the corner.

[Eliford]

Thank you for improving the SIMETA functionality (now with plugin). It is more user friendly. However,I am yet to start using it.

[Eliford]

@kylebmetrum: This is just a feedback. SIMETA plugin is working great for me. I started using it over the weekend and I get what I expect.
Thank you sooo Much!