$MIX option in mrgsolve?

[JasonW000]

Is there an option in mrgsolve similar to $MIX in NONMEM, with the goal of simulating from different populations? My current workaround is simply to run two separate simulations with the number of individuals based on the probability of being in population 1 vs 2 and then merging them into a single dataframe. This works, but if the option exists I would like to do this within the mrgsolve model.

thanks

Jason

[kylebaron]

I'm familiar with $MIX, but I'm trying to figure out what exactly needs to be done here.

Do you want mrgsolve to "randomize" people to this population or that population as the run proceeds?

[JasonW000]

Yes, randomization based on a probability that is defined as a model parameter rather than being defined outside of the model. Since it is all in R it is not to difficult to make it work using the latter approach. If this option is not currently available that is okay, but I wanted to check.

Jason

Sent from my iPhone

On May 25, 2017, at 3:04 PM, Kyle Baron notifications@github.com wrote:

I'm familiar with $MIX, but I'm trying to figure out what exactly needs to be done here.

Do you want mrgsolve to "randomize" people to this population or that population as the run proceeds?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub, or mute the thread.

[kylebaron]

There isn't a pre-wrapped $MIX per se. But it's not too hard to code your own.

I put this in $MAIN ... but I might think about wrapping it up in a function for heavy project work.

For two groups, just call rbinom. For three groups you can do it with uniform distribution? Or maybe there's a more-clever way to do it.

If you have more groups and there's an R function that will so this really neatly, I can help you implement that as well.

