在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)。然后，我使用模型差异和残余自由度差异运行了测试。

这两种方法是否可信？如果不是，请帮助我编写实际的代码来确定是否接受或拒绝比例几率假设！

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

我通过观看Chris Bilder的YouTube课程学习了基础知识；他在第42分钟这里谈到了平行图形。

感谢您的帮助！

- Nameless

---

这个问题更像是一个统计学问题，而不是一个编程问题。你应该在[stats.se]寻求统计建议，而不是在Stack Overflow上提问。 - MrFlick

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 regression\n")
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 violated\n")
cat("P value for difference in AIC between POLR and Multinomial Logit model\n")
# 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()

- abcihep