Wie kann man die Ausgabe von Rs polr-Funktion (geordnete logistische Regression) verstehen?

Ich bin neu in R, bestellt logistische Regression und polr.

Der Abschnitt "Beispiele" unten auf der Hilfeseite für polr (der ein logistisches oder Probit-Regressionsmodell an eine geordnete Faktorantwort anpasst ) zeigt

options(contrasts = c("contr.treatment", "contr.poly"))
house.plr <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)
pr <- profile(house.plr)
plot(pr)
pairs(pr)

Welche Informationen prenthält? Die Hilfeseite zum Profil ist allgemein gehalten und enthält keine Anleitung für polr.
Was plot(pr)zeigt? Ich sehe sechs Grafiken. Jede hat eine numerische X-Achse, obwohl die Bezeichnung eine Indikatorvariable ist (sieht aus wie eine Eingabevariable, die ein Indikator für einen Ordnungswert ist). Dann ist die Y-Achse "Tau", was völlig unerklärt ist.
Was pairs(pr)zeigt? Es sieht aus wie ein Diagramm für jedes Paar von Eingabevariablen, aber ich sehe auch hier keine Erklärung für die X- oder Y-Achse.
Wie kann man verstehen, ob das Modell gut passt? summary(house.plr)zeigt Residual Deviance 3479.149 und AIC (Akaike Information Criterion?) von 3495.149. Ist das gut? Was ist ein gutes absolutes Maß für den Fall, dass diese nur als relative Maße (dh zum Vergleich mit einer anderen Modellanpassung) nützlich sind? Ist die Restabweichung ungefähr im Chi-Quadrat verteilt? Kann man "% korrekt vorhergesagt" für die Originaldaten oder eine Kreuzvalidierung verwenden? Wie geht das am einfachsten?
Wie wendet man anovadieses Modell an und wie interpretiert man es? In den Dokumenten heißt es: "Es gibt Methoden für die Standardfunktionen zur Modellanpassung, einschließlich" Vorhersagen "," Zusammenfassung "," VCOV "und" ANOVA "." Laufen anova(house.plr)führt jedoch zuanova is not implemented for a single "polr" object
Wie interpretiert man die t-Werte für jeden Koeffizienten? Im Gegensatz zu einigen Modellanpassungen gibt es hier keine P-Werte.

Mir ist klar, dass dies eine Menge Fragen sind, aber es ist für mich sinnvoll, als ein Bündel zu fragen ("Wie verwende ich dieses Ding?") Und nicht als 7 verschiedene Fragen. Alle Informationen geschätzt.

r logistic dfrankow
quelle

@dfrankow Etwas grobe und sicherlich sehr partielle Hilfe für Ihre ersten beiden Fragen, aber methods("profile")Sie erhalten die (in diesem Fall S3) Methoden, die einem R- profileObjekt zugeordnet sind, dann sehen Sie polr, dass es eine dedizierte Methode für Ergebnisse gibt, die Sie durchsuchen können Online durch Eingabe getAnywhere("profile.polr")an der R-Eingabeaufforderung.

chl

Vielen Dank! Der Quellcode ist gut. Erklärung wäre noch besser. :)

Dfrankow

Jemand hat mich auf "Modern Applied Statistics with S" von Venables und Ripley hingewiesen. Abschnitt 7.3 enthält "Ein Beispiel für eine Vier-Wege-Frequenztabelle", die dieses Hausmodell ausführlich behandelt. Reading ..

dfrankow

Tatsächlich ist der Abschnitt "ein proportionales Gewinnchancenmodell"

dfrankow

Antworten:

Ich würde vorschlagen, dass Sie sich Bücher zur kategorialen Datenanalyse ansehen (vgl. Alan Agrestis Categorical Data Analysis, 2002), um die geordnete logistische Regression besser zu erklären und zu verstehen . Alle Fragen, die Sie stellen, werden im Grunde durch ein paar Kapitel in solchen Büchern beantwortet. Wenn Sie nur an Rverwandten Beispielen interessiert sind , ist die Erweiterung linearer Modelle in R von Julian Faraway (CRC Press, 2008) eine hervorragende Referenz.

