Jak rozumieć dane wyjściowe z funkcji polr R (uporządkowana regresja logistyczna)?

Jestem nowy w R, uporządkowałem regresję logistyczną i polr.

Sekcja „Przykłady” u dołu strony pomocy dla polr (która pasuje do modelu regresji logistycznej lub probitowej do uporządkowanej odpowiedzi czynnikowej) pokazuje

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)

Jakie informacje przawierają? Strona pomocy w profilu jest ogólna i nie zawiera wskazówek dla polr.
Co plot(pr)pokazuje Widzę sześć wykresów. Każda z nich ma oś X, która jest numeryczna, chociaż etykieta jest zmienną wskaźnikową (wygląda jak zmienna wejściowa, która jest wskaźnikiem wartości porządkowej). Wtedy oś Y jest „tau”, co jest całkowicie niewyjaśnione.
Co pairs(pr)pokazuje Wygląda jak wykres dla każdej pary zmiennych wejściowych, ale znowu nie widzę wyjaśnienia osi X lub Y.
Jak zrozumieć, czy model dobrze pasował? summary(house.plr)pokazuje Residual Deviance 3479.149 i AIC (Akaike Information Criterion?) z 3495.149. Czy to dobrze? W przypadku, gdy są one użyteczne jedynie jako miary względne (tj. Do porównania z dopasowaniem innego modelu), jaka jest dobra miara bezwzględna? Czy odchylenie resztkowe jest w przybliżeniu rozkładem chi-kwadrat? Czy można użyć „% poprawnie przewidywanych” na oryginalnych danych lub jakiejś weryfikacji krzyżowej? Jak najłatwiej to zrobić?
Jak stosować i interpretować anovaten model? Dokumenty mówią: „Istnieją metody dla standardowych funkcji dopasowania modelu, w tym przewidywania, podsumowania, vcov, anova”. Jednak uruchomienie anova(house.plr)powodujeanova is not implemented for a single "polr" object
Jak interpretować wartości t dla każdego współczynnika? W przeciwieństwie do niektórych modeli, nie ma tutaj wartości P.

Zdaję sobie sprawę, że jest to wiele pytań, ale sensownie jest zadawać je jako jeden pakiet („jak korzystać z tej rzeczy?”), A nie 7 różnych pytań. Wszelkie informacje są mile widziane.

r logistic

— dfrankow
źródło

@dfrankow Nieco prymitywna i na pewno bardzo częściowa pomoc na pierwsze dwa pytania, ale methods("profile")da ci metody (w tym przypadku S3) związane z profileobiektem R , wtedy zobaczysz, że istnieje specjalna metoda na polrwyniki, którą możesz przeglądać online, pisząc getAnywhere("profile.polr")w wierszu polecenia R.

— chl

Dzięki! Kod źródłowy jest dobry. Wyjaśnienie byłoby jeszcze lepsze. :)

— dfrankow

Ktoś wskazał mi „Modern Applied Statistics with S” Venables i Ripley. Sekcja 7.3 zawiera „Przykład czterokierunkowej tabeli częstotliwości”, który w szerokim zakresie obejmuje ten model domu. Czytanie ..

— dfrankow

W rzeczywistości sekcja jest „proporcjonalnym modelem kursów”

— dfrankow

Odpowiedzi:

Proponuję przyjrzeć się książkom na temat analizy danych kategorycznych (por. Analiza danych kategorycznych Alana Agrestiego, 2002), aby lepiej wyjaśnić i zrozumieć uporządkowaną regresję logistyczną . Na wszystkie zadane pytania odpowiada w zasadzie kilka rozdziałów takich książek. Jeśli jesteś zainteresowany tylko w Rpowiązanych przykładach Rozszerzanie modele liniowe w R przez Juliana recz (CRC Press, 2008) jest doskonałym odniesienia.

