Summary and Analysis of Extension Program Evaluation in R

Salvatore S. Mangiafico

One-way Repeated Ordinal Regression with CLMM

Appropriate data

•  Two-way data with one treatment variable and one blocking variable

•  Dependent variable is ordered factor

•  Treatment variable is a factor with at least two levels or groups

•  Blocking variable is a factor with at least two levels

•  Observations can be paired or repeated measures within blocks

Interpretation

A significant result can be interpreted as, “There was a significant difference among groups.”  Or, “There was a significant effect of Independent Variable.”

A significant post-hoc analysis indicates e.g. “There was a significant difference between Group A and Group B”, and so on.

Packages used in this chapter

The packages used in this chapter include:

•  psych

•  FSA

•  lattice

•  ordinal

•  car

•  RVAideMemoire

•  multcomp

•  emmeans

The following commands will install these packages if they are not already installed:

if(!require(psych)){install.packages("psych")}
if(!require(FSA)){install.packages("FSA")}

if(!require(lattice)){install.packages("lattice")}
if(!require(ordinal)){install.packages("ordinal")}
if(!require(car)){install.packages("car")}
if(!require(RVAideMemoire)){install.packages("RVAideMemoire")}
if(!require(multcomp)){install.packages("multcomp")}
if(!require(emmeans)){install.packages("emmeans")}

One-way repeated ordinal regression example

The following example revisits the Belcher family data from the Friedman Test chapter.

In this example, we want to determine if there is a difference in scores among five instructors from the Belcher family.  Because we know which rater gave each score, we will use the Rater as a blocking variable, with the suspicion that one rater might rate consistently higher than another one.  If we hadn’t recorded this information about raters, we could use one-way ordinal regression.

