I ran an ordinal logistic regression (using the function clmm from the R package ordinal) with a two-factor interaction and a random effect. The response is a factor w/ 5 levels (Liker scale: 1 2 3 4 5), the independent variables are a factor w/ 2 levels (time) and a factor w/ 3 levels (group)

The code looks like this:

library(ordinal)

# dataset 

ID         time   group  response
person1      1      a       3
person2      1      a       5
person3      1      c       5
person4      1      b       2
person5      1      c       2
person6      1      a       4
person1      2      a       2
person2      2      a       2
person3      2      c       1
person4      2      b       4
person5      2      c       3
person6      2      a       4
...         ...    ...     ...

# model  

model <- clmm(response ~ time*group + (1|ID))

# model results

formula: response ~ time * group + (1 | ID)
data:    dataset

 link  threshold nobs logLik  AIC    niter     max.grad cond.H 
 logit flexible  168  -226.76 473.52 508(4150) 9.42e-05 1.8e+02

Random effects:
 Groups Name        Variance Std.Dev.
 ID     (Intercept) 5.18     2.276   
Number of groups:  ID 84 

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
time2                0.2837     0.5289   0.536  0.59170    
group_b              1.8746     0.6946   2.699  0.00695 ** 
group_c              4.0023     0.9383   4.265    2e-05 ***
time2:group_b       -0.5100     0.7294  -0.699  0.48447    
time2:group_c       -0.8830     0.9749  -0.906  0.36508    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Threshold coefficients:
    Estimate Std. Error z value
1|2  -2.6223     0.6440  -4.072
2|3   0.2474     0.5427   0.456
3|4   2.5384     0.5824   4.359
4|5   4.6786     0.7143   6.550

As you can see, the model results show only whether there are differences compared to the intercept (i.e. time1:group_a). However, what I am also interested in is to check if the difference between time1:group_b and time2:group_b is statistically significant, same for group_c.

Since I have to account for the random effect, I cannot use a simple chi-square test to check for statistically significant differences between groups. I therefore tried to run the function contrast from the R package emmeans, which uses the output of the function emmeans, see the code below:

library(emmeans)

em <- emmeans(model, ~ time | group) #calculates the estimated marginal means

contrast(em, "consec", simple = "each")

# contrast results

$`simple contrasts for time`
group = a:
 contrast estimate    SE  df z.ratio p.value
 2 - 1       0.284 0.529 Inf  0.536  0.5917 

group = b:
 contrast estimate    SE  df z.ratio p.value
 2 - 1      -0.226 0.482 Inf -0.470  0.6386 

group = c:
 contrast estimate    SE  df z.ratio p.value
 2 - 1      -0.599 0.816 Inf -0.734  0.4629 

Note: contrasts are still on the as.factor scale 

$`simple contrasts for group`
time = 1:
 contrast estimate    SE  df z.ratio p.value
 b - a        1.87 0.695 Inf 2.699   0.0137 
 c - b        2.13 0.871 Inf 2.443   0.0284 

time = 2:
 contrast estimate    SE  df z.ratio p.value
 b - a        1.36 0.687 Inf 1.986   0.0897 
 c - b        1.75 0.838 Inf 2.095   0.0695 

Note: contrasts are still on the as.factor scale 
P value adjustment: mvt method for 2 tests 

My questions are: a) Is this a correct and valid method to check whether the differences are significant? b) If not, what is the correct way to do this?

Of course any other suggestion is extremely welcome! Thanks a lot.

0

There are 0 best solutions below