Zanim odpowiem na pytania, uporządkowana regresja logistyczna jest przypadkiem wielomianowych modeli logitowych, w których kategorie są uporządkowane. Załóżmy, że mamy uporządkowane kategorie i indywidualnego , z porządkowej reakcji , dla . Z uporządkowaną odpowiedzią często łatwiej jest pracować ze skumulowanymi prawdopodobieństwami, $J$ $i$ $Y_i$ $p_{ij}=P(Yi=j)$ $j=1,..., J$ . Skumulowane prawdopodobieństwa są rosnące i niezmienne dla łączenia sąsiednich kategorii. Ponadto, , więc potrzebujemy tylkoprawdopodobieństwmodelu . $\gamma_{ij}=P(Y_i \le j)$ $\gamma_{iJ}=1$ $J–1$

Teraz chcemy połączyć z współzmiennymi . W twoim przypadku, ma 3 poziomy zamówione: , , . Bardziej sensowne jest traktowanie ich jako uporządkowanych niż nieuporządkowanych. Pozostałe zmienne są twoimi współzmiennymi. Konkretny model, który rozważasz, jest modelem proporcjonalnych szans i jest matematycznie równoważny z: $\gamma_{ij}$ $x$ Satlowmediumhigh

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

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

where γ_{j} (x_{i}) = P (Y_{i} \leq j | x_{i})

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

Jest to tak zwane, ponieważ szanse względne dla porównaniu i wynoszą: $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))$

$j$

Teraz odpowiem na niektóre (1, 2, 4) pytania.

Jak zrozumieć, czy model dobrze pasował? Podsumowanie (house.plr) pokazuje Residual Deviance 3479.149 i AIC (Akaike Information Criterion?) z 3495.149. Czy to dobrze? W przypadku, gdy są one użyteczne jedynie jako miary względne (tj. Do porównania z dopasowaniem innego modelu), jaka jest dobra miara bezwzględna? Czy odchylenie resztkowe jest w przybliżeniu rozkładem chi-kwadrat? Czy można użyć „% poprawnie przewidywanych” na oryginalnych danych lub jakiejś weryfikacji krzyżowej? Jak najłatwiej to zrobić?

Dopasowany model polrjest wyjątkowy glm, więc wszystkie założenia, które dotyczą tradycyjnego glmtrzymania się tutaj. Jeśli odpowiednio zadbasz o parametry, możesz ustalić rozkład. W szczególności, aby sprawdzić, czy model jest dobry, czy nie, możesz wykonać test dobroci dopasowania , który przetestuje następujący zerowy (zauważ, że jest to subtelny, przeważnie chcesz odrzucić zerowy, ale tutaj nie chcesz odrzuć, aby uzyskać dobre dopasowanie):

H_{o} : current model is good enough

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

Można by użyć testu chi-kwadrat dla tego produktu. Wartość p otrzymuje się jako:

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

Przez większość czasu chciałbyś uzyskać wartość p większą niż 0,05, aby nie odrzucać wartości zerowej, aby stwierdzić, że model jest dobrze dopasowany (poprawność filozoficzna jest tutaj ignorowana).

AIC powinno być wysokie, aby dobrze pasować, a jednocześnie nie chcesz mieć dużej liczby parametrów. stepAICto dobry sposób na sprawdzenie tego.

Tak, zdecydowanie możesz użyć weryfikacji krzyżowej, aby sprawdzić, czy prognozy się utrzymują. Patrz predictfunkcja (opcja:) type = "probs"w ?polr. Trzeba tylko zadbać o zmienne towarzyszące.

Jakie informacje zawiera pr? Strona pomocy w profilu jest ogólna i nie zawiera wskazówek dla polr

Jak wskazał @chl i inni, przawiera wszystkie informacje potrzebne do uzyskania CI oraz inne informacje związane z prawdopodobieństwem polr fit. Wszystkie glms są dopasowane przy użyciu iteracyjnie ważonej metody oszacowania najmniejszych kwadratów dla prawdopodobieństwa dziennika. W tej optymalizacji uzyskujesz wiele informacji (zobacz referencje), które będą potrzebne do obliczenia macierzy kowariancji wariancji, CI, wartości t itp. Obejmuje to wszystko.

Jak interpretować wartości t dla każdego współczynnika? W przeciwieństwie do niektórych modeli> pasuje, nie ma tutaj wartości P.

