Ordinal Data
Dr Lauren M Ackerman
02 April 2019
What is “ordinal” data?
Data like reaction times is continuous because there are no discrete boundaries. RTs can be seconds, fractions of sections, fractions of milliseconds, or even finer grained than that.
Data like counts of words in corpora are discrete, but hypothetically infinite. You could have 0 attestations, 1, or one million… but you could not have 13.5 attestations of a word.
Binary data is also discrete, and can come from looking at hits vs misses (e.g., correct or incorrect responses) among other things. This data type is binary because there are only two values you can observe: 1 (hit) or 0 (miss).
But there is another type of data we use often, too. Ordinal data refers to data that is discrete but also ranked like in a Likert scale task. What makes this type of data unique is that it has a maximum and minimum endpoint, but also a limited number of (ordered) values in between. In this respect, binary data and ordinal data are similar, but when you have more than two possible responses, analysis can get quite complicated.
In particular, what makes analysis of ordinal data different from continuous, count, and binary data is that the perceptual size of each level may vary. That is, on a 7-point scale, the difference between a rating of 1 and a rating of 2 may not be the same distance as the difference between a rating of 2 and a rating of 3, or between a rating of 6 and a rating of 7. The distance between each level of the scale can vary across participants and items, so when we analyse this data, we really should take into account this dimension of variability.
That’s where the ordinal package comes in.
library(ordinal)
dataOrd <- read.csv("data/ordinal_data.csv",header=TRUE)
head(dataOrd)## subject item factor1 factor2 rating
## 1 1 1 Baseline Absent 5
## 2 1 2 Baseline Absent 6
## 3 1 3 Baseline Absent 7
## 4 1 4 Baseline Absent 7
## 5 1 5 Baseline Absent 7
## 6 1 6 Baseline Absent 7When we graph our rating data in a box plot, we see there is likely to be a main effect of Factor 1 and possibly an interaction between Factor 1 and Factor 2. It’s hard to tell more at this point.
ggplot(dataOrd, aes(x=factor1,y=rating)) + geom_boxplot(aes(fill=factor2)) + scale_y_continuous(breaks=c(1:7)) + ggtitle("Latin square rating data")
First, we can look at only the main effect of Factor 1. Visually speaking, it looks likely to be significant because the interquartile ranges do not overlap (although this is just a rule of thumb).
ggplot(dataOrd, aes(x=factor1,y=rating)) + geom_boxplot(aes(fill=factor1)) + scale_y_continuous(breaks=c(1:7)) + ggtitle("Main effect of Factor 1")
Turning to Factor 2, it’s not clear whether these two levels differ or not. The interquartile ranges largely overlap and the median bars appear to overlap the other group’s interquartile range. Simply put, this suggests the distributions don’t differ very much at this level of granularity. (Still, this is a rule of thumb.)
ggplot(dataOrd, aes(x=factor2,y=rating)) + geom_boxplot(aes(fill=factor2)) + scale_y_continuous(breaks=c(1:7)) + ggtitle("Main effect of Factor 2")
However, we have seen in the interaction plot that the influence of Factor 2 appears to be larger within the Treatment condition of Factor 1, so we might anticipate there being an interaction between the two main effects. We should conduct ordinal regression to see if this is the case.
Ordinal regressions
When we summarise a linear mixed effect model of this rating data, this is what we get.
lmer.rate.max <- lmer(rating ~ factor1*factor2 + (1|subject)+(1|item),data=dataOrd, REML=FALSE)
summary(lmer.rate.max)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: rating ~ factor1 * factor2 + (1 | subject) + (1 | item)
## Data: dataOrd
##
## AIC BIC logLik deviance df.resid
## 2892.6 2926.7 -1439.3 2878.6 953
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4690 -0.8577 -0.0005 0.8743 3.6434
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.071283 0.26699
## item (Intercept) 0.001825 0.04272
## Residual 1.128253 1.06219
## Number of obs: 960, groups: subject, 40; item, 24
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 6.02083 0.08099 74.342
## factor1Treatment -1.65833 0.09696 -17.102
## factor2Present -0.07083 0.09696 -0.731
## factor1Treatment:factor2Present -1.65417 0.13713 -12.063
##
## Correlation of Fixed Effects:
## (Intr) fctr1T fctr2P
## fctr1Trtmnt -0.599
## factr2Prsnt -0.599 0.500
## fctr1Trt:2P 0.423 -0.707 -0.707As in typical LMERs that one might produce using lme4, we can see the types of output, but it’s not looking good for Factor 2. We’re presented with various residuals and a summary of the random effects (which are important but I will skip for now). If we conducted a proper comparison of nested models with this dataset, we would get singular fits for each of the depleted models, but we would see that Factor 2 does appear to have a significant contribution. But, in any case, let’s see what happens when we analyze it with a CLMM instead of LMER.
The following model is a cumulative link mixed model (CLMM), which operates in a very similar way to a linear mixed effect model on the surface, but takes into account the possibility that each level of the rating scale could be a different size or distance from its neighbors.
clmm.max <- clmm(as.factor(rating) ~ factor1*factor2 + (1|subject)+(1|item),data=dataOrd)## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.summary(clmm.max)## Cumulative Link Mixed Model fitted with the Laplace approximation
##
## formula: as.factor(rating) ~ factor1 * factor2 + (1 | subject) + (1 |
## item)
## data: dataOrd
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 960 -1377.05 2776.10 1269(3796) 9.11e-04 6.4e+01
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.17920 0.4233
## item (Intercept) 0.01514 0.1230
## Number of groups: subject 40, item 24
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## factor1Treatment -2.7711 0.1974 -14.040 <2e-16 ***
## factor2Present -0.1735 0.1672 -1.038 0.299
## factor1Treatment:factor2Present -2.3701 0.2564 -9.243 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2 -6.7035 0.2714 -24.703
## 2|3 -5.4889 0.2400 -22.867
## 3|4 -4.1191 0.2131 -19.332
## 4|5 -2.6917 0.1826 -14.745
## 5|6 -0.9258 0.1480 -6.254
## 6|7 0.5231 0.1458 3.588You can see this in the section labeled Threshold coefficients, which shows the boundary estimates for each of the ordinal levels. (Note: you could probably plot these, too, just as I plot the estimates later on, but I won’t do that here.)
Just like you might with a LMER, you can do a nested model comparison to determine the contribution of each main effect and their interaction.
clmm.int <- clmm(as.factor(rating) ~ factor1+factor2 + (1|subject)+(1|item),data=dataOrd)## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.clmm.one <- clmm(as.factor(rating) ~ factor2 + (1|subject)+(1|item),data=dataOrd)
clmm.two <- clmm(as.factor(rating) ~ factor1 + (1|subject)+(1|item),data=dataOrd)First, by removing the interaction effect and comparing the depleted model to the maximal model, we find that the interaction contributes significantly to the overall model fit. This is expected, since we could visually detect something that looks like an interaction.
anova(clmm.max,clmm.int)## Likelihood ratio tests of cumulative link models:
##
## formula:
## clmm.int as.factor(rating) ~ factor1 + factor2 + (1 | subject) + (1 | item)
## clmm.max as.factor(rating) ~ factor1 * factor2 + (1 | subject) + (1 | item)
## link: threshold:
## clmm.int logit flexible
## clmm.max logit flexible
##
## no.par AIC logLik LR.stat df Pr(>Chisq)
## clmm.int 10 2863.5 -1421.8
## clmm.max 11 2776.1 -1377.0 89.439 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Next, we can look for the main effect of Factor 1. Note, however, that we compared this depleted model not to the maximal model but to the model with the interaction removed. This is because we should not check for a main effect by removing the main effect of one factor while also comparing it to the model with both the main effect and interaction. That would be counting the main effect twice, effectively. Technically, it is possible to do so, but the degrees of freedom with be 2. This is less informative because it includes two types of effects in the one model comparison.
anova(clmm.int,clmm.one)## Likelihood ratio tests of cumulative link models:
##
## formula:
## clmm.one as.factor(rating) ~ factor2 + (1 | subject) + (1 | item)
## clmm.int as.factor(rating) ~ factor1 + factor2 + (1 | subject) + (1 | item)
## link: threshold:
## clmm.one logit flexible
## clmm.int logit flexible
##
## no.par AIC logLik LR.stat df Pr(>Chisq)
## clmm.one 9 3536.0 -1759.0
## clmm.int 10 2863.5 -1421.8 674.43 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1As expected, Factor 1 contributes significantly to the model’s fit. But, if you look below, you will see that Factor 2 also contributes significantly (although the Likelihood Ratio statistic is smaller for Factor 2 than Factor 1).
anova(clmm.int,clmm.two)## Likelihood ratio tests of cumulative link models:
##
## formula:
## clmm.two as.factor(rating) ~ factor1 + (1 | subject) + (1 | item)
## clmm.int as.factor(rating) ~ factor1 + factor2 + (1 | subject) + (1 | item)
## link: threshold:
## clmm.two logit flexible
## clmm.int logit flexible
##
## no.par AIC logLik LR.stat df Pr(>Chisq)
## clmm.two 9 2970.1 -1476.0
## clmm.int 10 2863.5 -1421.8 108.53 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1That solves our problem of the graphs! But maybe you’re not convinced. After all, the box and whisker plot for Factor 2 aren’t very different…
The problem here is that box plots assume your data are continuous. This is why the y-axis is continuous and why the boxes are filled in even between the discrete numbers on the y-axis. But ordinal data isn’t continuous. There is a better way to visualise it, but it takes a little explanation:
Visualising the ratings
If we visualise the ratings as probabilities, we can see subtle differences that aren’t visible in box plots or bar charts. The x-axis is the value of \(\beta\), or the estimate of the coefficients (with the baseline set to 0). This means that if two conditions are very similar in terms of how participants used the ordinal scale to rate these stimuli, they will appear close together on the x-axis. The y-axis is the probability with which participants used each level of the scale, from 0 (extremely unlikely) to 1 (extremely likely). Crucially, this allows us to visualise how different levels vary within a condition, as well as between, using probability curves.
In this graph below, we can see the probability curves for each of the seven ratings. In particular, we can now see that the Baseline+Present condition and the Baseline+Absent condition are not rated in precisely the same manner: Baseline+Absent is more likely to receive a rating of 7 and less like to receive a rating of 5 than Baseline+Present.
# Script written by Lauren M Ackerman
# 26 February 2019
# Contact: l.m.ackerman AT gmail DOT com
# Adapted from Jalal Al-Tamimi
clmm.graph <- clmm(as.factor(rating) ~ interaction(factor1,factor2) +
(1|subject) +
(1|item),
data=dataOrd) # ordinal regression with INTERACTION BUILT INTO MAIN EFFECT## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.vlines <- c(0, # intercept, 0 = Baseline.Absent
clmm.graph$beta[[1]], # interaction(factor1, factor2) = Treatment.Absent
clmm.graph$beta[[2]], # interaction(factor1, factor2) = Baseline.Present
clmm.graph$beta[[3]]) # interaction(factor1, factor2) = Treatment.Present
xaxis <- seq(min(vlines-.5), max(vlines+.5), length.out = 100) # create 100 steps
yaxis <- rep(c(0,1),50) # fill in 0s and 1s for y-axis
colors <- c("one" = "orangered",
"two" = "orange2",
"three" = "yellow",
"four" = "green",
"five" = "turquoise",
"six" = "blue",
"seven" = "purple") # establish the colour-rating level correspondence
tibble(xaxis,yaxis) %>% # baseline tibble for plot dimensions
mutate(one= plogis(clmm.graph$Theta[1] - xaxis), # add column for rating levels (must be adapted for larger scales)
two= plogis(clmm.graph$Theta[2] - xaxis) - plogis(clmm.graph$Theta[1] - xaxis),
three=plogis(clmm.graph$Theta[3] - xaxis) - plogis(clmm.graph$Theta[2] - xaxis),
four= plogis(clmm.graph$Theta[4] - xaxis) - plogis(clmm.graph$Theta[3] - xaxis),
five= plogis(clmm.graph$Theta[5] - xaxis) - plogis(clmm.graph$Theta[4] - xaxis),
six= plogis(clmm.graph$Theta[6] - xaxis) - plogis(clmm.graph$Theta[5] - xaxis),
seven=1 - (plogis(clmm.graph$Theta[6] - xaxis))) %>%
gather(rating,probability,3:9) %>% # make long data (change right-hand number to suit the ordinal levels)
mutate(rating=factor(rating,levels=c("one","two","three","four","five","six","seven"))) %>% # make factor and relevel
ggplot(aes(x=xaxis,y=yaxis)) + # set up ggplot
geom_hline(yintercept=0,lty="dotted") + # add lower horizontal line
geom_hline(yintercept=1,lty="dotted") + # add upper horizontal line
geom_line(aes(y=probability,colour=rating),lwd=1,alpha=.8) + # add predicted curves
annotate("segment", # type of annotation = line segments
x=vlines, y=0, xend=vlines, yend=1, # add estimates
lty="solid", alpha=.75) + # visual properties of vertical lines
annotate("text", # type of annotation = text
x=vlines,y=.75, # location of labels
label=c("Baseline-Absent",
"Treatment-Absent",
"Baseline-Present",
"Treatment-Present"), # label names aligned with vlines[1:4]
angle=90,vjust=-0.2) + # visual properties of text labels
scale_x_continuous(breaks=c(min(xaxis-.5),max(xaxis+.5))) + # expand x axis horizontally
scale_y_continuous(breaks=c(0,.25,.5,.75,1)) + # expand y axis with consistent breaks
ylab("Probability") + xlab("") + ggtitle("Predicted curves") + # label plot properties
scale_colour_manual(values = colors) + # apply colours manually
theme_bw() # improve visibility with white background
However, they are still very similar in terms of their mean rating (which is why they are grouped so closely). In some cases, when this difference ends up being significant, it will not be visible in a box plot but will be visible in this sort of probability curve plot.
One could also use this method to visualise the main effect of Factor 2, in order to get an idea of whether its two levels really are similar, or whether they just have similar median values (which is what the box plot shows us).
clmm.graph <- clmm(as.factor(rating) ~ factor2 +
(1|subject) +
(1|item),
data=dataOrd) # ordinal regression with ONLY FACTOR 2
vlines <- c(0, # intercept, 0 = Absent
clmm.graph$beta[[1]]) # interaction(factor1, factor2) = Present
xaxis <- seq(min(vlines-.5), max(vlines+.5), length.out = 100) # create 100 steps
yaxis <- rep(c(0,1),50) # fill in 0s and 1s for y-axis
colors <- c("one" = "orangered",
"two" = "orange2",
"three" = "yellow",
"four" = "green",
"five" = "turquoise",
"six" = "blue",
"seven" = "purple") # establish the colour-rating level correspondence
tibble(xaxis,yaxis) %>% # baseline tibble for plot dimensions
mutate(one= plogis(clmm.graph$Theta[1] - xaxis), # add column for rating levels (must be adapted for larger scales)
two= plogis(clmm.graph$Theta[2] - xaxis) - plogis(clmm.graph$Theta[1] - xaxis),
three=plogis(clmm.graph$Theta[3] - xaxis) - plogis(clmm.graph$Theta[2] - xaxis),
four= plogis(clmm.graph$Theta[4] - xaxis) - plogis(clmm.graph$Theta[3] - xaxis),
five= plogis(clmm.graph$Theta[5] - xaxis) - plogis(clmm.graph$Theta[4] - xaxis),
six= plogis(clmm.graph$Theta[6] - xaxis) - plogis(clmm.graph$Theta[5] - xaxis),
seven=1 - (plogis(clmm.graph$Theta[6] - xaxis))) %>%
gather(rating,probability,3:9) %>% # make long data (change right-hand number to suit the ordinal levels)
mutate(rating=factor(rating,levels=c("one","two","three","four","five","six","seven"))) %>% # make factor and relevel
ggplot(aes(x=xaxis,y=yaxis)) + # set up ggplot
geom_hline(yintercept=0,lty="dotted") + # add lower horizontal line
geom_hline(yintercept=1,lty="dotted") + # add upper horizontal line
geom_line(aes(y=probability,colour=rating),lwd=1,alpha=.8) + # add predicted curves
annotate("segment", # type of annotation = line segments
x=vlines, y=0, xend=vlines, yend=1, # add estimates
lty="solid", alpha=.75) + # visual properties of vertical lines
annotate("text", # type of annotation = text
x=vlines,y=.75, # location of labels
label=c("Baseline-Absent",
"Baseline-Present"), # label names aligned with vlines[1:4]
angle=90,vjust=-0.2) + # visual properties of text labels
scale_x_continuous(breaks=c(min(xaxis-.5),max(xaxis+.5))) + # expand x axis horizontally
scale_y_continuous(breaks=c(0,.25,.5,.75,1)) + # expand y axis with consistent breaks
ylab("Probability") + xlab("") + ggtitle("Predicted curves") + # label plot properties
scale_colour_manual(values = colors) + # apply colours manually
theme_bw() # improve visibility with white background
This is an extremely informative, if somewhat boring-looking graph. In the box plot, it was hard to tell whether these two groups different because their interquartile ranges were so similar. Here, we get a clearer view of why they do indeed differ. In the Baseline+Present condition, ratings of seven (and to a lesser extent six) are much less likely than in the Baseline+Absent condition. The inverse is true for ratings of one and two (and to a lesser extent, three and four). Ratings of five appear to be approximately equally likely for both.
Some visualisations under construction
I have wondered whether it is possible to visualise the variable “sizes” of each of the levels of the ordinal scale. This is an attempt to plot them to visualise how they vary, including standard error confidence intervals. If you have suggestions for improvements, amendments, or uses, please email me at lauren [dot] ackerman [at] ncl [dot] ac [dot] uk.
clmm.graph <- clmm(as.factor(rating) ~ interaction(factor1,factor2) +
(1|subject) +
(1|item),
data=dataOrd) # ordinal regression with INTERACTION BUILT INTO MAIN EFFECT## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.levellines <- c(summary(clmm.max)$coefficient[,1][[1]],
summary(clmm.max)$coefficient[,1][[2]],
summary(clmm.max)$coefficient[,1][[3]],
summary(clmm.max)$coefficient[,1][[4]],
summary(clmm.max)$coefficient[,1][[5]],
summary(clmm.max)$coefficient[,1][[6]])
levellinesSE <- c(summary(clmm.max)$coefficient[,2][[1]],
summary(clmm.max)$coefficient[,2][[2]],
summary(clmm.max)$coefficient[,2][[3]],
summary(clmm.max)$coefficient[,2][[4]],
summary(clmm.max)$coefficient[,2][[5]],
summary(clmm.max)$coefficient[,2][[6]])
tibble(xaxis,yaxis) %>% # baseline tibble for plot dimensions
mutate(one= plogis(clmm.graph$Theta[1] - xaxis), # add column for rating levels (must be adapted for larger scales)
two= plogis(clmm.graph$Theta[2] - xaxis) - plogis(clmm.graph$Theta[1] - xaxis),
three=plogis(clmm.graph$Theta[3] - xaxis) - plogis(clmm.graph$Theta[2] - xaxis),
four= plogis(clmm.graph$Theta[4] - xaxis) - plogis(clmm.graph$Theta[3] - xaxis),
five= plogis(clmm.graph$Theta[5] - xaxis) - plogis(clmm.graph$Theta[4] - xaxis),
six= plogis(clmm.graph$Theta[6] - xaxis) - plogis(clmm.graph$Theta[5] - xaxis),
seven=1 - (plogis(clmm.graph$Theta[6] - xaxis))) %>%
gather(rating,probability,3:9) %>% # make long data (change right-hand number to suit the ordinal levels)
mutate(rating=factor(rating,levels=c("one","two","three","four","five","six","seven"))) %>% # make factor and relevel
ggplot(aes(x=xaxis,y=yaxis)) + # set up ggplot
geom_hline(yintercept=0,lty="dotted") + # add lower horizontal line
geom_hline(yintercept=1,lty="dotted") + # add upper horizontal line
#geom_line(aes(y=probability,colour=rating),lwd=1,alpha=.8) + # add predicted curves
annotate("segment", # type of annotation = line segments
x=levellines, y=0, xend=levellines, yend=1, # add estimates
lty="solid") +
annotate("rect", # type of annotation = rectangle
xmin=levellines-levellinesSE, ymin=0, xmax=levellines+levellinesSE, ymax=1, # add estimates
lty="dotted", alpha=.5,fill="grey") +
annotate("text", # type of annotation = text
x=levellines,y=.75, # location of labels
label=c("1|2","2|3","3|4","4|5","5|6","6|7")) + # label names aligned with vlines[1:4]
xlab("Thresholds Beta value and SE") +
ylab("") +
theme_bw()
One interesting insight this plot provides is that the lower levels of the ordinal scale are both smaller and more variable (or less certain).
I would like to be able to plot confidence intervals around the predicted curves. This is my attempt to do so:
clmm.graph <- clmm(as.factor(rating) ~ interaction(factor1,factor2) +
(1|subject) +
(1|item),
data=dataOrd) # ordinal regression with INTERACTION BUILT INTO MAIN EFFECT## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.vlines <- c(0, # intercept, 0 = Baseline.Absent
clmm.graph$beta[[1]], # interaction(factor1, factor2) = Treatment.Absent
clmm.graph$beta[[2]], # interaction(factor1, factor2) = Baseline.Present
clmm.graph$beta[[3]]) # interaction(factor1, factor2) = Treatment.Present
xaxis <- seq(min(vlines-.5), max(vlines+.5), length.out = 100) # create 100 steps
yaxis <- rep(c(0,1),50) # fill in 0s and 1s for y-axis
colors <- c("one" = "orangered",
"two" = "orange2",
"three" = "yellow",
"four" = "green",
"five" = "turquoise",
"six" = "blue",
"seven" = "purple") # establish the colour-rating level correspondence
tibble(xaxis,yaxis) %>% # baseline tibble for plot dimensions
mutate(one.base= plogis(clmm.graph$Theta[1] - xaxis), # add column for rating levels (must be adapted for larger scales)
two.base= plogis(clmm.graph$Theta[2] - xaxis) -
plogis(clmm.graph$Theta[1] - xaxis),
three.base=plogis(clmm.graph$Theta[3] - xaxis) -
plogis(clmm.graph$Theta[2] - xaxis),
four.base= plogis(clmm.graph$Theta[4] - xaxis) -
plogis(clmm.graph$Theta[3] - xaxis),
five.base= plogis(clmm.graph$Theta[5] - xaxis) -
plogis(clmm.graph$Theta[4] - xaxis),
six.base= plogis(clmm.graph$Theta[6] - xaxis) -
plogis(clmm.graph$Theta[5] - xaxis),
seven.base=1 - (plogis(clmm.graph$Theta[6] - xaxis))) %>%
mutate(one.add= plogis(clmm.graph$Theta[1]+summary(clmm.max)$coefficient[,2][[1]] - xaxis), # add column for upper SE CI (must be adapted for larger scales)
two.add= plogis(clmm.graph$Theta[2]+summary(clmm.max)$coefficient[,2][[2]] - xaxis) -
plogis(clmm.graph$Theta[1]+summary(clmm.max)$coefficient[,2][[1]] - xaxis),
three.add=plogis(clmm.graph$Theta[3]+summary(clmm.max)$coefficient[,2][[3]] - xaxis) -
plogis(clmm.graph$Theta[2]+summary(clmm.max)$coefficient[,2][[2]] - xaxis),
four.add= plogis(clmm.graph$Theta[4]+summary(clmm.max)$coefficient[,2][[4]] - xaxis) -
plogis(clmm.graph$Theta[3]+summary(clmm.max)$coefficient[,2][[3]] - xaxis),
five.add= plogis(clmm.graph$Theta[5]+summary(clmm.max)$coefficient[,2][[5]] - xaxis) -
plogis(clmm.graph$Theta[4]+summary(clmm.max)$coefficient[,2][[4]] - xaxis),
six.add= plogis(clmm.graph$Theta[6]+summary(clmm.max)$coefficient[,2][[6]] - xaxis) -
plogis(clmm.graph$Theta[5]+summary(clmm.max)$coefficient[,2][[5]] - xaxis),
seven.add=1 - (plogis(clmm.graph$Theta[6]+summary(clmm.max)$coefficient[,2][[6]] - xaxis))) %>%
mutate(one.sub= plogis(clmm.graph$Theta[1]-summary(clmm.max)$coefficient[,2][[1]] - xaxis), # add column for lower SE CI (must be adapted for larger scales)
two.sub= plogis(clmm.graph$Theta[2]-summary(clmm.max)$coefficient[,2][[2]] - xaxis) -
plogis(clmm.graph$Theta[1]-summary(clmm.max)$coefficient[,2][[1]] - xaxis),
three.sub=plogis(clmm.graph$Theta[3]-summary(clmm.max)$coefficient[,2][[3]] - xaxis) -
plogis(clmm.graph$Theta[2]-summary(clmm.max)$coefficient[,2][[2]] - xaxis),
four.sub= plogis(clmm.graph$Theta[4]-summary(clmm.max)$coefficient[,2][[4]] - xaxis) -
plogis(clmm.graph$Theta[3]-summary(clmm.max)$coefficient[,2][[3]] - xaxis),
five.sub= plogis(clmm.graph$Theta[5]-summary(clmm.max)$coefficient[,2][[5]] - xaxis) -
plogis(clmm.graph$Theta[4]-summary(clmm.max)$coefficient[,2][[4]] - xaxis),
six.sub= plogis(clmm.graph$Theta[6]-summary(clmm.max)$coefficient[,2][[6]] - xaxis) -
plogis(clmm.graph$Theta[5]-summary(clmm.max)$coefficient[,2][[5]] - xaxis),
seven.sub=1 - (plogis(clmm.graph$Theta[6]-summary(clmm.max)$coefficient[,2][[6]] - xaxis))) %>%
gather(level,probability,3:23) %>% # make long data (change right-hand number to suit the ordinal levels)
separate(level,c("rating","type")) %>% # tidy values in "rating" column for next step
spread(type,probability) %>% # spread so it's slightly wider data, keeping baseline and SE CI (upper and lower) values seperate
mutate(rating=factor(rating,levels=c("one","two","three","four","five","six","seven"))) %>% # make factor and relevel
ggplot(aes(x=xaxis,y=yaxis)) + # set up ggplot
geom_hline(yintercept=0,lty="dotted") + # add lower horizontal line
geom_hline(yintercept=1,lty="dotted") + # add upper horizontal line
geom_ribbon(aes(ymin=sub,ymax=add,fill=rating),lwd=1,alpha=.2) + # add predicted curves confidence intervals *** NEEDS WORK **
geom_line(aes(y=base,colour=rating),lwd=1,alpha=.8) + # add predicted curves
annotate("segment", # type of annotation = line segments
x=vlines, y=0, xend=vlines, yend=1, # add estimates
lty="solid", alpha=.75) + # visual properties of vertical lines
annotate("text", # type of annotation = text
x=vlines,y=.75, # location of labels
label=c("Baseline-Absent",
"Treatment-Absent",
"Baseline-Present",
"Treatment-Present"), # label names aligned with vlines[1:4]
angle=90,vjust=-0.2) + # visual properties of text labels
scale_x_continuous(breaks=c(min(xaxis-.5),max(xaxis+.5))) + # expand x axis horizontally
scale_y_continuous(breaks=c(0,.25,.5,.75,1)) + # expand y axis with consistent breaks
ylab("Probability") + xlab("") + ggtitle("Predicted curves") + # label plot properties
scale_colour_manual(values = colors) + # apply colours manually
theme_bw() # improve visibility with white background
I think this is neither necessarily the correct nor the easiest way to do this, and I cannot recommend using this code to produce professional visualisations of the predicted curves’ standard error confidence intervals. I am providing this code so that others may build on it and feed back to me so that I can update it and complete the tutorial properly. ¡Buena suerte!
Summary
Ordinal regressions are more appropriate for ordinal data than linear regressions are for a number of reasons. You don’t need to transform your data, you can account for the levels not being “sized” linearly and you can account for subjects treating the individual levels differently than each other.
If you analyse an ordinal regression, however, it is also important to visualise your data in an appropriate way. Box plots (and similarly, violin plots) are not optimal because they assume continuous data. Ordinal data is discrete and its levels are independent, so it makes more sense to visualise the probability of each level being selected, rather than the mean or median values for each condition.