Testen der proportionalen Quoten Annahme in R

Question

Ich arbeite in R mit einer Antwortvariablen, die die Buchstabenklasse ist, die der Schüler in einem bestimmten Kurs erhielt. Die Antwort ist ordinal und scheint meiner Meinung nach logisch proportional zu sein. Mein Verständnis ist, dass ich testen muss, dass es proportional ist, bevor ich polr() anstelle von multinom() verwenden kann.

für einen meiner Kurse von Daten, ich die Proportionalität wie diese „getestet“:

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 

Für einen zweiten Weg, um die proportional odds Annahme zu testen, lief ich auch zwei vglm Modelle, eines mit family=cumulative(parallel =TRUE) der andere mit family=cumulative(parallel =FALSE). Ich führte dann einen pchisq() Test mit dem Unterschied der Abweichungen der Modelle und der Unterschiede der restlichen Freiheitsgrade durch.

Ist einer dieser Wege respektabel? Wenn nicht, würde ich gerne Hilfe mit der tatsächlichen Codierung für die Bestimmung, ob die Annahme der proportionalen Quoten zu akzeptieren oder abzulehnen!

Zusätzlich zu den beiden obigen Tests habe ich meine kumulativen Wahrscheinlichkeiten einzeln für jeden der Prädiktoren grafisch dargestellt. Ich habe gelesen, dass ich möchte, dass diese Zeilen parallel sind. Was ich nicht verstehe, ist, mit polr() Ihre Ausgabe ist eine einzige Steigung für jede unabhängige Variable (ein Koeffizient) und dann ein bestimmter Abschnitt abhängig davon, mit welcher kumulativen Wahrscheinlichkeit Sie arbeiten (zB: P (Y < = A), P (Y < = B) usw.). Wenn also Ihre Neigungskoeffizienten für jede der Gleichungen gleich sind, wie könnten dann die Linien nicht parallel sein?

Ich habe die Grundlagen meines Wissens in Chris Bilder's YouTube-Klasse gelernt; er spricht über die parallelen Graphen here at minute 42.

Jede Hilfe wird geschätzt! Vielen Dank!

Answer 1

Ihre Vorgehensweise ist im Wesentlichen richtig. Ich habe den folgenden Code von Fox "An R und S-PLUS Begleiter Applied Regression" inspiriert. Kapitel 5: Anpassen von linearen Modellen. Seiten 155-189. Bitte geben Sie das Buchkapitel an, wenn Sie den Code verwenden. Dieses Kapitel enthält auch einen Abschnitt zur grafischen Darstellung.

rm(list = ls()) 
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()