检验 R 中的比例赔率假设

Question

我在R中使用响应变量，该变量是学生在特定课程中获得的字母等级。 反应是有序的，在我看来，在逻辑上似乎是相称的。我的理解是，我需要先测试它是成比例的，然后才能使用 polr(( 而不是 multinom((。

对于我的一个数据课程，我像这样"测试"了比例：

M1 <- logLik(polrModel)  #'log Lik.' -1748.180691 (df=8)
M2 <- logLik(multinomModel)  #'log Lik.' -1734.775727 (df=20)
G <- -2*(M1$1 - M2$2)) #I used a block bracket here in the real code
# 26.8099283
pchisq(G,12,lower.tail = FALSE) #DF is #of predictors
#0.008228890393     #THIS P-VAL TELLS ME TO REJECT PROPORTIONAL

对于测试比例赔率假设的第二种方法，我还运行了两个 vglm 模型，一个使用 family=cumulative(parallel =TRUE) 另一个使用 family=cumulative(parallel =FALSE) 。 然后，我用模型偏差的差异和剩余自由度的差异进行了pchisq()测试。

这两种方式都值得尊敬吗？ 如果没有，我很乐意帮助实际编码，以确定是接受还是拒绝比例赔率假设！

除了上述两个测试之外，我还分别绘制了每个预测因子的累积概率图。 我读到我希望这些线是平行的。 我不明白的是，polr()你的输出是每个自变量(系数(的单个斜率，然后是一个特定的截距，具体取决于你正在使用的累积概率(例如：P(Y<=A(、P(Y<=B(等(。 那么，如果每个方程的斜率系数都相同，那么这些线怎么可能不平行呢？

我在克里斯·比尔德(Chris Bilder(的YouTube课上学习了我的基本知识;他在第42分钟谈到了并行图。

任何帮助不胜感激！谢谢！

Answer 1

你的方法基本上是正确的。 我有以下代码，灵感来自 Fox 的"应用回归的 R 和 S-PLUS 伴侣"。第 5 章：拟合广义线性模型。第155-189页。使用代码时，请引用书籍章节。本章也有关于绘图的部分。

library(car)
library(nnet)
library(xlsx)
library(MASS)
options(warn=1)
options(digits = 3)
#
Trial <- read.xlsx("Trial.xls", "Sheet 1")
# Set up an out file structure
sink("Testing_adequacy_of_Prop_odds.txt")
# Trial$Outcome is assessed on a six point scale 0-5
schtyp_M_M.f <- factor(Trial$Outcome, labels = c("M0", "M1", "M2", "M3", "M4", "M5"))
#
cat("Multinomial logistic regression n")
# Assign takes on a value of 1 (Treatment) or 0 (Control) 
mod.multinom <-multinom(schtyp_M_M.f~Assign, data = Trial)
print(summary(mod.multinom, cor=F, Wald=T))
x1<-logLik(mod.multinom)
cat("Degrees of freedom Multinomial logistic regression n")
print(df_of_multinom_model <- attributes(x1)$df)
cat("Proportional odds logistic regressionn")
mod.polr <- polr(schtyp_M_M.f ~ Assign, data=Trial)
print(summary(mod.polr))
x2<-logLik(mod.polr)
cat("Degrees of freedom Proportional Odds Logistic Regression n")
print(df_of_polr_model <- attributes(x2)$df)
cat("Answering the question: Is proportional odds model assumption violatedn")
cat("P value for difference in AIC between POLR and Multinomial Logit modeln")
# abs since the values could be negative. That is negative difference of degrees of freedom would produce p=NaN
print(1-pchisq(abs(mod.polr$deviance-mod.multinom$deviance),   abs(df_of_multinom_model-df_of_polr_model)))
sink()