Bevor ich Ihre Fragen beantworte, ist die geordnete logistische Regression ein Fall von multinomialen Logit-Modellen, in denen die Kategorien geordnet sind. Angenommen , wir haben geordnete Kategorien und dass für einzelne , mit Ordnungs Antwort , für . Bei einer geordneten Antwort ist es oft einfacher, mit den kumulativen Wahrscheinlichkeiten . Die kumulativen Wahrscheinlichkeiten nehmen zu und sind für die Kombination benachbarter Kategorien unveränderlich. Außerdem ist , daher benötigen wir nur Modell- Wahrscheinlichkeiten. $J$ $i$ $Y_i$ $p_{ij}=P(Yi=j)$ $j=1,..., J$ $\gamma_{ij}=P(Y_i \le j)$ $\gamma_{iJ}=1$ $J–1$

Jetzt wollen wir s mit Kovariaten verknüpfen . In Ihrem Fall hat 3 bestellt Ebene: , , . Es ist sinnvoller, sie als bestellt und nicht als ungeordnet zu behandeln. Die restlichen Variablen sind Ihre Kovariaten. Das spezifische Modell, das Sie in Betracht ziehen, ist das Proportional-Odds-Modell und ist mathematisch äquivalent zu: $\gamma_{ij}$ $x$ Satlowmediumhigh

logit γ_{j} (x_{ich}) = θ_{j} - β^{T} x_{ich}, j = 1 \dots J - 1

$\mbox{logit } \gamma_j(x_i) = \theta_j - \beta^T x_i, j = 1 \ldots J-1$

woher γ_{j} (x_{ich}) = P ({Y.}_{ich} \leq j | x_{ich})

$\mbox{where }\gamma_j(x_i)=P(Y_i \le j | x_i)$

Es wird so genannt, weil die relativen Quoten für , das und sind: $Y \le j$ $x_1$ $x_2$

(\frac{γ_{j} (x_{1})}{1 - γ_{j} (x_{1})}) / (\frac{γ_{j} (x_{2})}{1 - γ_{j} (x_{2})}) = \exp (- β^{T} (x_{1} - x_{2}))

$\left(\frac {\gamma_j(x_1)}{1-\gamma_j(x_1)}\right) / \left(\frac {\gamma_j(x_2)}{1-\gamma_j(x_2)}\right)=\exp(-\beta^T (x_1-x_2))$

Beachten Sie, dass der obige Ausdruck nicht von abhängt . Natürlich muss die Annahme proportionaler Quoten für einen bestimmten Datensatz überprüft werden. $j$

Nun beantworte ich einige (1, 2, 4) Fragen.

Wie kann man verstehen, ob das Modell gut passt? Zusammenfassung (house.plr) zeigt Residual Deviance 3479.149 und AIC (Akaike Information Criterion?) von 3495.149. Ist das gut? Was ist ein gutes absolutes Maß für den Fall, dass diese nur als relative Maße (dh zum Vergleich mit einer anderen Modellanpassung) nützlich sind? Ist die Restabweichung ungefähr im Chi-Quadrat verteilt? Kann man "% korrekt vorhergesagt" für die Originaldaten oder eine Kreuzvalidierung verwenden? Was ist der einfachste Weg das zu tun?

Ein passendes Modell polrist ein besonderes glm, also alle Annahmen, die für ein traditionelles glmhier gelten. Wenn Sie sich richtig um die Parameter kümmern, können Sie die Verteilung herausfinden. Um zu testen, ob das Modell gut ist oder nicht , möchten Sie möglicherweise einen Anpassungstest durchführen , bei dem die folgende Null überprüft wird (beachten Sie, dass dies subtil ist. Meistens möchten Sie die Null ablehnen, aber dies ist hier nicht der Fall lehne es ab, um eine gute Anpassung zu erhalten):

H_{O} : aktuelles Modell ist gut genug

$H_o: \mbox{ current model is good enough }$

