Interpreting the results of a clmm model with one independent variable - Some questions

Question

I'm just beginning to use R and learning about statistics, so please bear with me. Questions are at the end, but I'll give my interpretation whilst I show the code.

I ran an perception experiment in which 21 subjects had to listen to some sentences and had to judge from 0 to 3 the degree of proeminence on each syllable of said sentence. Each of these syllables are associated with a type of accent (NONE, AI, AF, AFD).

I am interested in knowing if the response (ordinal variable with 4 levels: 0, 1, 2, 3, which could also be read as "no accent", "weak accent", "medium accent", "strong accent") is influenced by the type of accent (categorical variable with 4 levels: NONE, AI, AF, AFD).

I decided to run a cummulative mixed effects model where my variables are:

Dependent variable: score (this is the response)
Independent variable: accent
Random factors: auditeur (listener), item, locuteur (speaker)

This is my code:

my_data$score <- as.ordered(my_data$score)
my_data$accent_position <- as.factor(my_data$accent)
my_data$auditeur <- as.factor(my_data$auditeur)
my_data$item <- as.factor(my_data$item)
my_data$locuteur <- as.factor(my_data$locuteur)

library(ordinal)

m1 <- clmm(score ~ accent + (1|auditeur) + (1|item) + (1|locuteur), data = my_data)

This is the output:

> summary(m1)
Cumulative Link Mixed Model fitted with the Laplace approximation

formula: score ~ accent + (1 | auditeur) + (1 | item) + (1 | locuteur)
data:    my_data

 link  threshold nobs logLik   AIC      niter     max.grad cond.H 
 logit flexible  6468 -6541.36 13100.72 798(3196) 1.99e-03 6.2e+02

Random effects:
 Groups   Name        Variance Std.Dev.
 item     (Intercept) 0.06828  0.2613  
 auditeur (Intercept) 1.57589  1.2553  
 locuteur (Intercept) 0.32135  0.5669  
Number of groups:  item 36,  auditeur 21,  locuteur 6 

Coefficients:
           Estimate Std. Error z value Pr(>|z|)    
accentAFD  -0.72078    0.09166  -7.864 3.73e-15 ***
accentAI    0.26449    0.06770   3.907 9.36e-05 ***
accentNONE -2.95239    0.07725 -38.219  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Threshold coefficients:
    Estimate Std. Error z value
0|1  -1.8949     0.3669  -5.165
1|2   0.1063     0.3659   0.291
2|3   2.0375     0.3670   5.551

From that summary, I interpret that the three types of accents (AFD, AI, NONE ) are significantly different from the intercept which is AF. In order to observe the interactions, I run a post-hoc test with emmeans:

library(emmeans)

pairwise <- emmeans(m1, pairwise~accent)

> pairwise
$emmeans
 accent emmean    SE  df asymp.LCL asymp.UCL
 AF     -0.083 0.365 Inf    -0.799    0.6333
 AFD    -0.804 0.369 Inf    -1.527   -0.0805
 AI      0.182 0.364 Inf    -0.533    0.8956
 NONE   -3.035 0.366 Inf    -3.754   -2.3171

Confidence level used: 0.95 

$contrasts
 contrast   estimate     SE  df z.ratio p.value
 AF - AFD      0.721 0.0917 Inf   7.864 <.0001 
 AF - AI      -0.264 0.0677 Inf  -3.907 0.0005 
 AF - NONE     2.952 0.0773 Inf  38.219 <.0001 
 AFD - AI     -0.985 0.0850 Inf -11.590 <.0001 
 AFD - NONE    2.232 0.0904 Inf  24.677 <.0001 
 AI - NONE     3.217 0.0742 Inf  43.338 <.0001 

P value adjustment: tukey method for comparing a family of 4 estimates 

The means seem to indicate that subjects perceived (from 0 to 3) the accents as follow: AI > AF > AFD > NONE.

Then, the contrasts show me that these differences below are statistically significant:

AF > AFD
AI > AF
AI > AFD
AF > NONE
AFD > NONE
AI > NONE

Now my questions :

1. Could someone please tell me if this interpretation is correct and if the model is well made for what I am trying to show? I'm having trouble understanding for instance why the means are negative and positive. Am I missing something?

2. I read a lot about using mean.class and polr but I don't know if that's necessary for my model and for what I am trying to show. I feel a little bit overwhelmed by the amount of different models you can apply when having an ordinal dependent variable and I am not sure about having understood everything.

Thank you so much in advance.

Answer 1

The ordinal model you fitted is based on the idea that there is a variable $y$ having a logistic distribution (slightly heavier-tailed than normal), whose mean depends on the predictors. But instead of being able to observe $y$ directly, we only observe which interval $(-\infty, c_1], (c_1, c_2], (c_2, c_3], (c_3, \infty)$ it falls in. That makes $y$ a latent variable (exists but not observed). The cut points $c_j$ are estimated by the model-fitting software, in this case they are about $-1.89, 0.11, 2.04$ respectively as shown in the "thresholds" part of the output.

Notice that for any constant $a$, if we were to replace $y$ by $y+a$, and $c_j$ by $c_j+a$, then nothing would change in terms of the probabilities of falling in the respective intervals. So what is done (I think) for sake of identifiability of the parameters is to assume that the overall average of the $y$ values is zero. And that makes about half of the predicted $y$ values negative and the other half positive.

By default, emmeans() estimates the marginal means of this latent $y$ variable; so, as just explained, some of those estimates are negative and some are positive.

Another option is to use emmeans(..., mode = "mean.class"). In that case, we obtain the estimated distributions of the four ordinal levels for each accent and obtain the mean of that distribution. So this is something like what you would obtain by fitting a regression model to the ordinal response as a numerical variable. The catch is that this presumes the levels are numbered 1, 2, 3, and 4; and so those means will be one higher than you might expect, given that you numbered them 0, 1, 2, and 3. I'd suggest you stick with the latent-variable parameterization anyway.

The mixed model you have fitted also assumes there are random effects for listener, item, and speaker. Those random effects are on the latent scale; that is, the latent variable $y$ is subject to additional random variations.