Jak mogę obliczyć Pearsona

10

Wskaźnik prawdopodobieństwa (inaczej dewiacja) Statystyka i test braku dopasowania (lub dobroci dopasowania) jest dość prosty do uzyskania dla modelu regresji logistycznej (dopasowanie przy użyciu funkcji) w R. Jednak może być łatwe jest, aby niektóre liczby komórek były wystarczająco niskie, aby test był niewiarygodny. Jednym ze sposobów weryfikacji wiarygodności testu współczynnika wiarygodności dla braku dopasowania jest porównanie jego statystyki testowej i wartości P z wynikami testu chi-Pearsona (lub ) braku dopasowania.G2glm(..., family = binomial)χ2

Ani glmobiekt, ani jego summary()metoda nie podają statystyki testu testu chi kwadrat Pearsona z powodu braku dopasowania. W moich poszukiwaniach jedyną rzeczą, jaką wymyśliłem, jest chisq.test()funkcja (w statspakiecie): w jej dokumentacji jest napisane: „ chisq.testwykonuje testy tabeli kontyngencji chi-kwadrat i testy dopasowania”. Jednak dokumentacja jest ograniczona, jak przeprowadzić takie testy:

Jeśli xjest to macierz z jednym rzędem lub kolumną, lub jeśli xjest wektorem i ynie jest podana, wówczas przeprowadzany jest test zgodności pod względem dopasowania ( xjest traktowany jako jednowymiarowa tabela kontyngencji). Wpisy xmuszą być liczbami całkowitymi nieujemnymi. W tym przypadku testowana hipoteza sprawdza, czy prawdopodobieństwo populacji jest równe tym p, czy są równe, jeśli pnie zostanie podane.

Wyobrażam sobie, że możesz użyć ykomponentu glmobiektu jako xargumentu chisq.test. Nie można jednak użyć fitted.valueskomponentu glmobiektu jako pargumentu chisq.test, ponieważ pojawi się błąd: „ probabilities must sum to 1.

Jak mogę (w R) przynajmniej obliczyć statystykę testową Pearson kątem braku dopasowania bez konieczności ręcznego wykonywania kroków?χ2

Ogień Ogień
źródło

Odpowiedzi:

13

Suma kwadratów reszt Pearsona jest dokładnie równa statystyce testu Pearsona z powodu braku dopasowania. Jeśli więc zostanie wywołany dopasowany model (tj. Obiekt) , poniższy kod zwróci statystyki testowe:χ2glmlogistic.fit

sum(residuals(logistic.fit, type = "pearson")^2)

Zobacz dokumentację, residuals.glmaby uzyskać więcej informacji, w tym jakie inne pozostałości są dostępne. Na przykład kod

sum(residuals(logistic.fit, type = "deviance")^2)

dostaniesz G2statystyki testowe, tak samo jak deviance(logistic.fit)zapewnia.

Ogień Ogień
źródło
Zgadzam się z Macro. Jeśli chcesz uzyskać odpowiedzi od grupy, powinieneś poczekać, aby usłyszeć, co inni mają do powiedzenia w pierwszej kolejności. Wszystko, co możesz teraz uzyskać, jest stronnicze, gdy zobaczysz swoją odpowiedź. Ponadto, jeśli znasz odpowiedź, co próbujesz udowodnić, robiąc to?
Michael R. Chernick
@Macro - Firefeather opublikował cztery pytania na tej stronie (w tym na tym) i udzielił odpowiedzi na trzy z nich (w tym na tę), akceptując jedną ze swoich odpowiedzi raz. Jeszcze kilka takich i mogę zacząć widzieć wzór!
jbowman
@jbowman, mogę sobie wyobrazić odpowiedź na twoje własne pytanie, jeśli wymyślisz je samodzielnie, zanim ktokolwiek inny opublikuje odpowiedź, ale firefeather opublikował odpowiedź dosłownie mniej niż 5 minut po opublikowaniu pytania, wygląda na to, że nie potrzebował pomocy , co sprawiło, że zapytałem, dlaczego ... Naprawdę nie rozumiem motywacji ...
Makro
3
@Macro, proszę zobaczyć ten oficjalny link: blog.stackoverflow.com/2011/07/... ( link znajduje się na stronie Zadaj pytanie w etykiecie pola wyboru u dołu: „Odpowiedz na własne pytanie - podziel się swoją wiedzą, pytaniami i odpowiedziami w stylu „). Miałem to pytanie podczas odrabiania lekcji (wybrałem użycie R zamiast Minitabu, chociaż Minitab został zademonstrowany na zajęciach), ale nie miałem wystarczająco dużo czasu, aby napisać pytanie i czekać na odpowiedź. Zrozumiałem to obejście i postanowiłem podzielić się nim ze społecznością.
Firefeather
2
@Macro, nie ma za co! Chciałbym móc zadawać więcej pytań, na które nie udzielam odpowiedzi, i odpowiadać na pytania, których nie zadałem. Ale jbowman „prawda o wzór: moje składki do społeczności skłaniając się ku mówię do siebie. :) (Przynajmniej w jakiś sposób przyczyniam się do społeczności, prawda?)
Firefeather
10

