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

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)

Czy to już mi coś mówi?

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.

$(\mbox{logistic}(x))^\lambda$$\lambda$

## 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
}

$\lambda$

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

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.

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ą.

Wykres dopasowany względem reszt nie powinien wykazywać żadnego (wyraźnego) wzoru. Wykres pokazuje, że model nie działa dobrze z danymi. Zobacz http://www.r-bloggers.com/model-validation-interpreting-residual-plots/