Using a single mixture weight parameter in a multivariate mixture model

Unfortunately, this is not supported and might never be, because it introduces some difficult issues - I discussed this a bit more broadly at Specifying joint mixture model

Your actual use case, where you know that both components are larger at the same time actually avoids the particular problems (and even better if you actually know exactly the amount of increase) but there is no straightforward way to tell this to brms. I think that you should be able to achieve this particular case with the custom_family feature.

A sketch of the code (I didn’t test it, please double check) - it is a bit hacky, so feel free to ask for clarifications if it is hard to follow. The biggest hack is that we treat this as a univariate custom distribution and abuse the vreal addition term to let us specify the second response. Second, we enforce global ordering by modelling the mixture as a mean of the first component and the positive shift for the second component. For simplicity we assume the same sigma for all components, but the extension to varying should be straightforward.

my_mix <- custom_family(
  "my_mix", dpars = c("mu1", "shift1", "mu2", "shift2", "sigma", "theta"),
  links = c("identity", "log", "identity", "log", "log", "logit"), 
  lb = c(NA, 0, NA, 0, 0, 0),
  ub = c(NA, NA, NA, NA, NA, 1),
  type = "real", vars = "vreal1[n]"
)

stan_funs <- "
  real my_mix_lpdf(real y1, real mu1, real shift1, real mu2, real shift2, real sigma, real theta, real y2) {
     real lp1 = log_mix(theta, normal_lpdf(y1 | mu1, sigma), normal_lpdf(y1 | mu1 + shift1, sigma));
     real lp2 = log_mix(theta, normal_lpdf(y2 | mu2, sigma), normal_lpdf(y2 | mu2 + shift2, sigma));
     return lp1 + lp2;
  }
"

stanvars <- stanvar(scode = stan_funs, block = "functions")

fit <- brm(
  y1 | vreal(y2) ~ 1, data = dat,
  family = my_mix, stanvars = stanvars
)

Best of luck with your model and feel free to ask for clarifications (this stuff can be quite intimidating)

1 Like