Sie würden dafür den Chi-Quadrat-Test verwenden . Der p-Wert wird erhalten als:

1-pchisq(deviance(house.plr),df.residual(house.plr))

In den meisten Fällen möchten Sie einen p-Wert von mehr als 0,05 erhalten, damit Sie die Null nicht ablehnen, um zu dem Schluss zu gelangen, dass das Modell gut passt (philosophische Korrektheit wird hier ignoriert).

Der AIC sollte für eine gute Anpassung hoch sein, während Sie keine große Anzahl von Parametern haben möchten. stepAICist ein guter Weg, um dies zu überprüfen.

Ja, Sie können definitiv Kreuzvalidierung verwenden, um festzustellen, ob die Vorhersagen zutreffen. Siehe predictFunktion (Option:) type = "probs"in ?polr. Sie müssen sich nur um die Kovariaten kümmern.

Welche Informationen enthält pr? Die Hilfeseite zum Profil ist allgemein gehalten und enthält keine Anleitung für polr

Enthält, wie von @chl und anderen angegeben, pralle Informationen, die zum Abrufen von CIs und anderer wahrscheinlichkeitsbezogener Informationen des polr fit. Alle glms werden unter Verwendung einer iterativ gewichteten Methode zur Schätzung der kleinsten Quadrate für die logarithmische Wahrscheinlichkeit angepasst. Bei dieser Optimierung erhalten Sie eine Vielzahl von Informationen (siehe Referenzen), die für die Berechnung der Varianz-Kovarianz-Matrix, des CI, des t-Werts usw. benötigt werden. Dies beinhaltet alles.

Wie interpretiert man die t-Werte für jeden Koeffizienten? Im Gegensatz zu einigen Modellanpassungen gibt es hier keine P-Werte.

Im Gegensatz zu einem normalen linearen Modell (speziell glm) haben andere glms nicht die schöne t-Verteilung für die Regressionskoeffizienten. Daher können Sie nur die Parameterschätzungen und ihre asymptotische Varianz-Kovarianz-Matrix unter Verwendung der Max-Likelihood-Theorie erhalten. Deshalb:

Varianz (\hat{β}) = (X^{T} W X)^{- 1} \hat{ϕ}

$\text{Variance}(\hat \beta) = (X^T W X)^{-1}\hat \phi$

Die Schätzung geteilt durch den Standardfehler ist das, was BDR und WV als t-Wert bezeichnen (ich gehe MASShier von einer Konvention aus). Es entspricht dem t-Wert aus der normalen linearen Regression, folgt jedoch keiner t-Verteilung. Mit CLT wird es asymptotisch normal verteilt. Aber sie bevorzugen es, diese ungefähre Zahl nicht zu verwenden (ich vermute), daher keine p-Werte. (Ich hoffe, dass ich mich nicht irre, und wenn, dann hoffe ich, dass der BDR nicht in diesem Forum ist. Ich hoffe außerdem, dass mich jemand korrigiert, wenn ich mich irre.)

suncoolsu
quelle

Ich werde mehr hinzufügen.

Suncoolsu

Danke dafür. Ich habe es mehrmals gelesen. Es bleiben noch viele Fragen offen. 1. Wie testet man funktional in R die Proportional-Odds-Annahme? 2. Sind Sie sicher, dass der Chi-Quadrat-Test richtig ist? In diesem Beispiel wird 0 zurückgegeben, was bedeutet .. crappy fit? Einige der t-Werte sind jedoch ziemlich hoch (InflHigh 10.1, InflMedium 5.4, ContHigh 3.7). 3. Was zeigen die Handlungen oder Paare?

Dfrankow

Danke für deine ausführliche Antwort suncoolsu. Ich bin in einer ähnlichen Situation und habe ein paar Fragen. 1. Ich erhalte auch 0 für jedes Modell, das Ihre Chi-Quadrat-Testgleichung verwendet. 2. Auf der Wikipedia-Seite zu AIC heißt es "Das bevorzugte Modell ist das mit dem minimalen AIC-Wert", aber Sie sagten, "AIC sollte für eine gute Anpassung hoch sein." Ich versuche, diese Konten in Einklang zu bringen.

