Jak ocenić dopasowanie dwumianowego GLMM wyposażonego w lme4 (> 1.0)?

19

Mam GLMM z rozkładem dwumianowym i funkcją linku logit i mam wrażenie, że ważny aspekt danych nie jest dobrze reprezentowany w modelu.

Aby to sprawdzić, chciałbym wiedzieć, czy dane są dobrze opisane przez funkcję liniową w skali logit. Dlatego chciałbym wiedzieć, czy reszty są dobrze wychowane. Nie mogę jednak dowiedzieć się, na których wykresach pozostały wykresy i jak interpretować wykresy.

Zauważ, że używam nowej wersji lme4 ( wersja rozwojowa od GitHub ):

packageVersion("lme4")
## [1] ‘1.1.0’

Moje pytanie brzmi: w jaki sposób mogę sprawdzić i zinterpretować pozostałości dwumianowych uogólnionych liniowych modeli mieszanych z funkcją logit link?

Następujące dane stanowią tylko 17% moich rzeczywistych danych, ale dopasowanie zajmuje już około 30 sekund na moim komputerze, więc zostawiam to w ten sposób:

require(lme4)
options(contrasts=c('contr.sum', 'contr.poly'))

dat <- read.table("http://pastebin.com/raw.php?i=vRy66Bif")
dat$V1 <- factor(dat$V1)

m1 <- glmer(true ~ distance*(consequent+direction+dist)^2 + (direction+dist|V1), dat, family = binomial)

Najprostszy wykres ( ?plot.merMod) daje następujące wyniki:

plot(m1)

wprowadź opis zdjęcia tutaj

Czy to już mi coś mówi?

binomial residuals glmm lme4-nlme Henrik
źródło
1
I może znaleźć czas, aby wrócić i podjąć pęknięcie na to, ale myślę, że ogólna odpowiedź jest taka, że jest to trudne do zrobienia, wiele z tych reszt z modeli binarnych. Moim głównym odkryciem tak daleko od powiększania nieco na działce masz powyżej i dodanie wygładzoną linię (stosując type=c("p","smooth")się plot.merModlub porusza się ggplot, jeśli chcesz przedziały ufności) jest to, że wygląda na to, że jest mały, ale znaczący wzór, który cię może być w stanie naprawić, przyjmując inną funkcję łącza. To wszystko na razie ...
Ben Bolker
@BenBolker Thanks. I czy możesz nie tylko opublikować to i link do freakonomics jako odpowiedź na pytanie? Wtedy przynajmniej zdobędziesz 150 punktów.
Henrik
3
Uważam, że ten wątek CV, stats.stackexchange.com/questions/63566/… , jest bardzo pomocny. W poście wyjaśniono, jak utworzyć skumulowany wykres resztek w R.
Nova
@Henrik Czy mógłbyś mi wyjaśnić, jak działa model true ~ distance*(consequent+direction+dist)^2 + (direction+dist|V1)? Czy oszacowanie dać model współdziałania distance*consequent, distance*direction, distance*disti nachylenie directiona dist , który zmienia się z V1? Co oznacza kwadrat (consequent+direction+dist)^2?
ABC
@Henrik Uruchomiłem twój kod i pokazuje on Warning message: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.123941 (tol = 0.001, component 1). Dlaczego ?
ABC

Odpowiedzi:

18

Krótka odpowiedź, ponieważ nie mam czasu na lepsze: jest to trudny problem; dane binarne prawie zawsze wymagają pewnego rodzaju binowania lub wygładzania w celu oceny dopasowania. Przydatne było użycie fortify.lmerMod(z lme4, eksperymentalnego) w połączeniu z, ggplot2a zwłaszcza geom_smooth()narysowanie zasadniczo tego samego wykresu resztkowo-dopasowanego, który masz powyżej, ale z przedziałami ufności (ja również nieco zawęziłem granice y, aby powiększyć ( -5,5) region). Sugerowało to pewne systematyczne zmiany, które można poprawić, modyfikując funkcję link. (Próbowałem też wykreślić wartości resztkowe w stosunku do innych predyktorów, ale nie było to zbyt przydatne).

Próbowałem dopasować model do wszystkich interakcji 3-kierunkowych, ale nie było to znacznej poprawy ani w odchyleniu, ani w kształcie wygładzonej krzywej resztkowej.

(logistyka(x))λλ

## uses (fragile) internal C calls for speed; could use plogis(),
##  qlogis() for readability and stability instead
logitpower <- function(lambda) {
    L <- list(linkfun=function(mu)
              .Call(stats:::C_logit_link,mu^(1/lambda),PACKAGE="stats"),
              linkinv=function(eta)
              .Call(stats:::C_logit_linkinv,eta,PACKAGE="stats")^lambda,
              mu.eta=function(eta) {
                  mu <-  .Call(stats:::C_logit_linkinv,eta,PACKAGE="stats")
                  mu.eta <-  .Call(stats:::C_logit_mu_eta,eta,PACKAGE="stats")
                  lambda*mu^(lambda-1)*mu.eta
              },
              valideta = function(eta) TRUE ,
              name=paste0("logit-power(",lambda,")"))
    class(L) <- "link-glm"
    L
}

λ

Zobacz także: http://freakonometrics.hypotheses.org/8210

Ben Bolker
źródło
3

Jest to bardzo powszechny temat na kursach biostatystycznych / epidemiologicznych i nie ma na to bardzo dobrych rozwiązań, zasadniczo ze względu na charakter modelu. Często rozwiązaniem było uniknięcie szczegółowej diagnostyki z wykorzystaniem pozostałości.

Ben już napisał, że diagnostyka często wymaga binowania lub wygładzania. Binning reszt jest dostępny (lub był) w ramieniu pakietu R, patrz np. Ten wątek . Ponadto wykonano pewne prace, które wykorzystują przewidywane prawdopodobieństwa; jedną z możliwości jest wykres separacji omówiony wcześniej w tym wątku . Mogą one pomóc lub nie bezpośrednio w twoim przypadku, ale mogą pomóc w interpretacji.

JTT
źródło
-1

Możesz użyć AIC zamiast resztkowych wykresów, aby sprawdzić dopasowanie modelu. Polecenie w R: AIC (model1) da ci liczbę ... więc musisz porównać to z innym modelem (na przykład z większą liczbą predyktorów) - AIC (model2), co da inną liczbę. Porównaj dwa wyjścia, a będziesz chciał modelu o niższej wartości AIC.

Nawiasem mówiąc, rzeczy takie jak AIC i współczynnik wiarygodności dziennika są już wymienione, gdy otrzymasz podsumowanie swojego modelu glitter, i oba dostarczą użytecznych informacji na temat dopasowania modelu. Chcesz, aby duża liczba ujemna dla współczynnika prawdopodobieństwa dziennika odrzuciła hipotezę zerową.

użytkownik108972
źródło
3
Byłoby to bardziej przydatne, gdyby OP próbował porównać konkurencyjne modele, ale nie wydaje się, że tak właśnie chcą, a AIC nie można użyć do oceny absolutnego dopasowania modelu.
Patrick Coulombe,