Post-hoc testing for cumulative link mixed-effects model with interactions in R

Question

I'm a resident physician working on my doctor's thesis and I'm trying to analyse data from a survey with R. I have basic mathematic and rookie statistics skills.

Participants of the survey have been shown four pictures of people (average man, average woman, attractive man, attractive woman) with four disfigurements (strabism, acne, piercing, tattoo), crossed in a latin square OR the control group (every face without disfigurement), randomized equally (20% each group). Participants were told they are gonna have a surgical treatment or a medical checkup (randomized) by the physician shown and they had to answer on a Likert-Scale from 0-10 (0 very improbable, 10 very probable) how much they would like to get the medical treatment from this physician. They also had to rate the physicians regarding to attractivity, competence, honesty, intelligence, kindness and reliability again on a Likert-Scale from 0-10 and also how important these attributes for the participants are.

Summary

participant_id      character       unique number
answer              ordered factor  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
attractivity        ordered factor  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
competence          ordered factor  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
honesty             ordered factor  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
intelligence        ordered factor  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
kindness            ordered factor  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
reliability         ordered factor  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
weight_attractivity ordered factor  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
weight_competence   ordered factor  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
weight_honesty      ordered factor  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
weight_intelligence ordered factor  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
weight_kindness     ordered factor  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
weight_reliability  ordered factor  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
situation           factor          Internist, Surgeon
face                factor          Man_Average, Woman_Average, Man_Attractive, Woman_Attractive
disfigurement       factor          None, Strabism, Acne, Piercing, Tattoo

I have 4009 unique participants with four answer blocks each, so 16036 rows in wide format in total.

I'm beginning with answer as dependent variable and situation, face and disfigurement as independent variables as well as participant_id as a random factor in a cumulative link mixed-effect model with clmm from the ordinal package for R. There seem to be significant interactions between situation, face and disfigurement, so I implemented the model as follows:

model.clmm <- clmm(answer ~ situation * disfigurement * face + (1 | participant_id), data = answers_full, Hess = TRUE, threshold = "symmetric")

I would like to do post-hoc tests with the emmeans package.

Question One: Am I (in general) on the right way with this strategy or totally wrong?

Question Two: Is it acceptable to do pairwise comparisons given the fact that there are interactions or what are the alternatives?

Question Three: To analyse all of the given answers, is it acceptable to do multiple univariate ordinal regressions or should I switch to multivariate ordinal regression (e.g. with the mvord package)?

Thank you very much in advance! Kindest regards, Pascal

Answer 1

Not sure I can adequately answer the stated questions, especially (3). But a few comments:

1. Perhaps you should have a group factor in your model to distinguish the control group (disfigurement of "none") from the other levels; because the design is different for the control group than for the others. emmeans() should detect a nesting structure of disfigurement nested in group. It may be OK to not have a group factor, but I think it'd help with interpretation.
2. You state that face and disfigurement are arranged in a Latin square, for the subjects in the non-control group. But in a Latin square, you give up the ability to estimate the interaction fully. Thus, I think you should exclude face : disfigurement and situation : face : disfigurement in the model.
3. Before plunging headlong into a sea of P values do some diagnostics to see how well the model fits the data. I do not have much experience with cumulative-link models, so I'm shaky on what those diagnostics might be; but you owe it to yourself to make sure the model you use actually explains the data and there are no glaring shortcomings; if the model is bad, then any post hoc analysis is garbage.
4. Again, before plunging into comparisons and $P$ values, look at the results graphically. For example, emmip(model.clmm, situation ~ disfigurement | face). This could help a lot in understanding which pairwise comparisons are actually going to be of interest. Some comparisons may be skipped simply because there are no practical differences in what you see.
5. If there are strong interactions, you should not do marginal comparisons of each factor. Those don't make sense because it doesn't make sense to average over factors that interact with the one of interest. On the other hand, if you do emm = emmeans(model.clmm, ~ situation * disfigurement * face); emm you will obtain estimates and CIs for all factor combinations; then if you subsequently do pairs(emm), you'll get pairwise comparisons of all of those estimates, which is an overwhelming number. You can reduce the number of comparisons considerably by doing separate comparisons of one factor, holding the others fixed: pairs(emm, simple = "each").
6. I'm no expert on ordinal models, but is the "symmetric" threshold justified? This is a pretty restrictive assumption, and it may not be what you think it is. It specifies that the cut points you estimate are arranged symmetrically; it does not specify that the underlying latent variable is symmetric; that requirement is already implicitly stated in your use of the default logistic link function.