Sam Swift

@dfrankow und @Sam Swift. Es tut mir leid, ich war ein bisschen beschäftigt mit dem Schreiben von Papieren. Ok - wenn Sie einen p-Wert = 0 erhalten, bedeutet dies, dass das Modell NICHT gut passt, da der Test der Anpassungsgüte fehlschlägt. In Bezug auf das Problem von AIC verwenden Wikipedia und ich eine andere Konvention. Ich benutze die, die von BDR und WV verwendet wird. (vgl. Erweiterung linearer Modelle in R, von Dr. Julian Faraway)

suncoolsu

Es gibt einige spezielle Fragen für 0/1 p - Werte und AIC Interpretation Sie hilfreich finden könnte: stats.stackexchange.com/questions/15223/... stats.stackexchange.com/questions/81427/...

Scott

Ich habe das Gespräch hier sehr genossen, aber ich habe das Gefühl, dass die Antworten nicht alle (sehr guten) Komponenten der von Ihnen gestellten Frage korrekt angesprochen haben. In der zweiten Hälfte der Beispielseite für polrdreht sich alles um die Profilerstellung. Eine gute technische Referenz sind hier Venerables und Ripley, die sich mit der Profilerstellung und ihren Funktionen befassen. Dies ist eine wichtige Technik, wenn Sie die Komfortzone der Anpassung exponentieller Familienmodelle mit voller Wahrscheinlichkeit (reguläre GLMs) verlassen.

Die Hauptabweichung hier ist die Verwendung von kategorialen Schwellenwerten. Sie werden feststellen, dass POLR keine übliche Intercept-Laufzeit einschätzt. Vielmehr gibt es Belästigungsparameter: Schwellenwerte, bei denen das angepasste Risiko tendenziell in eine bestimmte Summe der möglichen Kategorien fällt . Da diese Schwellenwerte niemals gemeinsam geschätzt werden, ist ihre Kovarianz mit Modellparametern unbekannt. Im Gegensatz zu GLMs können wir einen Koeffizienten nicht um einen bestimmten Betrag "stören" und sicher sein, wie er sich auf andere Schätzungen auswirken könnte. Wir verwenden die Profilerstellung, um die Belästigungsschwellen zu berücksichtigen. Profilierung ist ein immenses Thema, aber im Grunde ist das Ziel robust die Kovarianz von Regressionskoeffizienten zu messen , wenn das Modell einer irreguläre Wahrscheinlichkeit maximiert wird, wie mit , , , und $k-1$ $k$ lmernlspolrglm.nb.

Die Hilfeseite für ?profile.glmsollte von Nutzen sein, da es sich bei den polrObjekten im Wesentlichen um GLMs handelt (zuzüglich der kategorialen Schwellenwerte). Zuletzt können Sie den Quellcode, wenn er von Nutzen ist, mithilfe von erreichen getS3method('profile', 'polr'). Ich benutze diese getS3methodFunktion häufig, weil R darauf zu bestehen scheint, dass viele Methoden ausgeblendet werden sollten, aber man kann überraschend viel über Implementierung und Methoden lernen, indem man den Code überprüft.

• Welche Informationen enthält pr? Die Hilfeseite zum Profil ist allgemein gehalten und enthält keine Anleitung für polr.

prist ein profile.polr, profileObjekt (geerbte Klasse profile). Für jede Kovariate gibt es einen Eintrag. Der Profiler durchläuft jede Kovariate und berechnet die optimale Modellanpassung mit dieser Kovariate neu, die auf einen geringfügig anderen Wert festgelegt ist. Die Ausgabe zeigt den festen Wert der Kovariate, der als skalierte "Z-Score" -Differenz von ihrem geschätzten Wert und den resultierenden festen Effekten in anderen Kovariaten gemessen wird. Wenn Sie sich zum Beispiel anschauen pr$InflMedium, werden Sie feststellen, dass die anderen festen Effekte, wenn "z" 0 ist, dieselben sind wie in der ursprünglichen Anpassung.

