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