W przeciwieństwie do normalnego modelu liniowego (specjalnego glm) inne glmnie mają ładnego rozkładu t dla współczynników regresji. Dlatego wszystko, co możesz uzyskać, to oszacowania parametrów i ich asymptotyczna macierz wariancji kowariancji przy użyciu teorii maksymalnego prawdopodobieństwa. W związku z tym:

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

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

Oszacowanie podzielone przez błąd standardowy to, co BDR i WV nazywają wartością t ( MASStutaj zakładam konwencję). Jest to równoważne wartości t z normalnej regresji liniowej, ale nie podlega rozkładowi t. Używając CLT, jest on asymptotycznie normalnie rozłożony. Ale wolą nie używać tego przybliżonego (tak sądzę), stąd brak wartości p. (Mam nadzieję, że się nie mylę, a jeśli tak, mam nadzieję, że BDR nie jest na tym forum. Mam nadzieję, że ktoś mnie poprawi, jeśli się mylę).

— suncoolsu
źródło

Dodam więcej.

— suncoolsu,

Dzięki za to. Przeczytałem to kilka razy. Pozostaje mnóstwo pytań. 1. Jak funkcjonalnie w R testuje się założenie o proporcjonalnych szansach? 2. Czy jesteś pewien, że test chi-kwadrat jest prawidłowy? W tym przykładzie zwraca 0, co oznacza… pasujące dopasowanie? Ale niektóre z wartości t są dość wysokie (InflHigh 10.1, InflMedium 5.4, ContHigh 3.7). 3. Co pokazują wykresy lub pary?

— dfrankow

Dzięki za wyczerpującą odpowiedź suncoolsu. Jestem w podobnej sytuacji i mam kilka pytań. 1. Dostaję również 0 dla każdego modelu za pomocą równania testu chi-sq. 2. Strona wikipedii na AIC mówi „preferowany model to ten o minimalnej wartości AIC”, ale powiedziałeś: „AIC powinno być wysokie, aby dobrze pasowało”. Próbuję pogodzić te konta.

— Sam Swift,

@dfrankow i @Sam Swift. Przepraszam, byłem trochę zajęty pisaniem dokumentów. Ok - jeśli otrzymasz wartość p = 0, oznacza to, że model NIE jest dobrym dopasowaniem, ponieważ test poprawności dopasowania nie powiedzie się. Jeśli chodzi o problem AIC, wikipedia i ja używamy do tego innej konwencji. Używam tego, który jest używany przez BDR i WV. (por. Rozszerzanie modeli liniowych w R., dr. Julian Faraway)

— suncoolsu,

Istnieje kilka dedykowanych pytań dotyczących wartości 0/1 p i interpretacji AIC, które mogą Ci się przydać: stats.stackexchange.com/questions/15223/... stats.stackexchange.com/questions/81427/...

— Scott

Bardzo podobała mi się rozmowa tutaj, jednak uważam, że odpowiedzi nie odpowiedziały poprawnie na wszystkie (bardzo dobre) elementy postawionego pytania. Druga połowa przykładowej strony polrdotyczy profilowania. Dobre referencje techniczne tutaj to Czcigodni i Ripley, którzy omawiają profilowanie i jego działanie. Jest to technika krytyczna, gdy wychodzisz poza strefę komfortu dopasowania wykładniczych modeli rodzinnych z pełnym prawdopodobieństwem (zwykłe GLM).

$k-1$ $k$ lmernls, polri glm.nb.

Strona pomocy dla ?profile.glmpowinna być użyteczna, ponieważ polrobiekty są zasadniczo GLM (plus progi kategorialne). Wreszcie, możesz faktycznie osiągnąć szczyt w kodzie źródłowym, jeśli ma to jakiekolwiek zastosowanie, używając getS3method('profile', 'polr'). Często używam tej getS3methodfunkcji, ponieważ choć R wydaje się nalegać, aby wiele metod było ukrytych, można zaskakująco wiele nauczyć się o implementacji i metodach, przeglądając kod.

• Jakie informacje zawiera pr? Strona pomocy w profilu jest ogólna i nie zawiera wskazówek dla polr.