Data = read.table(header=TRUE, stringsAsFactors=TRUE, text="

 Instructor        Rater  Likert
 'Bob Belcher'        a      4
 'Bob Belcher'        b      5
 'Bob Belcher'        c      4
 'Bob Belcher'        d      6
 'Bob Belcher'        e      6
 'Bob Belcher'        f      6
 'Bob Belcher'        g     10
 'Bob Belcher'        h      6
 'Linda Belcher'      a      8
 'Linda Belcher'      b      6
 'Linda Belcher'      c      8
 'Linda Belcher'      d      8
 'Linda Belcher'      e      8
 'Linda Belcher'      f      7
 'Linda Belcher'      g     10
 'Linda Belcher'      h      9
 'Tina Belcher'       a      7
 'Tina Belcher'       b      5
 'Tina Belcher'       c      7
 'Tina Belcher'       d      8
 'Tina Belcher'       e      8
 'Tina Belcher'       f      9
 'Tina Belcher'       g     10
 'Tina Belcher'       h      9
 'Gene Belcher'       a      6
 'Gene Belcher'       b      4
 'Gene Belcher'       c      5
 'Gene Belcher'       d      5
 'Gene Belcher'       e      6
 'Gene Belcher'       f      6
 'Gene Belcher'       g      5
 'Gene Belcher'       h      5
 'Louise Belcher'     a      8
 'Louise Belcher'     b      7
 'Louise Belcher'     c      8
 'Louise Belcher'     d      8
 'Louise Belcher'     e      9
 'Louise Belcher'     f      9
 'Louise Belcher'     g      8
 'Louise Belcher'     h     10            
")


### Order levels of the factor; otherwise R will alphabetize them

Data$Instructor = factor(Data$Instructor,
                         levels=unique(Data$Instructor))


### Create a new variable which is the likert scores as an ordered factor

Data$Likert.f = factor(Data$Likert,
                       ordered=TRUE)


###  Check the data frame

library(psych)

headTail(Data)

str(Data)

summary(Data)

Summarize data treating Likert scores as factors

XT = xtabs( ~ Instructor + Likert.f,
           data = Data)

XT

                Likert.f
Instructor       4 5 6 7 8 9 10
  Bob Belcher    2 1 4 0 0 0  1
  Linda Belcher  0 0 1 1 4 1  1
  Tina Belcher   0 1 0 2 2 2  1
  Gene Belcher   1 4 3 0 0 0  0
  Louise Belcher 0 0 0 1 4 2  1

prop.table(XT,
           margin = 1)

                Likert.f
Instructor           4     5     6     7     8     9    10
  Bob Belcher    0.250 0.125 0.500 0.000 0.000 0.000 0.125
  Linda Belcher  0.000 0.000 0.125 0.125 0.500 0.125 0.125
  Tina Belcher   0.000 0.125 0.000 0.250 0.250 0.250 0.125
  Gene Belcher   0.125 0.500 0.375 0.000 0.000 0.000 0.000
  Louise Belcher 0.000 0.000 0.000 0.125 0.500 0.250 0.125

### Proportion with each row

Bar plots by group

Note that the bar plots don’t show the effect of the blocking variable.

library(lattice)

histogram(~ Likert.f | Instructor,
          data=Data,
          layout=c(1,5)      #  columns and rows of individual plots
          )

Summarize data treating Likert scores as numeric

library(FSA)

Summarize(Likert ~ Instructor,
          data=Data,
          digits=3)

      Instructor n  mean    sd min   Q1 median   Q3 max percZero
1    Bob Belcher 8 5.875 1.885   4 4.75      6 6.00  10        0
2  Linda Belcher 8 8.000 1.195   6 7.75      8 8.25  10        0
3   Tina Belcher 8 7.875 1.553   5 7.00      8 9.00  10        0
4   Gene Belcher 8 5.250 0.707   4 5.00      5 6.00   6        0
5 Louise Belcher 8 8.375 0.916   7 8.00      8 9.00  10        0

One-way repeated ordinal regression

The model is specified using formula notation with the clmm function.  Likert.f is the dependent variable, and an ordered factor in R, and Instructor is the independent variable.  Rater is the blocking variable, identified as the random variable in the mixed model.  The data= option indicates the data frame that contains the variables.  For the meaning of other options, see ?clmm.

is determined using the Anova function in the packages RVAideMemoire and car.

Define model

library(ordinal)

model = clmm(Likert.f ~ Instructor + (1|Rater),
             data = Data)

Analysis of deviance

anova(model, type="II")

Error in eval(predvars, data, env) : object 'Likert.f' not found

   ### I don’t know why this error occurs

Comparison of models analysis

model.null = clmm(Likert.f ~ 1 + (1|Rater),
                 data=Data)

anova(model, model.null)

Likelihood ratio tests of cumulative link models:

           no.par    AIC logLik LR.stat df Pr(>Chisq)   
model.null      7 163.00 -74.50                         
model          11 134.94 -56.47  36.059  4  2.814e-07 ***

### The p-value is for the effect of Instructor

Alternate analysis of deviance analysis

library(car)

library(RVAideMemoire)

Anova.clmm(model, type = "II")

Analysis of Deviance Table (Type II tests)

           LR Chisq Df Pr(>Chisq)   
Instructor   36.059  4  2.814e-07 ***

Determine the effect of the random variable

model.fixed = clm(Likert.f ~ Instructor,
                  data = Data)

anova(model,
       null = model.fixed)

Likelihood ratio tests of cumulative link models:

            no.par    AIC  logLik LR.stat df Pr(>Chisq)  
model.fixed     10 140.74 -60.368                        
model           11 134.94 -56.470   7.795  1   0.005239 **

   ### p-value is effect of Rater

Post-hoc test with emmeans

Because the main effect of Instructor in the model was significant, we want to know which instructor has statistically different scores than which other instructors.  Comparing EM means for each instructor will allow us to do this.

For clmm model objects, the emmean, SE, LCL, and UCL values should be ignored, unless specific options in emmeans are selected.

library(emmeans)

marginal = emmeans(model,
                   ~ Instructor)

pairs(marginal,
      adjust="tukey")

contrast                       estimate    SE  df z.ratio p.value
 Bob Belcher - Linda Belcher      -4.641 1.326 Inf  -3.500  0.0042
 Bob Belcher - Tina Belcher       -4.381 1.284 Inf  -3.413  0.0058
 Bob Belcher - Gene Belcher        1.036 1.026 Inf   1.010  0.8511
 Bob Belcher - Louise Belcher     -5.376 1.398 Inf  -3.845  0.0011
 Linda Belcher - Tina Belcher      0.260 0.914 Inf   0.284  0.9986
 Linda Belcher - Gene Belcher      5.677 1.434 Inf   3.958  0.0007
 Linda Belcher - Louise Belcher   -0.735 0.908 Inf  -0.809  0.9279
 Tina Belcher - Gene Belcher       5.417 1.385 Inf   3.911  0.0009
 Tina Belcher - Louise Belcher    -0.994 0.919 Inf  -1.082  0.8160
 Gene Belcher - Louise Belcher    -6.412 1.514 Inf  -4.235  0.0002

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

library(multcomp)

cld(marginal, Letters=letters)

Instructor     emmean    SE  df asymp.LCL asymp.UCL .group
 Gene Belcher    -3.59 1.104 Inf   -5.7539    -1.427  a   
 Bob Belcher     -2.55 1.038 Inf   -4.5878    -0.521  a   
 Tina Belcher     1.83 0.909 Inf    0.0456     3.608   b  
 Linda Belcher    2.09 0.937 Inf    0.2490     3.924   b  
 Louise Belcher   2.82 0.980 Inf    0.9011     4.742   b  

Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05

NOTE: Compact letter displays can be misleading
      because they show NON-findings rather than findings.
      Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.

   ### Groups sharing a letter in .group are not significantly different

   ### Remember to ignore “emmean”, “SE”, “LCL”, and “UCL” for CLM

Check model assumptions

At the time of writing, the nominal_test and scale_test functions don’t work with clmm model objects, so we will define a similar model with clm and no random effects, and test the proportional odds assumption on this model.  In this example, neither the nominal_test nor scale_test functions produced a p-value, so the results of these tests are not conclusive.

model.clm = clm(Likert.f ~ Instructor + Rater,
            data = Data)

nominal_test(model.clm)

Tests of nominal effects

           Df  logLik    AIC LRT Pr(>Chi)
<none>        -45.907 125.81             
Instructor                              
Rater

   ### No p-value produced for this example.

scale_test(model.clm)

Tests of scale effects

           Df  logLik    AIC LRT Pr(>Chi)
<none>        -45.907 125.81            
Instructor                               
Rater

Warning messages:
1: (-2) Model failed to converge: degenerate Hessian with 1 negative eigenvalues
In addition: maximum number of consecutive Newton modifications reached
2: (-1) Model failed to converge with max|grad| = 6.85558e-06 (tol = 1e-06)
In addition: maximum number of consecutive Newton modifications reached

   ### No p-value produced for this example.
   ### Errors produced for this example.

©2016 by Salvatore S. Mangiafico.
Rutgers Cooperative Extension, New Brunswick, NJ.
Updated: 2025.

Non-commercial reproduction of this content, with attribution, is permitted.
For-profit reproduction without permission is prohibited.

If you use the code or information in this site in a published work, please cite it as a source.  Also, if you are an instructor and use this book in your course, please let me know.   My contact information is on the About the Author of this Book page.

Citation

Mangiafico, S.S. 2016. Summary and Analysis of Extension Program Evaluation in R, version 1.23.0, revised 2025. rcompanion.org/handbook/. (Pdf version: rcompanion.org/documents/RHandbookProgramEvaluation.pdf.)