Jak mogę obliczyć Pearsona

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.$G^2$glm(..., family = binomial)$\chi^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?$\chi^2$

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:$\chi^2$glmlogistic.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 $G^2$statystyki testowe, tak samo jak deviance(logistic.fit)zapewnia.

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

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.

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