Das Überprüfen der Proportional-Odds-Annahme gilt für eine ordinale logistische Regression mit der Polr-Funktion

Ich habe die Funktion 'polr' im MASS-Paket verwendet, um eine ordinale logistische Regression für eine ordinale kategoriale Antwortvariable mit 15 kontinuierlichen erklärenden Variablen auszuführen.

Ich habe den Code (siehe unten) verwendet, um zu überprüfen, ob mein Modell die Annahme der proportionalen Gewinnchancen gemäß den Empfehlungen im UCLA-Leitfaden erfüllt . Ich bin jedoch ein wenig besorgt über die Ausgabe, die impliziert, dass die Koeffizienten über verschiedene Schnittpunkte nicht nur ähnlich sind, sondern genau gleich sind (siehe Grafik unten).

FGV1b <- data.frame(FG1_val_cat=factor(FGV1b[,"FG1_val_cat"]), 
                    scale(FGV1[,c("X","Y","Slope","Ele","Aspect","Prox_to_for_FG", 
                          "Prox_to_for_mL", "Prox_to_nat_border", "Prox_to_village", 
                          "Prox_to_roads", "Prox_to_rivers", "Prox_to_waterFG", 
                          "Prox_to_watermL", "Prox_to_core", "Prox_to_NR", "PCA1", 
                          "PCA2", "PCA3")]))
b     <- polr(FG1_val_cat ~ X + Y + Slope + Ele + Aspect + Prox_to_for_FG + 
                            Prox_to_for_mL + Prox_to_nat_border + Prox_to_village + 
                            Prox_to_roads + Prox_to_rivers + Prox_to_waterFG + 
                            Prox_to_watermL + Prox_to_core + Prox_to_NR, 
              data=FGV1b, Hess=TRUE)

Eine Zusammenfassung des Modells anzeigen:

summary(b)
(ctableb <- coef(summary(b)))
q        <- pnorm(abs(ctableb[, "t value"]), lower.tail=FALSE) * 2
(ctableb <- cbind(ctableb, "p value"=q))

Und jetzt können wir uns die Konfidenzintervalle für die Parameterschätzungen ansehen:

(cib <- confint(b)) 
confint.default(b)

Diese Ergebnisse sind jedoch immer noch schwer zu interpretieren. Lassen Sie uns die Koeffizienten in Quotenverhältnisse umwandeln

exp(cbind(OR=coef(b), cib))

Überprüfung der Annahme. Der folgende Code schätzt also die zu grafischen Werte. Zunächst werden die Logit-Transformationen der Wahrscheinlichkeiten gezeigt, die größer oder gleich jedem Wert der Zielvariablen sind

FG1_val_cat <- as.numeric(FG1_val_cat)
sf <- function(y) {
  c('VC>=1' = qlogis(mean(FG1_val_cat >= 1)),
    'VC>=2' = qlogis(mean(FG1_val_cat >= 2)),
    'VC>=3' = qlogis(mean(FG1_val_cat >= 3)),
    'VC>=4' = qlogis(mean(FG1_val_cat >= 4)),
    'VC>=5' = qlogis(mean(FG1_val_cat >= 5)),
    'VC>=6' = qlogis(mean(FG1_val_cat >= 6)),
    'VC>=7' = qlogis(mean(FG1_val_cat >= 7)),
    'VC>=8' = qlogis(mean(FG1_val_cat >= 8)))
}
(t <- with(FGV1b, summary(as.numeric(FG1_val_cat) ~ X + Y + Slope + Ele + Aspect + 
                             Prox_to_for_FG + Prox_to_for_mL + Prox_to_nat_border + 
                             Prox_to_village + Prox_to_roads + Prox_to_rivers + 
                             Prox_to_waterFG + Prox_to_watermL + Prox_to_core + 
                             Prox_to_NR, fun=sf)))

Die obige Tabelle zeigt die (linearen) vorhergesagten Werte, die wir erhalten würden, wenn wir unsere abhängige Variable nacheinander ohne die Annahme paralleler Steigungen von unseren Prädiktorvariablen zurückführen würden. Jetzt können wir eine Reihe von binären logistischen Regressionen mit unterschiedlichen Schnittpunkten für die abhängige Variable ausführen, um die Gleichheit der Koeffizienten über die Schnittpunkte hinweg zu überprüfen

par(mfrow=c(1,1))
plot(t, which=1:8, pch=1:8, xlab='logit', main=' ', xlim=range(s[,7:8]))

Entschuldigung, dass ich kein Statistikexperte bin und hier vielleicht etwas Offensichtliches fehlt. Ich habe jedoch lange Zeit versucht, herauszufinden, ob es ein Problem gibt, wie ich die Modellannahme getestet habe, und auch andere Möglichkeiten zu finden, um dieselbe Art von Modell auszuführen.

Zum Beispiel habe ich in vielen Hilfe-Mailinglisten gelesen, dass andere die vglm-Funktion (im VGAM-Paket) und die lrm-Funktion (im rms-Paket) verwenden (siehe hier: Proportional-Odds-Annahme in der ordinalen logistischen Regression in R mit den Paketen VGAM und rms ). Ich habe versucht, die gleichen Modelle zu verwenden, stoße aber ständig auf Warnungen und Fehler.

