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

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

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.