Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

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

enter image description here

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")
like image 844
Tom Wenseleers Avatar asked Oct 24 '15 09:10

Tom Wenseleers


1 Answers

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)

enter image description here

like image 185
Tom Wenseleers Avatar answered Nov 03 '22 09:11

Tom Wenseleers