Wykresy diagnostyczne dla regresji zliczania

Jakie wykresy diagnostyczne (i być może testy formalne) są najbardziej przydatne dla regresji, w których wynikiem jest zmienna licząca?

Szczególnie interesują mnie modele Poissona i modele dwumianowe ujemne, a także ich odpowiedniki zerowe i przeszkodowe. Większość źródeł, które znalazłem, po prostu kreśli wartości resztkowe w stosunku do dopasowanych wartości bez dyskusji o tym, jak te wykresy „powinny” wyglądać.

Mądrość i referencje bardzo mile widziane. Drugie pytanie dotyczy tego, dlaczego zadaję to pytanie, jeśli ma to znaczenie .

Powiązane dyskusje:

• Interpretowanie resztkowych wykresów diagnostycznych dla modeli GLM?
• Założenia uogólnionych modeli liniowych
• GLM - diagnostyka i która rodzina

Oto, co zazwyczaj lubię robić (dla ilustracji wykorzystuję rozproszone i niezbyt łatwo modelowane dane z czasów uczniów, których nie ma w szkole MASS):

1. Przetestuj i wykreśl oryginalne dane zliczania , wykreślając obserwowane częstotliwości i dopasowane częstotliwości (patrz rozdział 2 w Friendly ), który jest obsługiwany przez vcdpakiet Rw dużych częściach. Na przykład za pomocą goodfiti a rootogram:

   library(MASS)
   library(vcd)
   data(quine) 
   fit <- goodfit(quine$Days) 
   summary(fit) 
   rootogram(fit)

   lub z wykresami porządkowymi, które pomagają określić, który model danych zliczania leży u podstaw (np. tutaj nachylenie jest dodatnie, a punkt przecięcia jest dodatni, co przemawia za ujemnym rozkładem dwumianowym):

   Ord_plot(quine$Days)

   lub z wykresami „XXXXXXness”, gdzie XXXXX jest rozkładem wyboru, powiedzmy wykres Poissona (który przemawia przeciwko Poissonowi, spróbuj także type="nbinom"):

   distplot(quine$Days, type="poisson")

2. Sprawdź zwykłe miary dopasowania (takie jak statystyki współczynnika wiarygodności w porównaniu z modelem zerowym lub podobnym):

   mod1 <- glm(Days~Age+Sex, data=quine, family="poisson")
   summary(mod1)
   anova(mod1, test="Chisq")

3. Sprawdź nadmierną / niską dyspersję , patrząc na residual deviance/dfformalną statystykę testową lub na nią (np. Zobacz tę odpowiedź ). Tutaj mamy wyraźnie nadmierną dyspersję:

   library(AER)
   deviance(mod1)/mod1$df.residual
   dispersiontest(mod1)

4. Sprawdź, wpływowych i dźwigni punktów , np, przy czym influencePlotw caropakowaniu. Oczywiście tutaj wiele punktów ma duży wpływ, ponieważ Poisson jest złym modelem:

   library(car)
   influencePlot(mod1)

5. Sprawdź zerową inflację , dopasowując model danych zliczania i jego odpowiednik z zerową inflacją / przeszkodą i porównaj je (zwykle z AIC). Tutaj model z napompowaniem zerowym pasowałby lepiej niż zwykły Poisson (znowu prawdopodobnie z powodu nadmiernej dyspersji):

   library(pscl)
   mod2 <- zeroinfl(Days~Age+Sex, data=quine, dist="poisson")
   AIC(mod1, mod2)

6. Narysuj wartości resztkowe (surowe, odchylenie lub skalowane) na osi y w stosunku do (log) przewidywanych wartości (lub predyktora liniowego) na osi x. Widzimy tutaj bardzo duże reszty i znaczne odchylenie reszty dewiacji od normy (przemawiając przeciwko Poissonowi; edycja: @ FlorianHartig odpowiedź sugeruje, że nie należy oczekiwać normalności tych reszty, więc nie jest to rozstrzygająca wskazówka):

   res <- residuals(mod1, type="deviance")
   plot(log(predict(mod1)), res)
   abline(h=0, lty=2)
   qqnorm(res)
   qqline(res)

7. W razie zainteresowania wykreśl połowę normalnego wykresu prawdopodobieństwa reszt, wykreślając uporządkowane absolutne reszty w stosunku do oczekiwanych wartości normalnych Atkinsona (1981) . Specjalną cechą byłaby symulacja „linii” odniesienia i obwiedni z symulowanymi / ładowanymi przedziałami ufności (choć nie pokazano):

   library(faraway)
   halfnorm(residuals(mod1))

