# Specifying a mixture model

**URL:** <https://discourse.mc-stan.org/t/specifying-a-mixture-model/9724>\
**Category:** brms\
**Created:** [July 15, 2019, 2:05pm UTC](https://discourse.mc-stan.org/t/specifying-a-mixture-model/9724 "2019-07-15T14:05:17Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![petermacp](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/petermacp/32/3277_2.png) [@petermacp](https://discourse.mc-stan.org/u/petermacp)\
**Post date:** [July 15, 2019, 2:05pm UTC](https://discourse.mc-stan.org/t/specifying-a-mixture-model/9724/1 "2019-07-15T14:05:17Z")

</div>

Please also provide the following information in addition to your question:

- Operating System: OSX 10.14.5
- brms Version: 2.8.0

Hi,

I have data on the measurements of diameter (here in mm) of a skin test that look this this:

```
library(tidyverse)
library(brms)
#> Loading required package: Rcpp
#> Registered S3 method overwritten by 'xts':
#> method from
#> as.zoo.xts zoo
#> Loading 'brms' package (version 2.8.0). Useful instructions
#> can be found by typing help('brms'). A more detailed introduction
#> to the package is available through vignette('brms_overview').

set.seed(1234)

test <- data.frame(result = c(rpois(500, 0.25), rpois(150, 10)))

test %>%
  ggplot() +
  geom_bar(aes(x=result))

```

![](https://i.imgur.com/J9wjbfS.png)

Created on 2019-07-15 by the [reprex package](https://reprex.tidyverse.org) (v0.3.0)

There are two underlying processes to be estimated:

1. Between 0 and around 4, most people probably don’t have the disease, although some will
2. At 5 or above, most people probably do have the disease, although some people with a smaller result may in fact have disease

I have been reading around mixture models as implemented in `brms`, but have been struggling to specify correctly for my data. I have tried the following, with a message about not mixing families:

```
library(brms)
#> Loading required package: Rcpp
#> Registered S3 method overwritten by 'xts':
#> method from
#> as.zoo.xts zoo
#> Loading 'brms' package (version 2.8.0). Useful instructions
#> can be found by typing help('brms'). A more detailed introduction
#> to the package is available through vignette('brms_overview').
mix <- mixture(poisson, gaussian)
#> Error: Cannot mix families with real and integer support.

prior <- c(
  prior(poisson(0.25), Intercept, dpar = mu1),
  prior(gaussian(10, 4), Intercept, dpar = mu2)
)

fit <- brm(bf(result ~ 1), 
            data = test,
            family = mix, 
           prior = prior,
           chains = 1)
#> Error in validate_formula.brmsformula(formula, data = data, family = family, : object 'mix' not found

pp_check(fit)
#> Error in pp_check(fit): object 'fit' not found

```

Created on 2019-07-15 by the [reprex package](https://reprex.tidyverse.org) (v0.3.0)

Thanks.

---

<div class="post-metadata">

**Author:** ![dylanjm](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/dylanjm/32/5018_2.png) [@dylanjm](https://discourse.mc-stan.org/u/dylanjm)\
**Post date:** [July 15, 2019, 4:22pm UTC](https://discourse.mc-stan.org/t/specifying-a-mixture-model/9724/2 "2019-07-15T16:22:39Z")

</div>

I’m not an expert and fairly new to using `library(brms)` myself. Have you tried making a mixture family that’s both Poisson’s? Your example data was generated using 2 Poissons, so since (in this example) you are privy to the underlying process it would make sense to do that.

---

<div class="post-metadata">

**Author:** ![petermacp](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/petermacp/32/3277_2.png) [@petermacp](https://discourse.mc-stan.org/u/petermacp)\
**Post date:** [July 15, 2019, 6:51pm UTC](https://discourse.mc-stan.org/t/specifying-a-mixture-model/9724/3 "2019-07-15T18:51:54Z")

</div>

Thanks so much - yes, I did try that too, but without success:

```
library(brms)
#> Loading required package: Rcpp
#> Registered S3 method overwritten by 'xts':
#> method from
#> as.zoo.xts zoo
#> Loading 'brms' package (version 2.8.0). Useful instructions
#> can be found by typing help('brms'). A more detailed introduction
#> to the package is available through vignette('brms_overview').
mix <- mixture(poisson, poisson)
#> Setting order = 'mu' for mixtures of the same family.

prior <- c(
  prior(poisson(0.25), Intercept, dpar = mu1),
  prior(poisson(10), Intercept, dpar = mu2)
)

fit <- brm(bf(result ~ 1), 
            data = test,
            family = mix, 
           prior = prior,
           chains = 1)
#> Error in isTRUE(attr(data, "brmsframe")): object 'test' not found

pp_check(fit)
#> Error in pp_check(fit): object 'fit' not found

```

Created on 2019-07-15 by the [reprex package](https://reprex.tidyverse.org) (v0.3.0)

Gives the error message:

```
SYNTAX ERROR, MESSAGE(S) FROM PARSER:

No matches for: 

  poisson_lpdf(real, real)

Function poisson_lpdf not found.
  error in 'modelb44f45f46c52_fileb44f6d3f6de' at line 29, column 51
  -------------------------------------------------
    27: vector[N] mu2 = temp_mu2_Intercept + rep_vector(0, N);
    28: // priors including all constants
    29: target += poisson_lpdf(temp_mu1_Intercept | 0.25);
                                                          ^
    30: target += poisson_lpdf(temp_mu2_Intercept | 10);
  -------------------------------------------------

Error in stanc(model_code = paste(program, collapse = "\n"), model_name = model_cppname, : 
  failed to parse Stan model 'fileb44f6d3f6de' due to the above error.
 SYNTAX ERROR, MESSAGE(S) FROM PARSER:

No matches for: 

  poisson_lpdf(real, real)

Function poisson_lpdf not found.
  error in 'modelb44f71642f97_fileb44f6d3f6de' at line 29, column 51
  -------------------------------------------------
    27: vector[N] mu2 = temp_mu2_Intercept + rep_vector(0, N);
    28: // priors including all constants
    29: target += poisson_lpdf(temp_mu1_Intercept | 0.25);
                                                      ^
    30: target += poisson_lpdf(temp_mu2_Intercept | 10);
 -------------------------------------------------

Error in stanc(model_code = paste(program, collapse = "\n"), model_name = model_cppname, : 
  failed to parse Stan model 'fileb44f6d3f6de' due to the above error.
```

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [July 15, 2019, 8:38pm UTC](https://discourse.mc-stan.org/t/specifying-a-mixture-model/9724/4 "2019-07-15T20:38:35Z")

</div>

This looks strange. Do you have a full reproducible example for me?

---

<div class="post-metadata">

**Author:** ![dylanjm](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/dylanjm/32/5018_2.png) [@dylanjm](https://discourse.mc-stan.org/u/dylanjm)\
**Post date:** [July 15, 2019, 9:29pm UTC](https://discourse.mc-stan.org/t/specifying-a-mixture-model/9724/5 "2019-07-15T21:29:28Z")

</div>

Not OP, but I was testing his code and have a reprex demonstrating his exact problem:

```
library(tidyverse)
library(brms)
#> Loading required package: Rcpp
#> Loading 'brms' package (version 2.9.0). Useful instructions
#> can be found by typing help('brms'). A more detailed introduction
#> to the package is available through vignette('brms_overview').

set.seed(1234)

test <- data.frame(result = c(rpois(500, 0.25), rpois(150, 10)))

mix <- mixture(poisson, poisson)
#> Setting order = 'mu' for mixtures of the same family.

prior <- c(
    prior(poisson(0.25), class = "Intercept", dpar = "mu1"),
    prior(poisson(10), class = "Intercept", dpar = "mu2")
)

fit.brm <- brm(bf(result ~ 1), 
           data = test,
           family = mix, 
           prior = prior,
           chains = 1)
#> SYNTAX ERROR, MESSAGE(S) FROM PARSER:
#> 
#> No matches for:
#> 
#> poisson_lpdf(real, real)
#> 
#> Function poisson_lpdf not found.
#> error in 'model744555c541_file74455ad5a546' at line 29, column 52
#> -------------------------------------------------
#> 27: vector[N] mu2 = temp_mu2_Intercept + rep_vector(0, N);
#> 28: // priors including all constants
#> 29: target += poisson_lpdf(temp_mu1_Intercept | 0.25);
#> ^
#> 30: target += poisson_lpdf(temp_mu2_Intercept | 10);
#> -------------------------------------------------
#> 
#> Error in stanc(model_code = paste(program, collapse = "\n"), model_name = model_cppname, : 
#> failed to parse Stan model 'file74455ad5a546' due to the above error.
#> SYNTAX ERROR, MESSAGE(S) FROM PARSER:
#> No matches for:
#> 
#> poisson_lpdf(real, real)
#> 
#> Function poisson_lpdf not found.
#> error in 'model744533d0b466_file74455ad5a546' at line 29, column 52
#> -------------------------------------------------
#> 27: vector[N] mu2 = temp_mu2_Intercept + rep_vector(0, N);
#> 28: // priors including all constants
#> 29: target += poisson_lpdf(temp_mu1_Intercept | 0.25);
#> ^
#> 30: target += poisson_lpdf(temp_mu2_Intercept | 10);
#> -------------------------------------------------
#> 
#> Error in stanc(model_code = paste(program, collapse = "\n"), model_name = model_cppname, : failed to parse Stan model 'file74455ad5a546' due to the above error.

```

Created on 2019-07-15 by the [reprex package](https://reprex.tidyverse.org) (v0.3.0)

---

<div class="post-metadata">

**Author:** ![MauritsM](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mauritsm/32/430_2.png) [@MauritsM](https://discourse.mc-stan.org/u/MauritsM)\
**Post date:** [July 15, 2019, 10:18pm UTC](https://discourse.mc-stan.org/t/specifying-a-mixture-model/9724/6 "2019-07-15T22:18:38Z")

</div>

At the risk of stating the obvious, the Stan code should probably be `poisson_lpmf` instead of `_lpdf`. Hope this helps!

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [July 15, 2019, 10:34pm UTC](https://discourse.mc-stan.org/t/specifying-a-mixture-model/9724/7 "2019-07-15T22:34:22Z")

</div>

Ah, now that I look at this code, indeed, why did you specify a poison prior for the intercepts? This can’t possibly work as parameters in stan are continuous while the poisson is a discrete distribution.

---

<div class="post-metadata">

**Author:** ![petermacp](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/petermacp/32/3277_2.png) [@petermacp](https://discourse.mc-stan.org/u/petermacp)\
**Post date:** [July 16, 2019, 7:27am UTC](https://discourse.mc-stan.org/t/specifying-a-mixture-model/9724/8 "2019-07-16T07:27:10Z")

</div>

Okay - that makes sense. I guess I am still struggling to conceptualise how to specify priors to reflect my understanding of the data here. I think my difficultly is in ensuring that the first distribution doesn’t fall below zero, which wouldn’t be possible in the data:

```
library(brms)
#> Loading required package: Rcpp
#> Registered S3 method overwritten by 'xts':
#> method from
#> as.zoo.xts zoo
#> Loading 'brms' package (version 2.8.0). Useful instructions
#> can be found by typing help('brms'). A more detailed introduction
#> to the package is available through vignette('brms_overview').
test <- data.frame(result = c(rpois(500, 0.25), rpois(150, 10)))

mix <- mixture(gaussian, gaussian)
#> Setting order = 'mu' for mixtures of the same family.

prior <- c(
  prior(normal(0, 2), Intercept, dpar = mu1),
  prior(normal(0, 2), Intercept, dpar = mu2)
)

fit <- brm(bf(result ~ 1), 
            data = test,
            family = mix, 
           prior = prior,
           chains = 1)
#> Compiling the C++ model
#> Start sampling
#> 
#> SAMPLING FOR MODEL '7b517954ea2a96396e16f5dec78c77c2' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.000286 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.86 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration: 1 / 2000 [0%] (Warmup)
#> Chain 1: Iteration: 200 / 2000 [10%] (Warmup)
#> Chain 1: Iteration: 400 / 2000 [20%] (Warmup)
#> Chain 1: Iteration: 600 / 2000 [30%] (Warmup)
#> Chain 1: Iteration: 800 / 2000 [40%] (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [50%] (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [50%] (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [60%] (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [70%] (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [80%] (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [90%] (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 1: 
#> Chain 1: Elapsed Time: 1.04822 seconds (Warm-up)
#> Chain 1: 1.2112 seconds (Sampling)
#> Chain 1: 2.25942 seconds (Total)
#> Chain 1:

summary(fit)
#> Family: mixture(gaussian, gaussian) 
#> Links: mu1 = identity; sigma1 = identity; mu2 = identity; sigma2 = identity; theta1 = identity; theta2 = identity 
#> Formula: result ~ 1 
#> Data: test (Number of observations: 650) 
#> Samples: 1 chains, each with iter = 2000; warmup = 1000; thin = 1;
#> total post-warmup samples = 1000
#> 
#> Population-Level Effects: 
#> Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
#> mu1_Intercept 0.20 0.02 0.17 0.24 1479 1.00
#> mu2_Intercept 9.04 0.30 8.43 9.62 1197 1.00
#> 
#> Family Specific Parameters: 
#> Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
#> sigma1 0.40 0.01 0.38 0.43 1241 1.00
#> sigma2 3.56 0.22 3.17 4.02 958 1.00
#> theta1 0.74 0.02 0.71 0.78 1201 1.00
#> theta2 0.26 0.02 0.22 0.29 1201 1.00
#> 
#> Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
#> is a crude measure of effective sample size, and Rhat is the potential 
#> scale reduction factor on split chains (at convergence, Rhat = 1).
pp_check(fit)
#> Using 10 posterior samples for ppc type 'dens_overlay' by default.

```

![](https://i.imgur.com/MtcwwtQ.png)

Created on 2019-07-16 by the [reprex package](https://reprex.tidyverse.org) (v0.3.0)

I am struggling to work out how to specify this prior, and grateful for any suggestions.

Thanks

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [July 16, 2019, 8:41am UTC](https://discourse.mc-stan.org/t/specifying-a-mixture-model/9724/9 "2019-07-16T08:41:27Z")

</div>

What you are looking for is **not** a prior but a likelihood (specified via the `family` argument). In this case, a mixture of two poisson distributions. Use `mixture(poisson, poisson)` instead of `mixture(gaussian, gaussian)`.

---

<div class="post-metadata">

**Author:** ![petermacp](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/petermacp/32/3277_2.png) [@petermacp](https://discourse.mc-stan.org/u/petermacp)\
**Post date:** [July 16, 2019, 8:46am UTC](https://discourse.mc-stan.org/t/specifying-a-mixture-model/9724/10 "2019-07-16T08:46:03Z")

</div>

Thanks - I have tried that, but get an error message:

`Error: Family 'mixture(poisson, poisson)' requires response greater than or equal to 0.`

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [July 16, 2019, 9:46am UTC](https://discourse.mc-stan.org/t/specifying-a-mixture-model/9724/11 "2019-07-16T09:46:43Z")

</div>

But didn’t you create data according to a poisson distribution or is your actual data structured differently?

---

<div class="post-metadata">

**Author:** ![petermacp](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/petermacp/32/3277_2.png) [@petermacp](https://discourse.mc-stan.org/u/petermacp)\
**Post date:** [July 16, 2019, 9:54am UTC](https://discourse.mc-stan.org/t/specifying-a-mixture-model/9724/12 "2019-07-16T09:54:56Z")

</div>

You are correct. In trying to get working, I had changed my test data to see if that works. With this specified as:

```
library(tidyverse)
library(brms)
#> Loading required package: Rcpp
#> Registered S3 method overwritten by 'xts':
#> method from
#> as.zoo.xts zoo
#> Loading 'brms' package (version 2.8.0). Useful instructions
#> can be found by typing help('brms'). A more detailed introduction
#> to the package is available through vignette('brms_overview').

set.seed(1234)

test <- data.frame(result = c(rpois(500, 0.25), rpois(150, 10)))

test %>%
  ggplot() +
  geom_bar(aes(x=result))

```

![](https://i.imgur.com/GBpzI4J.png)

```

mix <- mixture(poisson, poisson)
#> Setting order = 'mu' for mixtures of the same family.

prior <- c(
  prior(cauchy(0,5), Intercept, dpar = mu1),
  prior(normal(0,2), Intercept, dpar = mu2)
)

fit <- brm(bf(result ~ 1), 
            data = test,
            family = mix, 
           prior = prior,
           chains = 1)
#> Compiling the C++ model
#> Start sampling
#> 
#> SAMPLING FOR MODEL '46b36408f139ed6444538aceb75a3d47' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.000296 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.96 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration: 1 / 2000 [0%] (Warmup)
#> Chain 1: Iteration: 200 / 2000 [10%] (Warmup)
#> Chain 1: Iteration: 400 / 2000 [20%] (Warmup)
#> Chain 1: Iteration: 600 / 2000 [30%] (Warmup)
#> Chain 1: Iteration: 800 / 2000 [40%] (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [50%] (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [50%] (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [60%] (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [70%] (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [80%] (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [90%] (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 1: 
#> Chain 1: Elapsed Time: 4.23044 seconds (Warm-up)
#> Chain 1: 3.74466 seconds (Sampling)
#> Chain 1: 7.9751 seconds (Total)
#> Chain 1:

summary(fit)
#> Family: mixture(poisson, poisson) 
#> Links: mu1 = log; mu2 = log; theta1 = identity; theta2 = identity 
#> Formula: result ~ 1 
#> Data: test (Number of observations: 650) 
#> Samples: 1 chains, each with iter = 2000; warmup = 1000; thin = 1;
#> total post-warmup samples = 1000
#> 
#> Population-Level Effects: 
#> Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
#> mu1_Intercept -1.47 0.10 -1.66 -1.29 286 1.00
#> mu2_Intercept 2.26 0.03 2.21 2.31 1060 1.00
#> 
#> Family Specific Parameters: 
#> Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
#> theta1 0.77 0.02 0.73 0.80 441 1.00
#> theta2 0.23 0.02 0.20 0.27 441 1.00
#> 
#> Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
#> is a crude measure of effective sample size, and Rhat is the potential 
#> scale reduction factor on split chains (at convergence, Rhat = 1).
pp_check(fit)
#> Using 10 posterior samples for ppc type 'dens_overlay' by default.

```

![](https://i.imgur.com/n49n79E.png)

Created on 2019-07-16 by the [reprex package](https://reprex.tidyverse.org) (v0.3.0)

This now works. Thanks so much for your help - has really helped my understanding here. Much appreciated.