- Operating System: Win 10.0.17134
- brms Version: 2.9.0
Greetings, I’m looking for advice on how to proceed with a current analytical challenge within the brms framework…
Background: analysing data from a single survey given to different sampling cohorts (
sample, n = 8) in which all participants (
subject, n = 500) answered the same question set under two fictional scenarios (
subject is nested within
scenario is fully-crossed within
subject. The question set contained 32 statements, each rated on the same 7pt ordinal scale (1-7, strongly disagree to strongly agree, although only verbal labels were displayed to participants). Most of the statements are indicators of five outcomes scales of interest (number of indicators ranges from 3-6 for a given scale), and the remaining nine statements comprise another “scale” that, for simplicity, is being scored via averaging a person’s responses across the items and treated as a standardised (z-score) continuous predictor (
pred_1). Of central interest is the effect of
pred_1 on a given outcome, and whether or not the effect of
pred_1 varies across
scenario. All outcomes are correlated with each other.
Current Approach: So far, and following Bürkner & Vuorre (2019), I’ve used separate multilevel ordinal regression models to predict scores for a each outcome (
resp ~ scenario * pred_1 + (1|item) + (1|subject) + (1|sample),
family = cumulative("logit")
prior <- prior("normal(0, 3)", class = "sd", group = "sample") + prior("normal(0, 3)", class = "sd", group = "subject") + prior("normal(0, 3)", class = "sd", group = "item")
While not ideal, since outcomes are correlated, this may ultimately suffice. However, I have two questions about such an approach. First, does the model syntax look correct? And second, is it possible to produce a scatterplot for a given model showing:
- predicted values of the latent outcome variable on the y-axis,
- predicted values of the continuous predictor,
pred_1, on the x-axis,
- colour the points by
- overlay lines representing the predicted interaction of
marginal_effects() shows the probability of responses on the observed metric (1-7 scale), which isn’t quite what I’m looking for.
Alternative Approach: Since the outcomes are correlated, might it make more sense to use a custom beta-binomial distribution to estimate all outcomes, as sum-scores, simultaneously? An example of this approach can be found here, although I’ve having difficulty mapping my situation to the model syntax.
Any insight, and/or additional approaches, would be greatly appreciated!