``` r
library(mrgsolve)

code <- '
$PLUGIN Rcpp 

$MAIN
if(NEWIND <=1) {
  int mix1 = R::rbinom(1,0.633);
  double mixu = R::runif(0,1);
  int mix2 = mixu < 0.33;
  int mix3 = mixu > 0.33 && mixu < 0.88;
  int mix4 = mixu > 0.88;
}
$CAPTURE mix1 mix2 mix3 mix4
'

mod <- mcode("foo", code)

## Compiling foo ...

## done.

out <- mrgsim(mod, nid=100000, end=-1)

head(out)

## Model:  foo

##   ID time mix1 mix2 mix3 mix4
## 1  1    0    0    0    1    0
## 2  2    0    1    0    1    0
## 3  3    0    1    1    0    0
## 4  4    0    1    1    0    0
## 5  5    0    1    0    0    1
## 6  6    0    0    1    0    0

mean(out$mix1)

## [1] 0.63241

mean(out$mix2)

## [1] 0.32859

mean(out$mix3)

## [1] 0.55308

mean(out$mix4)

## [1] 0.11833

[kylebaron]

For the previous example, I hard-coded the probabilities. You can make those "parameters" or other variables. Doesn't matter. Only call R::rbinom with 1 as the first argument and obviously use runif(0,1) for that approach.

[dpastoor]

well kyle beat me to it, but here was the example I cooked up as well - in this case adjusting CL by 50% for one group (CLGRP = 0.5) with a probability of 30% to be in the GRP

library(mrgsolve)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

mix_model <- "
$PLUGIN Rcpp

$PARAM CL = 1, CLGRP = 0.5, VC = 30, KA = 1.3

$CMT GUT CENT

$MAIN

if(NEWIND <=1) {
double GRP = R::rbinom(1,0.3);
}

$ODE
double CLi = CL*(1-GRP*CLGRP);
dxdt_GUT = -KA*GUT;
dxdt_CENT = KA*GUT - (CLi/VC)*CENT;

$TABLE
double CP = CENT/VC;

$CAPTURE CP GRP CL CLi
"

mod <- mcode_cache("mixture", mix_model)
#> Compiling mixture ...
#> done.

raw_ids <- mod %>% mrgsim(nid = 100) %>% as_data_frame %>% distinct(ID, .keep_all = T)

head(raw_ids, n = 10)
#> # A tibble: 10 × 8
#>       ID  time   GUT  CENT    CP   GRP    CL   CLi
#>    <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#>  1     1     0     0     0     0     0     1   1.0
#>  2     2     0     0     0     0     1     1   0.5
#>  3     3     0     0     0     0     1     1   0.5
#>  4     4     0     0     0     0     0     1   1.0
#>  5     5     0     0     0     0     1     1   0.5
#>  6     6     0     0     0     0     0     1   1.0
#>  7     7     0     0     0     0     0     1   1.0
#>  8     8     0     0     0     0     0     1   1.0
#>  9     9     0     0     0     0     0     1   1.0
#> 10    10     0     0     0     0     0     1   1.0

raw_ids %>%
  count(CLi)
#> # A tibble: 2 × 2
#>     CLi     n
#>   <dbl> <int>
#> 1   0.5    25
#> 2   1.0    75

[kylebaron]

@dpastoor Great minds think alike.

But somehow I get the feeling that I'm going to be reporting to you someday.

[dpastoor]

ha! you do all the hard work - I just get to have fun playing with it :-)

[kylebaron]

Captured some of this in a blog post: http://mrgsolve.github.io/2017/05/26/mixture-model-example/

[JasonW000]

Thank you both for the useful examples and the thorough explanation in the blog post. I have implemented this in my model and it is working out very well.

On Fri, May 26, 2017 at 10:39 PM, Kyle Baron notifications@github.com wrote:

Captured some of this in a blog post: http://mrgsolve.github.io/2017/05/26/mixture-model-example/

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/metrumresearchgroup/mrgsolve/issues/236#issuecomment-304429610, or mute the thread https://github.com/notifications/unsubscribe-auth/AUSKSW5v8r0LxSFrU1Ex804786_TjAUtks5r97cRgaJpZM4Nm8W4 .

[kylebaron]

Thanks for the good question, Jason.

[JasonW000]

Kyle, in your example (shown below), is there a reason you implement the ETA within the $MAIN block instead of using a standard ETA and $OMEGA ? It appears I am able to get satisfactory results using either approach. However, I like using the latter approach since I can easily add Covariance parameters using the 'cmat' function. Please let me know if I am overlooking something here.

thanks Jason

code <- ' $PLUGIN Rcpp $PARAM TVCL = 1, V = 30, KA=1.2, THETA1 = 0.5 $PKMODEL cmt="GUT CENT", depot=TRUE $MAIN if(NEWIND <=1) { int POP = 1 + R::rbinom(1,0.2); double myETA = R::rnorm(0,sqrt(0.09)); } double CL = TVCL; if(POP==2) CL = TVCL THETA1; double CLi = CLexp(myETA);

$CAPTURE POP CL CLi '

mod <- mcode_cache("C", code)

On Fri, May 26, 2017 at 10:39 PM, Kyle Baron notifications@github.com wrote:

Captured some of this in a blog post: http://mrgsolve.github.io/2017/05/26/mixture-model-example/

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/metrumresearchgroup/mrgsolve/issues/236#issuecomment-304429610, or mute the thread https://github.com/notifications/unsubscribe-auth/AUSKSW5v8r0LxSFrU1Ex804786_TjAUtks5r97cRgaJpZM4Nm8W4 .

[kylebaron]

Hi @JasonW000 -

I only did it here to illustrate another way to simulate random numbers. I think we use R::runif too later on in the example (or maybe that was on the blog post?) ... that's a pretty decent starter set.

But I agree with you that it's better / easier to stick with using $OMEGA and let mrgsolve simulate ETAs for you. I would definitely use $OMEGA for real work.

Kyle

[JasonW000]

Perfect, thank you for the clarification.

On Fri, Jun 23, 2017 at 11:10 AM, Kyle Baron notifications@github.com wrote:

Hi @JasonW000 https://github.com/jasonw000 -

I only did it here to illustrate another way to simulate random numbers. I think we use R::runif too later on in the example (or maybe that was on the blog post?) ... that's a pretty decent starter set.

But I agree with you that it's better / easier to stick with using $OMEGA and let mrgsolve simulate ETAs for you. I would definitely use $OMEGA for real work.

Kyle

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/metrumresearchgroup/mrgsolve/issues/236#issuecomment-310735528, or mute the thread https://github.com/notifications/unsubscribe-auth/AUSKSUm6eewiqD8gVnEL22Id0GINzRbMks5sG_-JgaJpZM4Nm8W4 .