Unequal Variance for Multivariate Multilevel Ordinal Models (Generalization)

Modeling
, ,
1

Hello Dear Expert Community

I hope my message finds you well

I am using the brms package to conduct a Bayesian multivariate multilevel model for my PhD thesis, a cross-national study in 7 countries. My aim is to discover how subjective well-being can affect the modes of unconventional political participation. In addition, I use the country’s score of voice and accountability as a moderator.

IVs : Life satisfaction (5 ordinal levels), financial satisfaction (5 ordinal levels), health (4 ordinal levels), happiness(4 ordinal levels), free choice (5 ordinal levels).

DVs: signing petition, joining peaceful demonstration, joining strike, joining boycott. all are 3 ordinal levels: (Would Never Do, Might do, Have Done)

Control Vs: Employed, Retired, Student, Unemployed, Male, careful_trust, age. (all are dummies, except age: interval)

Country-level Variable: country’s score of voice and accountability. (labeled as var)

Data: World Value Survey, I merged the last 2 waves.

Number of Observations is 13000.

Priors: non-informative priors, using the default settings.

The reason for using Bayesian approach is the small number of countries.

…

In order to get most reliable and certain results, I followed Prof. @paul.buerkner guidance in writing the code and model generalization
in his works:

  • Ordinal Regression Models in Psychology: A Tutorial
  • brms: An R Package for Bayesian Multilevel Models Using Stan

generalization for ordinal models according to these two works is done by fitting:

1- Adjacent category logit model.
2- Adjacent category logit specific effects model.
3- Cumulative logit model.
4- Cumulative logit with unequal variances model.

I have two issues relevant to this generalization, and I deeply seek help:

1- Cumulative model with unequal variances works with a single outcome variable when I add the threshold = “equidistant”, but this did not work with multivariate model, is there any special syntax for the multivariate model.

Here is my (not-working model):

library(brms)
fit_happiness2 ← brm(
mvbind(sign_pet, join_boycott, attend_demo, join_strike) ~
happiness * var,
data = scaled,
family = cumulative(link = “logit”, threshold = “equidistant”),
chains = 4,
cores = 4,
iter = 2000,
control = list(adapt_delta = 0.99, max_treedepth = 15)
)

Print the summary

summary(fit_happiness2)

2- for my multilevel multivariate model, do I have to use cse() or cs() for specific effects model.

library(brms)

fit_happiness3 ← brm(
mvbind(sign_pet, join_boycott, attend_demo, join_strike) ~
cs(happiness * var) +
(1 || country),
family = acat(“logit”),
data = scaled,
chains = 4,
cores = 4,
iter = 2000,
control = list(adapt_delta = 0.99, max_treedepth = 12)
)

Print the summary

summary(fit_happiness3)

I have used cs() function, it works pretty well, but I have a concern that I may should have used the cse() instead

Respectfully,
Mohammed Darkhalil

2

cse() is just an outdated name for cs(). So you don’t need to worry about it.

1 Like
3

Thank you so much Prof. @paul.buerkner , I highly appreciate your reply, this information is important.

I kindly ask about the first question, it is very critical for me, because I have a deep fear that excluding a model with unequal variances can negatively affect the validity and reliability of the model I choose.

1- Cumulative model with unequal variances works with a single outcome variable when I add the threshold = “equidistant”, but this did not work with multivariate model, is there any special syntax for the multivariate model.

4

Side question: what multivariate relationship does brm assume for the multiple outcome variables?

1 Like
5

None, unless you make your random effects correlated.

1 Like
6

I am afraid, I don’t fully understand. Where in your model is the “unequal variance” aspect? The way I understand modeling unequal variances in ordinal models is via modeling the disc discrimination parameter (see https://www.jstatsoft.org/article/view/v100i05 for details).

I don’t fully understand how that is related to your multivariate modeling aspect. The way you currently set up the multivariate model is equivalent to just 3 univariate models. I suggest you just start with the individual univariate models to make it easier to learn more about these models and also to spot problems in estimation.

1 Like
7

There is another paper (I’ve forgotten the reference) where the authors re-invent the constrained partial proportional odds model by modeling different variances for latent variables.

Related to the use of shared random effects across multivariate responses variables, have you seen a method that allows you to marginalize some way at the end? Example: a drug may reduce the chance of heart attack and the chance of a stroke. With a shared random effects the efficacy of the drug on one endpoint may be lessened because of the way it affects the other endpoint. The clinical trialists would want to know the effect on heart attack per se.

1 Like
8

I am not sure I would argue for the use of shared random effects (as in the same parameter appears in multiple responses) but for correlated random effects such that the information can be shared across responses. I don’t fully understand how this relates to the notion of marginalization here. Over what would we marginalize?

1 Like
9

Thanks for clarifying!

I am sorry for replying late, I just needed a time to read the valuable article you’ve shared, I also apologize for not being able to clarify what I am asking about.

After going through the article you mentioned, for ordinal models, I must add a discrimination parameter, for generalization purpose.

1- Is fitting the model this way correct considering adding the discrimination parameter?

library(brms)

fit_happiness5 ← brm(
bf(
mvbind(sign_pet, join_boycott, attend_demo, join_strike) ~ happiness * var + (1 |p| country),
disc ~ 1 + (1 |p| country)
),
data = scaled,
family = cumulative(“logit”),
chains = 4,
cores = 4,
iter = 2000,
control = list(adapt_delta = 0.99, max_treedepth = 12)
)

Note: I included correlations in this model

2- Am I correct about this? the generalization of multilevel ordinal models (particularly in case of my model) requires fitting:

  • Category-Specific Effects model “by using the function cs()” in my model I included the interaction between country level variable and individual independent variable, by writing this syntax “cs(financial_satis * var)” with acat family.
  • Adding discrimination parameters to the model this way “disc ~ 1 + (1 |p| country)”

I express my gratitude for the time you have given to answer my questions, this is an invaluable help

Respectfully,

10

Over the other endpoint, if possible. I’m thinking out loud but the goal is to estimate the effect on one endpoint ignoring the effect on the other endpoint. I’d like to have my cake and eat it too, i.e., be able to take between-endpoint correlations into account, but also have a way to estimate the effect on an endpoint as if the other endpoint was not in the model.

11

I am careful at calling any specific stats model “correct”.

  1. I think the model looks good.

  2. The acat family with cs() effects looks reasonable to me. The way you use the discrimination parameters looks reasonable too.

1 Like
12

I would have to think more about how to formalize this notion and how/if to express it in brms.

13

Thank you so much for such a helpful reply… brms package is a treasure as well

Respectfully,
MD