[banner]

Summary and Analysis of Extension Program Evaluation in R

Salvatore S. Mangiafico

Two-sample Paired Ordinal Test with CLMM

 

For cumulative link models with random effects, the clmm function is used instead of the clm function.  The clmm function specifies a mixed effects model.  This type of model is appropriate for paired and repeated measures analyses.

 

With these models, the order of the data doesn’t matter as it did in the paired signed-rank test example, because here the blocking variable, Student, is entered explicitly in the model.

 

Appropriate data

•  Two-sample data.  That is, one-way data with two groups only

•  Dependent variable is ordered factor

•  Independent variable is a factor with two levels.  That is, two groups

•  Observations between groups are paired

 

Interpretation

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

 

Packages used in this chapter

 

The packages used in this chapter include:

•  psych

•  ordinal

•  car

•  RVAideMemoire

 

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

 

if(!require(psych)){install.packages("psych")}
if(!require(ordinal)){install.packages("ordinal")}
if(!require(car)){install.packages("car")}
if(!require(RVAideMemoire)){install.packages("RVAideMemoire")}

 

Two-sample paired ordinal model example

 

The following example revisits the Pooh data from the Two-sample Paired Signed-rank Test chapter.

 

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

 Speaker  Time  Student  Likert
 Pooh      1     a        1
 Pooh      1     b        4
 Pooh      1     c        3
 Pooh      1     d        3
 Pooh      1     e        3
 Pooh      1     f        3
 Pooh      1     g        4
 Pooh      1     h        3
 Pooh      1     i        3
 Pooh      1     j        3
 Pooh      2     a        4
 Pooh      2     b        5
 Pooh      2     c        4
 Pooh      2     d        5
 Pooh      2     e        4
 Pooh      2     f        5
 Pooh      2     g        3
 Pooh      2     h        4
 Pooh      2     i        3
 Pooh      2     j        4
")


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


Plot the paired data

 

Scatter plot with one-to-one line

Paired data can be visualized with a scatter plot of the paired cases.  In the plot below, points that fall above and to the left of the blue line indicate cases for which the value for Time 2 was greater than for Time 1.

 

Note that the points in the plot are jittered slightly so that points that would fall directly on top of one another can be seen.

 

First, two new variables, Time.1 and Time.2, are created by extracting the values of Likert for observations with the Time variable equal to 1 or 2, respectively.

 

Note that for this plot, the data must be ordered so that the first observation where Time = 1 is paired to the first observation where Time = 2, and so on.

 

Time.1 = Data$Likert[Data$Time==1]
Time.2 = Data$Likert[Data$Time==2]
                      
plot(Time.1, jitter(Time.2),   # jitter offsets points so you can see them all
     pch = 16,                 # shape of points
     cex = 1.0,                # size of points
     xlim=c(1, 5.5),           # limits of x axis
     ylim=c(1, 5.5),           # limits of y axis
     xlab="Time 1",
     ylab="Time 2"
     )
abline(0,1, col="blue", lwd=2) # line with intercept of 0 and slope of 1

 

image

 

 

Bar plot of differences

Note that for this plot, the data must be ordered so that the first observation where Time = 1 is paired to the first observation where Time = 2, and so on.

 

Time.1 = Data$Likert[Data$Time==1]
Time.2 = Data$Likert[Data$Time==2]

Difference = Time.2 - Time.1

barplot(Difference,  
        col="dark gray",
        xlab="Observation",
        ylab="Difference (Time 2 – Time 1)")

 

image


Two-sample paired ordinal model

The model is specified using formula notation with the clmm function (not clm).  Here, Likert.f is the dependent variable and Time is the independent variable.  Student is specified as the blocking variable.  It is identified as a random variable in the mixed effects model.  The data= option indicates the data frame that contains the variables.  For the meaning of other options, see ?clmm.

 

Define model


library(ordinal)

model = clmm(Likert.f ~ Time + (1|Student),
             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|Student),
                 data=Data)

anova(model, model.null)


Likelihood ratio tests of cumulative link models:

           no.par    AIC  logLik LR.stat df Pr(>Chisq)  
model.null      4 54.445 -23.222                        
model           5 47.242 -18.621  9.2033  1   0.002416 **

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


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)  
Time   9.2033  1   0.002416 **


Determine the effect of the random variable

The effect of the random variable (Student, in this case) can be determined by comparing the full model to a model with only the fixed effects (Time, in this case). 


First, model.fixed is defined, and then the two models are passed to the anova function.


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

anova(model,
       null = model.fixed)


           no.par    AIC  logLik LR.stat df Pr(>Chisq)
model.null      4 45.262 -18.631                     
model           5 47.242 -18.621    0.02  1     0.8875

   ### p-value is effect of Student


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 ~ Time + Student,
            data = Data)

nominal_test(model.clm)


Tests of nominal effects

formula: Likert.f ~ Time + Student
        Df  logLik    AIC LRT Pr(>Chi)
<none>     -11.821 49.642            
Time                                  
Student

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



scale_test(model.clm)


Tests of scale effects

formula: Likert.f ~ Time + Student

        Df  logLik    AIC LRT Pr(>Chi)
<none>     -11.821 49.642            
Time                                 
Student                              

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.


Effect size statistics

 

Appropriate effect size statistics may be those used for the paired sign test or perhaps those used for the paired signed-rank test.