I’m new both to this forum and Bayesianism. And a humanities student.
I’ve used only non-Bayesian GL(M)Ms before and have now had to venture into Bayesian territory, after failing to find R packages capable of fitting frequentist multinomial GLMMs. After I fit my first Bayesian model with brm(), several problems arose immediately:
- I cannot find Deviance reported anywhere. I do see a function named log_lik(), but it returns a huge matrix rather than one statistic. Since its dimensions are 4000-by-N, I assume they are the observation-specific, “point-wise” log-likelihoods calculated for each MC sample. If so, does this mean Deviance is simply:
sum(colMeans(log_lik(fit)))  4172.319
If not, then how exactly does one compute good ole Deviance for this kind of model?
- I can’t access the fitted values. Both fitted(fit) and predict(fit, newdata = NULL, re_formula = NULL) report:
Error in Summary.factor(c(2L, 2L, 2L, 3L, 2L, 3L, 2L, 2L, 3L, 4L, 4L, : ‘max’ not meaningful for factors
- I don’t know how to access residuals of any kind – neither Pearson nor Deviance ones, neither regular nor standardized. And I can’t even begin to calculate them on my own without access to the fitted values.
- Trying to obtain the model’s WAIC() results in a warning:
> waic(fit) Computed from 4000 by 2321 log-likelihood matrix Estimate SE elpd_waic -2144.1 38.6 p_waic 109.5 4.1 waic 4288.1 77.2 Warning message: 32 (1.4%) p_waic estimates greater than 0.4. We recommend trying loo instead.
Can this test statistic still be used despite the warning? When I do what it says and try loo() instead, I get an error instead of a warning:
> loo(fit, fit2, compare = TRUE, reloo = TRUE) Error: At least three response categories are required.
What’s going on? There are 4 response categories.
Lastly, there are three priors whose purpose I don’t understand:
> prior_summary(fit) prior class coef group resp dpar nlpar bound 1 normal(0, 2.5) b muindicative 2 b actualized1_actualized muindicative 3 b adjPoly muindicative ... 58 Intercept muindicative 59 Intercept muother 60 Intercept mushould 61 student_t(3, 0, 10) sd muindicative 62 student_t(3, 0, 10) sd muother 63 student_t(3, 0, 10) sd mushould 64 cauchy(0, 3) sd trigger muindicative 65 sd Intercept trigger muindicative 66 cauchy(0, 3) sd trigger muother 67 sd Intercept trigger muother 68 cauchy(0, 3) sd trigger mushould 69 sd Intercept trigger mushould
What are these student priors and what do they do? I deliberately set a fairly skeptical prior on the betas (the normal one) and group-level effects (cauchy), and I deliberately refrained from setting any priors on the fixed intercepts. Yet there appear to be default student priors on the fixed intercepts? It doesn’t say “Intercept” in the class column, but the dpar column leads me to suspect that these priors have something to do with the intercepts because the three non-baseline response categories are named. The helpfile for set_prior() doesn’t seem to say anything about ‘class sd’ priors for fixed intercepts, nor have I ever heard of such a thing before, so I’m confused.
Here’s the code I used to fit the model:
priors <- prior( normal(0, 2.5), class = b) + #prior for betas prior(cauchy(0, 3), class = sd, group = "trigger", dpar = "muindicative") + #group-level prior, response category 2 prior(cauchy(0, 3), class = sd, group = "trigger", dpar = "mushould") + #group-level prior, response category 3 prior(cauchy(0, 3), class = sd, group = "trigger", dpar = "muother") #group-level prior, response category 4 fit <- brm(dep.var ~ (1|trigger) + syn.type + sem.stat + tense + ISCmod + should.merge + sub.verb + time.period + polarity + broadsheet2 + that + actualized + verbPoly + adjPoly, data = testdata, family = categorical(link = "logit"), prior = priors)