8. $\pm$

   plot(Days~Age, data=quine) 
   prs  <- predict(mod1, type="response", se.fit=TRUE)
   pris <- data.frame("pest"=prs[[1]], "lwr"=prs[[1]]-prs[[2]], "upr"=prs[[1]]+prs[[2]])
   points(pris$pest ~ quine$Age, col="red")
   points(pris$lwr  ~ quine$Age, col="pink", pch=19)
   points(pris$upr  ~ quine$Age, col="pink", pch=19)

To powinno dać ci wiele przydatnych informacji o twojej analizie i większość kroków działa dla wszystkich standardowych rozkładów danych zliczania (np. Poissona, Dwumianu ujemnego, COM Poissona, Prawa mocy).

Podejście polegające na stosowaniu standardowych wykresów diagnostycznych, ale chcących wiedzieć, jak powinny wyglądać, podoba mi się:

 Buja, A., Cook, D. Hofmann, H., Lawrence, M. Lee, E.-K., Swayne,
 D.F and Wickham, H. (2009) Statistical Inference for exploratory
 data analysis and model diagnostics Phil. Trans. R. Soc. A 2009
 367, 4361-4383 doi: 10.1098/rsta.2009.0120

Jednym z wymienionych tam podejść jest utworzenie kilku symulowanych zestawów danych, w których założenia są prawdziwe, i utworzenie wykresów diagnostycznych dla tych symulowanych zestawów danych, a także utworzenie wykresu diagnostycznego dla rzeczywistych danych. umieść wszystkie te wykresy na ekranie w tym samym czasie (losowo umieszczając wykres na podstawie rzeczywistych danych). Teraz masz wizualne odniesienie do tego, jak powinny wyglądać wykresy, a jeśli założenia dotyczą prawdziwych danych, to wykres ten powinien wyglądać tak samo jak inne (jeśli nie możesz powiedzieć, które są prawdziwe dane, wówczas testowane założenia są prawdopodobnie bliskie wystarczy prawda), ale jeśli wykres danych rzeczywistych wyraźnie różni się od drugiego, oznacza to, że przynajmniej jedno z założeń się nie sprawdza. vis.testFunkcji w pakiecie dla R TeachingDemos umożliwia wdrożenie tego jako testu.

To stare pytanie, ale pomyślałem, że przydałoby się dodać, że mój pakiet DHARMa R (dostępny z CRAN, patrz tutaj ) zapewnia teraz znormalizowane resztki dla GLM i GLMM, oparte na podejściu symulacyjnym podobnym do sugerowanego przez @GregSnow .

Z opisu pakietu:

Pakiet DHARMa wykorzystuje podejście oparte na symulacji do tworzenia łatwo interpretowalnych skalowanych reszt z dopasowanych uogólnionych liniowych modeli mieszanych. Obecnie obsługiwane są wszystkie klasy „merMod” z „lme4” („lmerMod”, „glmerMod”), „glm” (w tym „negbin” z „MASS”, ale z wyłączeniem quasi-dystrybucji) i klasy modeli „lm”. Alternatywnie można również wykonać symulacje stworzone zewnętrznie, np. Symulacje predykcyjne późniejsze z oprogramowania bayesowskiego, takie jak „JAGS”, „STAN” lub „BŁĘDY”. Otrzymane reszty są standaryzowane na wartości od 0 do 1 i mogą być interpretowane tak intuicyjnie, jak reszty z regresji liniowej. Pakiet zawiera również szereg funkcji wykresu i testowania typowego problemu z błędną specyfikacją modelu,

@Momo - możesz zaktualizować zalecenie 6, jest to mylące. Normalności odchyleń odchyleń zasadniczo nie oczekuje się w przypadku Poissona , jak wyjaśniono w winiecie DHARMa lub tutaj ; i dlatego wyszukiwanie odchyleń odchyleń (lub jakichkolwiek innych standardowych resztek), które różnią się od linii prostej na wykresie qqnorm, nie jest zatem generalnie żadnym problemem . Pakiet DHARMa zapewnia wykres qq, który jest niezawodny do diagnozowania odchyleń od Poissona lub innych rodzin GLM. Stworzyłem przykład ilustrujący problem z reszt dewiacji tutaj .

W glm.diag.plotspakiecie wywoływana jest funkcja bootdo generowania wykresów diagnostycznych dla GLM. Co to robi:

Tworzy wykres odchyleń odchylenia scyzoryka od predyktora liniowego, wykresy wyników normalnych odchyleń znormalizowanych, wykres przybliżonych statystyk Cooka względem dźwigni / (1-dźwigni) i wykres przypadku statystyki Cooka.