R: Plotting predictions of MASS polr ordinal model

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

1 Answers

Answers 1

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 Enjoyed This, Take 5 Seconds To Share It

• Facebook
• Twitter
• Google+
• StumbleUpon
• Digg
• Delicious
• LinkedIn
• Reddit
• Technorati