[banner]

Summary and Analysis of Extension Program Evaluation in R

Salvatore S. Mangiafico

One-way Ordinal Regression with CLM

Appropriate data

•  One-way data with two or more groups

•  Dependent variable is ordered factor

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

•  Observations between groups are not paired or repeated measures

 

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, “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

•  multcompView

•  lsmeans

•  rcompanion

 

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(multcompView)){install.packages("multcompView")}
if(!require(lsmeans)){install.packages("lsmeans")}
if(!require(rcompanion)){install.packages("rcompanion")}


One-way ordinal regression example


Input =("
 Speaker  Likert
 Pooh      3
 Pooh      5
 Pooh      4
 Pooh      4
 Pooh      4
 Pooh      4
 Pooh      4
 Pooh      4
 Pooh      5
 Pooh      5
 Piglet    2
 Piglet    4
 Piglet    2
 Piglet    2
 Piglet    1
 Piglet    2
 Piglet    3
 Piglet    2
 Piglet    2
 Piglet    3
 Tigger    4
 Tigger    4
 Tigger    4
 Tigger    4
 Tigger    5
 Tigger    3
 Tigger    5
 Tigger    4
 Tigger    4
 Tigger    3
")

Data = read.table(textConnection(Input),header=TRUE)


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


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


### 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)


### Remove unnecessary objects

rm(Input)


Summarize data treating Likert scores as factors


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


        Likert.f
Speaker  1 2 3 4 5
  Pooh   0 0 1 6 3
  Piglet 1 6 2 1 0
  Tigger 0 0 2 6 2


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

prop.table(XT,
           margin = 1)


        Likert.f
Speaker    1   2   3   4   5
  Pooh   0.0 0.0 0.1 0.6 0.3
  Piglet 0.1 0.6 0.2 0.1 0.0
  Tigger 0.0 0.0 0.2 0.6 0.2


Bar plots by group


library(lattice)

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

image


Summarize data treating Likert scores as numeric

 

library(FSA)

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


  Speaker  n mean    sd min Q1 median   Q3 max percZero
1    Pooh 10  4.2 0.632   3  4      4 4.75   5        0
2  Piglet 10  2.3 0.823   1  2      2 2.75   4        0
3  Tigger 10  4.0 0.667   3  4      4 4.00   5        0


One-way ordinal regression

The model is specified using formula notation.  Here, Likert.f is the dependent variable and Speaker is the independent variable.  The data= option indicates the data frame that contains the variables.  For the meaning of other options, see ?clm.

 

The p-value for the effect of Speaker is determined using the Anova function in the packages RVAideMemoire and car.

 

 

Define model

 

library(ordinal)

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


Analysis of deviance


library(car)

library(RVAideMemoire)

Anova(model,
      type = "II")


Analysis of Deviance Table (Type II tests)

        LR Chisq Df Pr(>Chisq)   
Speaker   23.395  2  8.315e-06 ***


Post-hoc test with lsmeans

In the lsmeans function, model specifies the model object that was previously fitted.  Note the specialized formula where pairwise indicates that all pairwise comparisons should be conducted, and Instructor indicates the variable whose levels will be compared.

Here we’ll create an object of the lsmeans output called leastsquare.  Then we’ll pass this object to the cld function to create a compact letter display.

 

Be sure to read the Least Square Means for Multiple Comparisons chapter for correct interpretation of least square means.  For clm model objects, the lsmean, SE, LCL, and UCL values should be ignored, unless specific options in lsmeans are selected.


library(multcompView)

library(lsmeans)

leastsquare = lsmeans(model,
                      pairwise ~ Speaker,
                      adjust="tukey")      ###  Tukey-adjusted comparisons

leastsquare


cld(leastsquare,
    alpha=0.05,
    Letters=letters,      ### Use lower-case letters for .group
    adjust="tukey")       ### Tukey-adjusted comparisons


Speaker    lsmean        SE df  asymp.LCL  asymp.UCL .group
 Piglet  -1.783599 0.7755958 NA -3.6355182 0.06832068  a   
 Tigger   2.526492 0.8413200 NA  0.5176407 4.53534421   b  
 Pooh     3.160223 0.8912968 NA  1.0320404 5.28840661   b  

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

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

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


Check model assumptions

 

nominal_test(model)


Tests of nominal effects

formula: Likert.f ~ Speaker
        Df  logLik    AIC LRT Pr(>Chi)
<none>     -30.149 72.298            
Speaker   

   ###  No p-value


scale_test(model)

Tests of scale effects

formula: Likert.f ~ Speaker
        Df  logLik    AIC     LRT Pr(>Chi)
<none>     -30.149 72.298                
Speaker  2 -29.839 75.677 0.62096   0.7331

   ###  No violation in assumptions


Optional analyses

 

Post-hoc test: pairwise ordinal tests for multiple comparisons of groups

 

Table format

Post-hoc testing can be conducted with the pairwiseOrdinalTest function, which produces a table of p-values comparing each pair of groups, or with the pairwiseOrdinalMatrix function in the rcompanion package, which produces a matrix of p-values.  This can then be converted into a compact letter display.

 

To prevent the inflation of type I error rates, adjustments to the p-values can be made using the p.adjust.method option. See ?p.adjust for details on available p-value adjustment methods.


### Order groups by median

Data$Speaker = factor(Data$Speaker,
                      levels=c("Pooh", "Tigger", "Piglet"))


### Pairwise ordinal regression tests

library(rcompanion)

PT = pairwiseOrdinalTest(Likert.f ~ Speaker,
                         data   = Data,
                         method = "fdr")
                              # Adjusts p-values for multiple comparisons;
                              # See ?p.adjust for options

PT


           Comparison   p.value  p.adjust
1   Pooh - Tigger = 0    0.4763 4.763e-01
2   Pooh - Piglet = 0 2.866e-05 8.598e-05
3 Tigger - Piglet = 0 9.477e-05 1.422e-04


#### Compact letter display

library(rcompanion)

cldList(p.adjust  ~ Comparison,
        data       = PT,
        threshold  = 0.05)


   Group Letter MonoLetter
1   Pooh      a         a
2 Tigger      a         a
3 Piglet      b          b

Groups sharing a letter are not significantly different (alpha = 0.05).


Matrix format and compact letter display

Usually you would want to first order your treatment variable (Speaker) by median or other location statistic for the letters in the compact letter display to be in their proper order.


### Order groups by median

Data$Speaker = factor(Data$Speaker,
                      levels=c("Pooh", "Tigger", "Piglet"))


### Pairwise ordinal regression tests

library(rcompanion)

PM = pairwiseOrdinalMatrix(Likert.f ~ Speaker,
                           data   = Data,
                           method = "fdr")
                               # Adjusts p-values for multiple comparisons;
                               # See ?p.adjust for options

PM


$Adjusted

            Pooh    Tigger    Piglet
Pooh   1.000e+00 0.4763000 8.598e-05
Tigger 4.763e-01 1.0000000 1.422e-04
Piglet 8.598e-05 0.0001422 1.000e+00


library(multcompView)

multcompLetters(PM$Adjusted,
                compare="<",
                threshold=0.05,  # p-value to use as significance threshold
                Letters=letters,
                reversed = FALSE)


  Pooh Tigger Piglet
   "a"    "a"    "b"

   ### Values sharing a letter are not significantly different