Hello all,

I’m trying to fit a distributional model in brms. To start, I just want to fit a model where the mean varies by group. So I simulated some simple data using the following code with fixed values:

Int <- c(-0.5, 0, 0.5)
B <- 0.01
N <- 75
Y <- rnorm(N, mean=0, sd=0.5)
X <- runif(N, min=4, max=6.5)

Res <- 0.01
Y <- c()
CG <- c()
for (i in 1:N) {
CG[i] <- sample(Int, 1)
Y[i] <- B * X[i] + CG[i] + rnorm(1, 0, Res)
}

testdf <- data.frame(Y = Y, X = X, CG = CG)

Then I am trying to fit a simple model as shown here:

bprior <- prior(student_t(5, 0, 5), class = b) +
prior(normal(0, 0.5), class = Intercept) +
prior(cauchy(0, 2.5), class = sd)

fit <- brm( Y ~ X + ( 1 | CG ), data = testdf, prior = bprior, family = gaussian(), control = list(adapt_delta = 0.85, max_treedepth = 10) )

I thought this model matched the generating simulation, but I am getting persistent warnings about post-warmup divergences (no matter how I tweak the adapt_delta). Using brms’ default priors also yields the persistent divergences. Thus I assume the model I specified does not match the simulation above, and so I was hoping for some advice on what’s awry.

Thank you,

–Jon

With only three groups, I can only think of two reasonable models for your data:

`fit <- brm( Y ~ X + CG, data = testdf, prior = bprior, family = gaussian(), control = list(adapt_delta = 0.85, max_treedepth = 10) )`

or

`fit <- brm( Y ~ X*CG, data = testdf, prior = bprior, family = gaussian(), control = list(adapt_delta = 0.85, max_treedepth = 10) )`

I played around with you model a little bit. I’m no expert. But you may have better luck having the model estimate sigma as well.

I don’t understand the variables you’ve created or how you want them to related to Y. the phrasing of your question would imply that you want to actually model Y ~ (1 | CG). (mean changes by group) I also did Y~ (X | CG ) but it was unclear in what way you want X to be present as a predictor. The latter formula represents the effect of X on Y varies by group CG

I tightened up the priors because your data varies very little. The estimated sigma is 0.01. The model predicts the data well [via: pp_check(fit, nsamples=25) ] but it looks super weird to me. The groups being defined as numbers may relate to this, but I don’t know because I deal almost entirely with factor data so this not something I usually see.

As far as too few levels in your grouping variable, one benefit to HMC is that having fewer groups can still be done and have models converge, whereas convergence in lme4 is difficult with fewer than 5 (maybe 4, maybe 3 if you have a LOT of samples)… at least it is a problem when you are trying to model your grouping factors as varying slopes and intercepts. (though if you dig DEEP in to Doug Bates talks and mailing list posts you will find a good approximation for this without breaking your model)

Don’t know if that helped in anyway.

Good luck. I should have been doing my own stuff! oops

4 Likes

In my experience, divergent transitions are common with so few groups. Consider tightening up the `class = sd` prior to something like `cauchy(0, 1)`, `normal(0, 1)`, or even tighter.

3 Likes