I fitted a proportional odds cumulative logit model on ordinal data using MASS
's polr
function using (in this case on data giving the preference for different kinds of cheese) :
data=read.csv("https://www.dropbox.com/s/psj74dx8ohnrdlp/cheese.csv?dl=1")
data$response=factor(data$response, ordered=T) # make response into ordered factor
head(data)
cheese response count
1 A 1 0
2 A 2 0
3 A 3 1
4 A 4 7
5 A 5 8
6 A 6 8
library(MASS)
fit=polr(response ~ cheese, weights=count, data=data, Hess=TRUE, method="logistic")
To plot the predictions of the model I made an effect plot using
library(effects)
library(colorRamps)
plot(allEffects(fit),ylab="Response",type="probability",style="stacked",colors=colorRampPalette(c("white","red"))(9))
I was wondering though if from the predicted means reported by the effects
package one could also plot something like the mean preference for each kind of cheese together with the 95% conf intervals on this?
EDIT: originally I also asked about how to obtain Tukey posthoc tests, but in the meantime I found that those can be obtained using
library(multcomp)
summary(glht(fit, mcp(cheese = "Tukey")))
or using package lsmeans
as
summary(lsmeans(fit, pairwise ~ cheese, adjust="tukey", mode = "linear.predictor"),type="response")
Russ Lenth just kindly pointed out to me that the mean preference and 95% confidence intervals can be obtained simply with lsmeans
with option mode="mean"
(?models
) (the same also applies to a clm
or clmm
model fit using package ordinal
):
df=summary(lsmeans(fit, pairwise ~ cheese, mode = "mean"),type="response")$lsmeans
cheese mean.class SE df asymp.LCL asymp.UCL
A 6.272828 0.1963144 NA 5.888058 6.657597
B 3.494899 0.2116926 NA 3.079989 3.909809
C 4.879440 0.2006915 NA 4.486091 5.272788
D 7.422159 0.1654718 NA 7.097840 7.746478
which gives me the plot I was looking for :
library(ggplot2)
library(ggthemes)
ggplot(df, aes(cheese, mean.class)) + geom_bar(stat="identity", position="dodge", fill="steelblue", width=0.6) +
geom_errorbar(aes(ymin=asymp.LCL, ymax=asymp.UCL), width=0.15, position=position_dodge(width=0.9)) +
theme_few(base_size=18) + xlab("Type of cheese") + ylab("Mean preference") +
coord_cartesian(ylim = c(1, 9)) + scale_y_continuous(breaks=1:9)
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With