Wenn ich beispielsweise versuche, das vglm-Modell mit dem Argument 'parallel = FALSE' anzupassen (wie im vorherigen Link erwähnt, ist dies wichtig, um die Annahme der proportionalen Gewinnchancen zu testen), tritt der folgende Fehler auf:

Fehler in lm.fit (X.vlm, y = z.vlm, ...): NA / NaN / Inf in 'y'
Zusätzlich: Warnmeldung:
In Deviance.categorical.data.vgam (mu = mu, y = y, w = w, Residuen = Residuen ,: angepasste Werte nahe 0 oder 1

Ich möchte bitte fragen, ob es jemanden gibt, der verstehen und mir erklären kann, warum das Diagramm, das ich oben erstellt habe, so aussieht, wie es aussieht. Wenn dies tatsächlich bedeutet, dass etwas nicht stimmt, können Sie mir bitte helfen, einen Weg zu finden, um die Proportional-Odds-Annahme zu testen, wenn Sie nur die polr-Funktion verwenden. Oder wenn dies einfach nicht möglich ist, werde ich versuchen, die vglm-Funktion zu verwenden, würde dann aber Hilfe benötigen, um zu erklären, warum ich den oben angegebenen Fehler immer wieder erhalte.

HINWEIS: Als Hintergrund gibt es hier 1000 Datenpunkte, die tatsächlich Standortpunkte in einem Untersuchungsgebiet sind. Ich suche nach Beziehungen zwischen der kategorialen Antwortvariablen und diesen 15 erklärenden Variablen. Alle diese 15 erklärenden Variablen sind räumliche Merkmale (z. B. Höhe, xy-Koordinaten, Nähe zum Wald usw.). Die 1000 Datenpunkte wurden zufällig mithilfe eines GIS zugewiesen, aber ich habe einen geschichteten Stichprobenansatz gewählt. Ich habe dafür gesorgt, dass 125 Punkte in jeder der 8 verschiedenen kategorialen Antwortstufen zufällig ausgewählt wurden. Ich hoffe, diese Informationen sind auch hilfreich.

Die abhängige Variable hat 8 geordnete Ebenen. Um die Proportional-Odds-Annahme zu testen, sollten Sie in der Grafik 8 verschiedene Symbole für jede unabhängige Variable sehen. Sie sehen nur 2 Symbole für jede unabhängige Variable, wahrscheinlich weil Sie ein zu kurzes Intervall für die Werte der x-Achse gewählt haben. Wenn meine Vermutung richtig ist, müssen Sie nur ein breiteres Intervall für die Werte der x-Achse verwenden. Versuchen Sie diesen Code:

par(mfrow=c(1,1))
plot(t, which=1:8, pch=1:8, xlab='logit', main=' ', xlim=range(s[,3:9]))

Also fand ich das durch googeln und ich denke, eine Antwort könnte aus diesem Grund immer noch nützlich sein. Ich denke, der Fehler liegt darin

sf <- function(y) {
  c('VC>=1' = qlogis(mean(FG1_val_cat >= 1)),
    'VC>=2' = qlogis(mean(FG1_val_cat >= 2)),
    'VC>=3' = qlogis(mean(FG1_val_cat >= 3)),
    'VC>=4' = qlogis(mean(FG1_val_cat >= 4)),
    'VC>=5' = qlogis(mean(FG1_val_cat >= 5)),
    'VC>=6' = qlogis(mean(FG1_val_cat >= 6)),
    'VC>=7' = qlogis(mean(FG1_val_cat >= 7)),
    'VC>=8' = qlogis(mean(FG1_val_cat >= 8)))
}

wo du FG1_val_cateher verwendest als y. Verwenden Sie das Beispiel aus Harrells Regressionsmodellierungsstrategien:

library(Hmisc)
getHdata(support)
support <- support[complete.cases(support[, c("sfdm2", "adlsc", "sex", "age", "meanbp")]), ]
sfdm <- as.integer (support$sfdm2 ) - 1

sf1 <- function (y) {
  c(' Y ≥ 1 ' = qlogis (mean(sfdm >= 1)), 
    ' Y ≥ 2 ' = qlogis (mean(sfdm >= 2)),
    ' Y ≥ 3 ' = qlogis (mean(sfdm >= 3))
  )
}

sf2 <- function (y) {
  c(' Y ≥ 1 ' = qlogis (mean(y >= 1)), 
    ' Y ≥ 2 ' = qlogis (mean(y >= 2)),
    ' Y ≥ 3 ' = qlogis (mean(y >= 3))
  )
}

s1 <- summary(sfdm ~ adlsc + sex + age + meanbp, fun=sf1,
              data = support)
s2 <- summary(sfdm ~ adlsc + sex + age + meanbp, fun=sf2,
              data = support)  

plot(s1, which =1:3, pch =1:3, xlab = ' logit ', main = ' ', width.factor = 1.4, cex.lab = 0.75)

plot(s2, which =1:3, pch =1:3, xlab = ' logit ', main = ' ', width.factor = 1.4, cex.lab = 0.75)

vs.