Statystyka Pearsona ma rozkład zdegenerowany, więc nie jest ogólnie zalecana dla dopasowania modelu logistycznego. Wolę testy strukturalne (liniowość, addytywność). Jeśli chcesz testu zbiorczego, zobacz pojedynczy stopień swobody le Cessie - van Houwelingen - Copas - Hosmer test nieważonej sumy kwadratów zaimplementowany w funkcji rmspakietu R.residuals.lrm

Frank Harrell
źródło
2
-1: Dziękuję za wgląd! To jednak nie odpowiada na moje pytanie. Ponieważ jest to odpowiedni komentarz / dyskusja na temat wypowiedzi, którą wypowiedziałem w tle mojego pytania, twoja odpowiedź prawdopodobnie należy do komentarza, a nie do odpowiedzi.
Firefeather
2
Myślę, że cztery osoby, które głosowały na moją odpowiedź, nie zgadzają się z tobą. I nie poradziłeś sobie z rozkładem zdegenerowanym.
Frank Harrell,
@FrankHarrell: Czy ta miara GOF różni się od testu GOF Hosmer-Lemeshow (HL)? Zakładając, że to z powodu nazwy, a także porównałem dwa: Przeprowadzony test HL GOF, jak znaleziono w ResourceSelectionpakiecie, a jego wynik różni się od tego, co otrzymuję resid(lrm_object, 'gof')po uruchomieniu mojego modelu regresji logistycznej jako lrm_object <- lrm(...). Jeśli rzeczywiście są różne, czy możesz skomentować, w jaki sposób test HL ma się do tego, o którym tu wspominasz? Dziękuję Ci!
Meg
1
Te dwa są bardzo różne. Statystyka HL (obecnie przestarzała) ma ustalony współczynnik df i zwykle opiera się na decylach ryzyka. HLχ2 Statystyka nie jest więc zdegenerowana jak N. Z drugiej strony, uważajcie na jakiekolwiekχ2 Statystyka, w której df stale się rozwija N.
Frank Harrell,
Bardzo chciałbym zobaczyć symulację pokazującą tę degenerację.
wdkrnls
0

Dzięki, nie zdawałem sobie sprawy, że to było tak proste jak: suma (reszty (f1, type = "pearson") ^ 2) Jednak zauważ, że reszta Pearsona różni się w zależności od tego, czy jest obliczana według grupy współzmiennej, czy indywidualnie. Prosty przykład:

m1 jest matrycą (ta jest głową większej matrycy):

m1 [1: 4,1: 8]

    x1 x2 x3 obs    pi   lev  yhat y
obs  1  1 44   5 0.359 0.131 1.795 2
obs  0  1 43  27 0.176 0.053 4.752 4
obs  0  1 53  15 0.219 0.062 3.285 1
obs  0  1 33  22 0.140 0.069 3.080 3

Gdzie x1-3 są predyktorami, obs jest nie. obserwacje w każdej grupie, pi jest prawdopodobieństwem przynależności do grupy (przewidywanym na podstawie równania regresji), lew jest dźwignią, przekątną macierzy kapelusza, czyli przewidywaną nie. (y = 1) w grupie iy faktyczna liczba nie.

To da ci Pearsona według grup. Zauważ, że jest inaczej, jeśli y == 0: ' 'fun1 <- function(j){        if (m1[j,"y"] ==0){ # y=0 for this covariate pattern     Pr1 <- sqrt( m1[i,"pi"] / (1-m1[i,"pi"]))     Pr2 <- -sqrt (m1[i,"obs"]) res <- round( Pr1 * Pr2, 3) return(res) } else {  Pr1 <- m1[j,"y"] - m1[j,"yhat"] Pr2 <- sqrt(   m1[j,"yhat"] * ( 1-(m1[j,"pi"]) )   ) res <- round( Pr1/Pr2, 3) return(res) }    }

A zatem

nr <-nrow(m1)
nr1 <- seq(1,nr)
Pr <- sapply(1:nrow[m1], FUN=fun1)
PrSj <- sum(Pr^2)

Jeśli istnieje duża liczba osobników o wzorcach współzmiennych y = 0, wówczas resztkowa wartość Pearons będzie znacznie większa, jeśli zostanie obliczona przy użyciu metody „według grupy”, a nie „według osoby”.

Patrz np. Hosmer & Lemeshow „Applied Logistic Regression”, Wiley, 200.

dardisco
źródło
0

Możesz także użyć c_hat(mod)tego, co da ten sam wynik, co sum(residuals(mod, type = "pearson")^2).

użytkownik54098
źródło
1
W której paczce znajduje c_hatsię?
Firefeather