prjest profile.polr, profileobiektem (dziedziczona klasa profile). Dla każdego współzmiennego znajduje się wpis. Profiler zapętla się nad każdą zmienną towarzyszącą, ponownie oblicza optymalne dopasowanie modelu z tą zmienną zmienną ustawioną na nieco inną wartość. Dane wyjściowe pokazują stałą wartość współzmiennej mierzoną jako przeskalowana różnica „z-score” od jej wartości szacunkowej i wynikające z tego ustalone efekty w innych współzmiennych. Na przykład, jeśli spojrzysz na to pr$InflMedium, zauważysz, że gdy „z” wynosi 0, inne ustalone efekty są takie same jak w oryginalnym dopasowaniu.

• Co pokazuje fabuła (pr)? Widzę sześć wykresów. Każda z nich ma oś X, która jest numeryczna, chociaż etykieta jest zmienną wskaźnikową (wygląda jak zmienna wejściowa, która jest wskaźnikiem wartości porządkowej). Wtedy oś Y jest „tau”, co jest całkowicie niewyjaśnione.

Ponownie ?plot.profilepodaje opis. Wykres w przybliżeniu pokazuje, w jaki sposób współczynniki regresji są kowalencyjne. tau to wyskalowana różnica, wcześniejszy wynik Z, więc jego wartość 0 daje optymalne współczynniki dopasowania, przedstawione za pomocą znaku podziałki. Nie powiedziałbyś, że to dopasowanie jest tak dobrze wychowane, ale te „linie” są tak naprawdę splajnami. Gdyby prawdopodobieństwo było bardzo nieregularnie zachowane przy optymalnym dopasowaniu, zaobserwowałbyś dziwne i nieprzewidywalne zachowanie na wykresie. To sprawi, że oszacujesz wynik za pomocą bardziej niezawodnego oszacowania błędu (bootstrap / jackknife), do obliczenia CI za pomocą method='profile', do przekodowania zmiennych lub do przeprowadzenia innej diagnostyki.

• Co pokazują pary (pr)? Wygląda jak wykres dla każdej pary zmiennych wejściowych, ale znowu nie widzę wyjaśnienia osi X lub Y.

Plik pomocy mówi: „Metoda par pokazuje, dla każdej pary parametrów x i y, dwie krzywe przecinające się przy oszacowaniu maksymalnego prawdopodobieństwa, które podają loci punktów, w których prawdopodobieństwo stycznych do konturów profilu dwuwymiarowego staje się pionowe i odpowiednio poziomo. W przypadku dokładnie dwuwymiarowego prawdopodobieństwa normalnego profilu te dwie krzywe byłyby liniami prostymi dającymi warunkowe środki y | x i x | y, a kontury byłyby dokładnie eliptyczne. " Zasadniczo ponownie pomagają w wizualizacji elips pewności. Osie nieortogonalne wskazują bardzo zmienne miary, takie jak InfMedium i InfHigh, które są intuicyjnie bardzo powiązane. Ponownie, nieregularne prawdopodobieństwa prowadziłyby do dość zaskakujących obrazów.

• Jak zrozumieć, czy model dobrze pasował? Podsumowanie (house.plr) pokazuje Residual Deviance 3479.149 i AIC (Akaike Information Criterion?) z 3495.149. Czy to dobrze? W przypadku, gdy są one użyteczne jedynie jako miary względne (tj. Do porównania z dopasowaniem innego modelu), jaka jest dobra miara bezwzględna? Czy odchylenie resztkowe jest w przybliżeniu rozkładem chi-kwadrat? Czy można użyć „% poprawnie przewidywanych” na oryginalnych danych lub jakiejś weryfikacji krzyżowej? Jak najłatwiej to zrobić?

Jednym z założeń, które warto ocenić, jest założenie o proporcjonalnych szansach. Znajduje to odzwierciedlenie w teście globalnym (który ocenia polr względem nasyconego modelu loglinearnego). Ograniczeniem jest tutaj to, że przy dużych danych globalne testy zawsze kończą się niepowodzeniem. W związku z tym dobrym pomysłem jest używanie grafiki i sprawdzanie szacunków (beta) i precyzji (SE) dla modelu logicznego i dopasowania polr. Jeśli masowo się nie zgadzają, być może coś jest nie tak.

