![[banner]](images/banner.jpg "Cohansey River, Fairfield, New Jersey")
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
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)")
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.