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