15

MASSの関数を使用して、序数データに比例オッズ累積ロジット モデルを当てはめましたpolr(この場合、さまざまな種類のチーズを優先するデータに) :

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% conf 間隔と共にプロットできるかどうか疑問に思っていましたか?

編集: もともと、テューキー ポストホック テストを取得する方法についても尋ねましたが、その間に、それらを使用して取得できることがわかりました

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

またはパッケージlsmeansを次のように使用する

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

1 に答える 1

2

Russ Lenth は親切にも、平均選好と 95% 信頼区間は( ) を使用lsmeansして簡単に取得できることを指摘してくれました (同じことが、またはpackage を使用したモデル フィットにも適用されます)。option mode="mean"?modelsclmclmmordinal

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)

ここに画像の説明を入力

于 2016-02-29T15:39:53.870 に答える