R：绘制MASS polr序数模型的预测

Question

我使用 MASS 的 polr 函数，对有序数据进行了比例几率累积logit模型拟合（例如，在描述不同种类奶酪喜好的数据上）。

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

为了绘制模型预测结果，我使用了一个效果图来展示。

library(effects)
library(colorRamps)
plot(allEffects(fit),ylab="Response",type="probability",style="stacked",colors=colorRampPalette(c("white","red"))(9))

我想知道是否可以通过 effects 软件包报告的预测均值来绘制每种奶酪的平均喜好程度以及其中的95%置信区间。

编辑：最初我也询问了如何获得Tukey事后检验，但与此同时，我发现可以使用下列方法获取：

library(multcomp)
summary(glht(fit, mcp(cheese = "Tukey")))

或者使用包lsmeans：

summary(lsmeans(fit, pairwise ~ cheese, adjust="tukey", mode = "linear.predictor"),type="response")

- Tom Wenseleers

---

有趣的问题。我假设（就像你一样）问题出现在你创建预测概率之后再取平均值。请参见这里和这里以获取更多关于此问题的信息。另外，由于有9个类别，我建议对响应变量进行OLS回归，这几乎可以得到相同的平均类别点估计值，并且具有合理的标准误差。但这是一个有趣的问题。 - Felix

是的，我认为这与在累积logit比例尺上进行平均与在最终反变换比例尺上进行平均有关。因此，基本上我想知道如何在链接比例尺上进行平均，然后将其反变换回原始序数比例尺。我知道对于9个类别，我也可以只使用OLS，但我也想为较少的类别（例如3或4）提供一般解决方案。 - Tom Wenseleers

“动态柱状图（那些条形图）只是糟糕的统计数据。你从中获得的见解并不比摘要统计表中的加权均值更多。由于这仅仅是一个摘要统计图，你会失去所有用于生成该图的数据。绘图应该显示数据而不是摘要统计数据。我认为这解决了你的问题，因为你一开始就不应该这样做。” - rawr

我的问题是关于如何正确计算我的wmeans表，而不是如何最好地显示它...我非常清楚那些反对条形图的人，说实话，我从来没有完全理解过，特别是在这种情况下，我将所有内容都显示在完整的响应比例上... - Tom Wenseleers

主要问题在于你试图基于需要正态性的假设来总结非正态数据。你可以像你建议的那样，在转换后的数据上创建置信区间，然后再反向转换。另一种选择是简单地使用非参数摘要。例如，也许你的误差条可以是第一和第三四分位数。 - Sam Dickson

Ha Russ Lenth刚刚给我发了一条消息，指出如何做到这一点，请见下文！ - Tom Wenseleers

Answer 1

Russ Lenth刚刚友好地指出，可以使用lsmeans和option mode="mean"（?models）轻松获得平均偏好和95％置信区间（对于使用ordinal包拟合的clm或clmm模型同样适用）：

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

这为我提供了我所寻找的情节：

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)

- Tom Wenseleers