Is there a way to constrain sampling of ETA values?

[Eliford]

I have been searching for a way to constrain sampled ETA to within chosen extremes; similar to how it is done in NONMEM. IF(ICALL.EQ.4) THEN DO WHILE(ABS(ETA(1)).GT.1.86) CALL SIMETA(ETA) END DO ENDIF

Is this currently possible with mrgsolve?

[kylebaron]

• Short answer is: no SIMETA type functionality right now; it's on our radar and I'd expect it to be available at some point
• Two workarounds you could try now
  • If no off diagonals are involved, you can try something like the DIST example below. This is only in the development version and not documented.
  • If there are covariances involved or you just don't want to play with DIST, you could simulate ETA in R and tack them on to the data set (remember to put those names in $PARAM to mrgsolve will pick them up. There is an example here

Obviously the workarounds entail more work than I'd like. But was sort of waiting for someone to ask about this ... will see if we can make this happen sooner rather than later.

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

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("test", code)

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

summary(out$ETA1)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -1.859000 -0.631800 -0.007315 -0.005622  0.626100  1.859000

devtools::session_info()

##  setting  value                       
##  version  R version 3.3.0 (2016-05-03)
##  system   x86_64, darwin13.4.0        
##  ui       X11                         
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  tz       America/Chicago             
##  date     2016-05-22                  
## 
##  package    * version    date       source                            
##  assertthat   0.1        2013-11-08 local                             
##  DBI          0.4-1      2016-05-08 CRAN (R 3.3.0)                    
##  devtools     1.10.0     2016-01-23 CRAN (R 3.2.1)                    
##  digest       0.6.9      2016-01-08 CRAN (R 3.2.1)                    
##  dplyr      * 0.4.3      2015-09-01 CRAN (R 3.2.1)                    
##  evaluate     0.8.3      2016-03-05 CRAN (R 3.2.3)                    
##  formatR      1.3        2016-03-05 CRAN (R 3.2.3)                    
##  htmltools    0.3.5      2016-03-21 CRAN (R 3.2.3)                    
##  knitr        1.12.27    2016-04-30 Github (yihui/knitr@77de0a4)      
##  lazyeval     0.1.10     2015-01-02 CRAN (R 3.1.2)                    
##  magrittr   * 1.5        2014-11-22 CRAN (R 3.1.2)                    
##  memoise      1.0.0      2016-01-29 CRAN (R 3.2.1)                    
##  mrgsolve   * 0.6.1.9000 2016-05-22 local                             
##  R6           2.1.2      2016-01-26 CRAN (R 3.2.3)                    
##  rbenchmark * 1.0.0      2012-08-30 CRAN (R 3.0.3)                    
##  Rcpp         0.12.4     2016-03-26 CRAN (R 3.2.3)                    
##  rmarkdown    0.9.6      2016-04-30 Github (rstudio/rmarkdown@e07c5f6)
##  stringi      1.0-1      2015-10-22 CRAN (R 3.2.1)                    
##  stringr      1.0.0      2015-04-30 CRAN (R 3.1.3)                    
##  yaml         2.1.13     2014-06-12 CRAN (R 3.0.2)

[Eliford]

Thanks @kylebmetrum . I have used approach 3. I wrote a function to sample the constrained etas. But it is messy and may not be reusable if number of etas needing to be constrained change.

#' @param n =  number of subjects
#' @param BSVCL, CLV_COV, BSVVC,CLKA_COV,VKA_COV,BSVKA = values for covariance matrix as in cmat(0.2, 0.7, 0.3, 0, 0, 0.3)
#' @param limeta1 = limit for BSVCL
#' @param limeta2 = limit for BSVVC
#' @param limeta3 = limit for BSVKA

sampling_PKETA_with_constrain<-function(n, BSVCL, CLV_COV, BSVVC,CLKA_COV,VKA_COV,BSVKA,limeta1, limeta2,limeta3){  
  suppressMessages(lapply(c("MASS","dplyr","mrgsolve"), require, character.only=TRUE))
  etas<-lapply(1:n, function(x,BSVCL,CLV_COV,BSVVC,CLKA_COV,VKA_COV,BSVKA,limeta1,limeta2,limeta3){
    myetas<-c(ETA1=1, ETA2=1, ETA3=1)
    while(abs(myetas[["ETA1"]])>limeta1|abs(myetas[["ETA2"]])>limeta2|abs(myetas[["ETA3"]])>limeta3){ 
      myetas <- mvrnorm(1, c(ETA1=0,ETA2=0, ETA3=0),cmat(BSVCL, CLV_COV, BSVVC,CLKA_COV,VKA_COV,BSVKA))

      if(abs(myetas[["ETA1"]])<limeta1&abs(myetas[["ETA2"]])<limeta2&abs(myetas[["ETA3"]])<limeta3){
        break;
      }
    }
    myetas<-data.frame(as.list(myetas))
    return(myetas)
  },  BSVCL,CLV_COV,BSVVC,CLKA_COV,VKA_COV,BSVKA,limeta1,limeta2,limeta3) 
  etas<-bind_rows(etas)
  return(etas)
}  

myetas<-sampling_PKETA_with_constrain(n=1000, BSVCL=0.3, CLV_COV=0.9, BSVVC=0.3,CLKA_COV=0,VKA_COV=0,BSVKA=0.3,limeta1=0.88, limeta2=0.88,limeta3=0.88)
summary(myetas$ETA1)
summary(myetas$ETA2)
summary(myetas$ETA3)

devtools::session_info()

I will continue this approach while waiting for your optimal solution. Thanks

[kylebaron]

Follow development here: https://github.com/metrumresearchgroup/mrgsolve/issues/53