IRT Mixture Bernoulli

I’m using the following model with one class. I would like to integrate mixture parameters to isolate specific groups.

Would you have an idea of how to change the code below so that it includes groups parameters?

Similar problematic is available here: Mixture Bernoulli

model <- "
  data {
      int<lower=0> n; // number of subjects
      int<lower=0> k; // number of items
      int<lower=0, upper=1> y1[n]; // manifestation variables
      int<lower=0, upper=1> y2[n];
      int<lower=0, upper=1> y3[n];
      int<lower=0, upper=1> y4[n];
      int<lower=0, upper=1> y5[n];
      int<lower=0, upper=1> y6[n];
      int<lower=0, upper=1> y7[n];
      int<lower=0, upper=1> y8[n];
  parameters {
      vector[k] alpha;
      real<lower=0> beta[k];
      vector[n] x;    // this is theta but its easier to see as x in the code
  model {
      // priors
      x ~ normal(0,1); 
      alpha ~ normal(0,4); 
      beta ~ gamma(4,3); 
      // items
  y1 ~ bernoulli_logit(alpha[1] + beta[1] * x);
  y2 ~ bernoulli_logit(alpha[2] + beta[2] * x);
  y3 ~ bernoulli_logit(alpha[3] + beta[3] * x);
  y4 ~ bernoulli_logit(alpha[4] + beta[4] * x);
  y5 ~ bernoulli_logit(alpha[5] + beta[5] * x);
  y6 ~ bernoulli_logit(alpha[6] + beta[6] * x);
  y7 ~ bernoulli_logit(alpha[7] + beta[7] * x);
  y8 ~ bernoulli_logit(alpha[8] + beta[8] * x);

Is to goal to model alpha and beta hierarchically? Is it you have these k things and you want to do partial pooling?

Or is it that you have g groups, and that’s a different thing, and you’d like to add it to the model? Like, would your desired model look like vector[k] alpha[g];?

Hello @bbbales2 ,

Yes, it’s more links to the second proposition.
So fare I identify n parameters alpha[n] and beta[n] like that: alpha[1] , beta[1]…
And I would like to integrate mixture parameter wich “g” groups to estimate alpha[n,g] and beta[n,g] like: alpha[item 1, group 1], alpha[item 1,group 2], beta[item 1,group 1], beta[item 1,group 2].

As you present, I think that I need to change both parameters and also loop around the likelihood but I’m confused about the implementation…

Okay so it sounds like you have N bernoulli observations of K items and you want to model the outcomes in each case as a mixture?

I think in a situation like this you just gotta futz with the math a while and generate data and fit models to double check yourself.

So if we’re just talking about one outcome, the density over our outputs and our unobserved mixture components looks like:

p(y_{nk}, g_k) = p(y_{nk} | g_k) p(g_k)

The n index is the observations index, the k index is the items index, and the g index is the mixture component index.

So we don’t know g_k so we integrate it out:

\sum_{g_k} p(y_{nk}, g_k) = p(y_{nk} | g_k = 1) p(g_k = 1) + p(y_{nk} | g_k = 2) p(g_k = 2)

And so if \lambda = p(g_k = 1), then I think the likelihood terms in @Random 's first post is correct (though I could easily be wrong), something like:

target += log_mix(lambda, bernoulli_logit_lpmf(y[n, k] | p[k, 1]), bernoulli_logit_lpmf(y[n, k] | p[k, 2]));

Also, define and pass in y like a two dimension array:

int<lower=0, upper=1> y[n, k];

I think it’ll make the code easier to work with.