• Was zeigt der Plot (pr)? Ich sehe sechs Grafiken. Jede hat eine numerische X-Achse, obwohl die Bezeichnung eine Indikatorvariable ist (sieht aus wie eine Eingabevariable, die ein Indikator für einen Ordnungswert ist). Dann ist die Y-Achse "Tau", was völlig unerklärt ist.

Wieder ?plot.profilegibt die Beschreibung. Das Diagramm zeigt ungefähr, wie die Regressionskoeffizienten kovary sind. tau ist die skalierte Differenz, davor der Z-Wert. Der Wert 0 gibt die optimalen Anpassungskoeffizienten an, die mit einem Häkchen gekennzeichnet sind. Sie würden nicht sagen, dass diese Passform so gut benommen ist, aber diese "Linien" sind tatsächlich Splines. Wenn sich die Wahrscheinlichkeit bei optimaler Anpassung sehr unregelmäßig verhält, würden Sie ein seltsames und unvorhersehbares Verhalten in der Handlung beobachten. Dies bedeutet, dass Sie die Ausgabe mit einer zuverlässigeren Fehlerschätzung (Bootstrap / Jackknife) abschätzen, CIs mit berechnen method='profile', Variablen neu codieren oder andere Diagnosen durchführen müssen.

• Was zeigen Paare (pr)? Es sieht aus wie ein Diagramm für jedes Paar von Eingabevariablen, aber ich sehe auch hier keine Erklärung für die X- oder Y-Achse.

In der Hilfedatei heißt es: "Die Methode" pairs "zeigt für jedes Parameterpaar x und y zwei Kurven, die sich bei der maximalen Wahrscheinlichkeitsschätzung schneiden und die Orte der Punkte angeben, an denen die Tangenten an die Konturen der Wahrscheinlichkeit des bivariaten Profils vertikal werden Im Fall einer genau bivariaten Normalprofilwahrscheinlichkeit wären diese beiden Kurven gerade Linien, die die bedingten Mittelwerte von y | x und x | y ergeben, und die Konturen wären genau elliptisch. Grundsätzlich helfen sie Ihnen wieder, die Vertrauensellipsen zu visualisieren. Nicht-orthogonale Achsen weisen auf stark kovariable Maße wie InfMedium und InfHigh hin, die intuitiv sehr verwandt sind. Auch hier würden unregelmäßige Wahrscheinlichkeiten zu Bildern führen, die hier ziemlich verwirrend sind.

• Wie kann man verstehen, ob das Modell gut passt? Zusammenfassung (house.plr) zeigt Residual Deviance 3479.149 und AIC (Akaike Information Criterion?) von 3495.149. Ist das gut? Was ist ein gutes absolutes Maß für den Fall, dass diese nur als relative Maße (dh zum Vergleich mit einer anderen Modellanpassung) nützlich sind? Ist die Restabweichung ungefähr im Chi-Quadrat verteilt? Kann man "% korrekt vorhergesagt" für die Originaldaten oder eine Kreuzvalidierung verwenden? Was ist der einfachste Weg das zu tun?

Eine gut zu bewertende Annahme ist die Proportional-Odds-Annahme. Dies spiegelt sich etwas im globalen Test wider (der polr gegen ein gesättigtes loglineares Modell bewertet). Eine Einschränkung hierbei ist, dass bei großen Datenmengen globale Tests immer fehlschlagen. Daher ist es eine gute Idee, Grafiken zu verwenden und Schätzungen (Betas) und Genauigkeiten (SEs) für das loglineare Modell und die Polr-Anpassung zu überprüfen. Wenn sie sich massiv widersprechen, stimmt vielleicht etwas nicht.

