title: R: Plotting predictions of MASS polr ordinal model
body: <p>I fitted a proportional odds cumulative logit model on ordinal data using <code>MASS</code>'s <code>polr</code> function using (in this case on data giving the preference for different kinds of cheese) :</p>
  
  <pre><code>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")
  </code></pre>
  
  <p>To plot the predictions of the model I made an effect plot using</p>
  
  <pre><code>library(effects)
  library(colorRamps)
  plot(allEffects(fit),ylab="Response",type="probability",style="stacked",colors=colorRampPalette(c("white","red"))(9))
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/HhetJ.png" rel="noreferrer"><img src="https://i.stack.imgur.com/HhetJ.png" alt="enter image description here"></a></p>
  
  <p>I was wondering though if from the predicted means reported by  the <code>effects</code> package one could also plot something like the mean preference for each kind of cheese together with the 95% conf intervals on this?</p>
  
  <p>EDIT: originally I also asked about how to obtain Tukey posthoc tests, but in the meantime I found that those can be obtained using</p>
  
  <pre><code>library(multcomp)
  summary(glht(fit, mcp(cheese = "Tukey")))
  </code></pre>
  
  <p>or using package <code>lsmeans</code> as</p>
  
  <pre><code>summary(lsmeans(fit, pairwise ~ cheese, adjust="tukey", mode = "linear.predictor"),type="response")
  </code></pre>