Sharing parameters across distributions in a mixture model

Interfaces brms
,
1

Hi,
I was wondering if I could share parameters across the distributions of a mixture in brms. Imagine a situation where a school grade depends on a predictor, but a school teacher might have favorite students who they don’t want to fail and so another distribution (of favorites) is truncated at 59.5:

library(extraDistr)
library(brms)
set.seed(123)
# number of students
N <- 400
alpha <- 55
beta <- 2
sigma <- 8
predictor <- rnorm(N)
# scores for "fair" grades
scores_1 <- rtnorm(N, alpha + beta * predictor,sigma, a =0,b=100)
# scores for grades that are always a pass
scores_2 <- rtnorm(N, alpha + beta * predictor,sigma, a =59.5,b=100)
# prob that the student is not a favorite:
p <- .7
# which students are a favorite
z <- rbern(N, p)
# distribution of scores:
scores <- ifelse(z, scores_1, scores_2)
hist(scores)

Created on 2023-05-03 with reprex v2.0.2

Now how do I code this in brms? I tried the following:

df <- data.frame(scores, predictor)

brm(bf(scores ~ predictor ,
       mu1 | rep_trunc(lb=0,ub=100)~ 0,
       mu2 | rep_trunc(lb=59.5,ub=100)~ 0
       )  , data= df,
    family = mixture(gaussian, gaussian, order = 'none'),
prior = ...
)

But I can see with make_stancode that the model is ignoring predictor. But if I move it to the lines of mu1 and mu2:

mu1 | rep_trunc(lb=0,ub=100)~ predictor,
mu2 | rep_trunc(lb=59.5,ub=100)~ predictor

I end up with two betas. (Ideally I would like to share sigma across models, but it’s less critical).

(I know how to do it in Stan, and I know how to edit the brms model, but I wonder if there’s a brms solution that I can use for teaching. And by the way this is a very realistic situation, most datasets with grades have a bump around 60).

  • Operating System: Ubuntu 20.04.5 LTS
  • brms Version: brms_2.19.0
2

I’m not totally sure I understand the model you’re trying to fit, or what sort of priors might be reasonable here, but I think you can do this in brms by using a non-linear model. For example:

mix <- mixture(gaussian, nmix=2, order = "none")
formula <- bf(scores | resp_trunc(lb=0, ub=100) ~ alpha + beta*predictor,
              nlf(mu2 | resp_trunc(lb=59.5, ub=100) ~ alpha + beta*predictor),
              nlf(sigma1 ~ x),
              nlf(sigma2 ~ x),
              alpha ~ 1,
              beta ~ 1,
              x ~ 1,
              nl = TRUE,
              family=mix)
get_prior(formula, data=df)
priors = c(set_prior("normal(50, 20)", nlpar="alpha", class="b"),
           set_prior("normal(0, 5)", nlpar="beta", class="b"),
           set_prior("normal(0, 10)", nlpar="x", class="b"))
fit <- brm(formula, data=df, prior = priors)

This seems to recover the parameter values okay, though not the mixing proportion.

1 Like