Bei geordneten Ergebnissen ist es schwierig, die prozentuale Übereinstimmung zu definieren. Wie wählen Sie einen Klassifikator basierend auf dem Modell aus und wie schätzen Sie die schlechte Leistung eines schlechten Klassifikators ein. modeist eine schlechte Wahl. Wenn ich 10 Kategorie-Logs habe und meine Vorhersage immer nur eine Kategorie abweicht, ist das vielleicht keine schlechte Sache. Außerdem kann mein Modell eine 40-prozentige Wahrscheinlichkeit einer 0-Antwort, aber auch eine 20-prozentige Wahrscheinlichkeit von 8, 9, 10 korrekt vorhersagen. Wenn ich also 9 beobachte, ist das gut oder schlecht? Wenn Sie die Übereinstimmung messen müssen, verwenden Sie einen gewichteten Kappa oder sogar MSE. Das loglineare Modell liefert immer die beste Übereinstimmung. Das macht die POLR nicht.

• Wie wendet man eine Anova auf dieses Modell an und interpretiert sie? In den Dokumenten heißt es: "Es gibt Methoden für die Standardfunktionen zur Modellanpassung, einschließlich" Vorhersagen "," Zusammenfassung "," VCOV "und" Anova "." Das Ausführen von anova (house.plr) führt jedoch dazu, dass anova nicht für ein einzelnes "polr" -Objekt implementiert wird

Sie können verschachtelte Modelle mit waldtestund lrtestim lmtestPaket in R testen . Dies entspricht ANOVA. Die Interpretation ist genauso wie bei GLMs.

• Wie interpretiert man die t-Werte für jeden Koeffizienten? Im Gegensatz zu einigen Modellanpassungen gibt es hier keine P-Werte.

Wiederum kann das POLR-Modell im Gegensatz zu linearen Modellen Probleme mit unregelmäßiger Wahrscheinlichkeit aufweisen, sodass eine auf dem Hessischen basierende Inferenz sehr instabil sein kann. Es ist analog zum Anpassen gemischter Modelle, siehe zum Beispiel die Hilfedatei confint.merModfür das lme4-Paket. Hier zeigen die mit der Profilerstellung vorgenommenen Bewertungen, dass sich die Kovarianz gut verhält. Die Programmierer hätten dies standardmäßig getan, mit der Ausnahme, dass die Profilerstellung sehr rechenintensiv sein kann und Sie es daher Ihren Händen überlassen. Wenn Sie die auf Wald basierende Schlussfolgerung sehen müssen, verwenden Sie sie coeftest(house.plr)aus dem lrtestPaket.

AdamO
quelle

Um die Proportionalitätsannahme in R zu 'testen' (dh auszuwerten), können Sie residuals.lrm () im Design-Paket von Frank Harrell Jr. verwenden. Wenn Sie? Residuals.lrm eingeben, gibt es ein Beispiel, wie Frank Harrell empfiehlt, die Proportionalitätsannahme zu bewerten (dh visuell anstatt durch einen Drucktastentest). Der Entwurf schätzt die geordneten logistischen Regressionen mithilfe von lrm (), das Sie durch polr () von MASS ersetzen können.

Ein formelleres Beispiel für das visuelle Testen der Proportionalitätsquotenannahme in R finden Sie in: Artikel: Regressionsmodelle für die ordinale Reaktion in der Ökologie Autor (en): Antoine Guisan und Frank E. Harrell Quelle: Journal of Vegetation Science, Vol. 3, No. 11, Nr. 5 (Okt. 2000), S. 617-626

mBrewster
quelle

Ich freue mich aufrichtig über Ihre Antwort. Der Zweck von StackExchange besteht jedoch darin, Antworten und keine Referenzen bereitzustellen. Statistiker scheinen für dieses Referenzproblem besonders anfällig zu sein. Gibt es Details, die Sie zur Verwendung von residuals.lrm hinzufügen können? Zum Beispiel ein Beispielbefehl und ein Beispiel für die Interpretation des Graphen für das house.plr-Beispiel?

Dfrankow

Update von der Website des Autors: "Das Design-Paket ist jetzt veraltet. R-Benutzer müssen stattdessen das rms-Paket verwenden". Mark, deine Antwort war sehr hilfreich für mich.

Tal Galili