Przy uporządkowanych wynikach trudno jest określić procentową zgodność. W jaki sposób wybierzesz klasyfikator na podstawie modelu, a jeśli tak, to w jaki sposób będziesz wykazywał słabą wydajność od słabego klasyfikatora. modeto zły wybór. Jeśli mam 10 logów kategorii, a moje prognozy są zawsze tylko jedna kategoria wyłączona, być może nie jest to zła rzecz. Co więcej, mój model może poprawnie przewidzieć 40% szansy na odpowiedź 0, ale także 20% szans na 8, 9, 10. Więc jeśli zauważę 9, czy to dobrze, czy źle? Jeśli musisz zmierzyć zgodność, użyj ważonej kappa lub nawet MSE. Model logiczny zawsze zapewni najlepszą zgodność. Nie to robi POLR.

• Jak stosować i interpretować anova w tym modelu? Dokumenty mówią: „Istnieją metody dla standardowych funkcji dopasowania modelu, w tym przewidywania, podsumowania, vcov, anova”. Jednak uruchomienie anova (house.plr) powoduje, że anova nie jest zaimplementowana dla pojedynczego obiektu „polr”

Można przetestować modele z zagnieżdżonych waldtestoraz lrtestw lmtestpakiet w R. Jest to odpowiednik ANOVA. Interpretacja jest dokładnie taka sama jak w przypadku GLM.

• Jak interpretować wartości t dla każdego współczynnika? W przeciwieństwie do niektórych modeli, nie ma tutaj wartości P.

Ponownie, w przeciwieństwie do modeli liniowych, model POLR może mieć problemy z nieregularnym prawdopodobieństwem, więc wnioskowanie oparte na Hesji może być bardzo niestabilne. Jest to analogiczne do dopasowania modeli mieszanych, patrz na przykład plik pomocy confint.merModdla pakietu lme4. Oceny dokonane za pomocą profilowania pokazują, że kowariancja jest dobrze zachowana. Programiści zrobiliby to domyślnie, z wyjątkiem tego, że profilowanie może być bardzo intensywne obliczeniowo, a zatem pozostawiają go w Twoich rękach. Jeśli musisz zobaczyć wnioskowanie oparte na Wald, skorzystaj coeftest(house.plr)z lrtestpakietu.

— AdamO
źródło

Aby „przetestować” (tj. Ocenić) założenie o proporcjonalnych szansach w R, możesz użyć residuals.lrm () w pakiecie projektowym Franka Harrella Jr. Jeśli wpiszesz? Residuals.lrm, istnieje szybki do odtworzenia przykład, w jaki sposób Frank Harrell zaleca ocenę założenia proporcjonalności szans (tj. Wizualnie, a nie za pomocą testu przycisku). Szacunki projektowe uporządkowane regresje logistyczne za pomocą lrm (), które można zastąpić polr () z MASS.

Bardziej formalny przykład tego, jak wizualnie przetestować założenie o proporcjonalnych szansach w R, patrz: Paper: Ordinal Response Regression Models in Ecology Autor (autorzy): Antoine Guisan i Frank E. Harrell Źródło: Journal of Vegetation Science, t. 11, nr 5 (październik 2000), str. 617-626

— mBrewster
źródło

Szczerze doceniam twoją odpowiedź. Jednak celem StackExchange jest udzielenie odpowiedzi, a nie referencji. Statystycy wydają się szczególnie podatni na ten problem referencyjny. Czy możesz dodać jakieś szczegóły dotyczące korzystania z pliku residuals.lrm? Na przykład przykładowe polecenie i przykład interpretacji wykresu dla przykładu house.plr?

— dfrankow

Aktualizacja ze strony autora: „Pakiet Design jest już nieaktualny. Użytkownicy R muszą zamiast tego użyć pakietu rms”. Mark, twoja odpowiedź była dla mnie bardzo pomocna.

— Tal Galili