Wie wird die Anpassung eines mit lme4 (> 1.0) ausgestatteten binomialen GLMM bewertet?

19

Ich besitze ein GLMM mit einer Binomialverteilung und einer Logit-Link-Funktion und habe das Gefühl, dass ein wichtiger Aspekt der Daten im Modell nicht gut dargestellt wird.

Um dies zu testen, möchte ich wissen, ob die Daten durch eine lineare Funktion auf der Logit-Skala gut beschrieben werden. Daher möchte ich wissen, ob die Residuen gut erzogen sind. Ich kann jedoch nicht herausfinden, bei welchen Residuen der Plot gezeichnet werden soll und wie der Plot zu interpretieren ist.

Beachten Sie, dass ich die neue Version von lme4 verwende ( die Entwicklungsversion von GitHub ):

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

Meine Frage ist: Wie überprüfe und interpretiere ich die Residuen eines binomial verallgemeinerten linearen gemischten Modells mit einer Logit-Link-Funktion?

Die folgenden Daten stellen nur 17% meiner realen Daten dar, aber das Anpassen dauert auf meinem Computer bereits rund 30 Sekunden, sodass ich es so belasse:

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)

Das einfachste plot ( ?plot.merMod) ergibt folgendes:

plot(m1)

Bildbeschreibung hier eingeben

Sagt mir das schon was?

Henrik
quelle
1
Ich könnte Zeit finden, um zurückzukommen und eine Pause einzulegen, aber ich denke, die allgemeine Antwort ist, dass es schwierig ist, mit den Residuen von Binärmodellen viel zu tun. Meine Haupt Entdeckung so weit davon entfernt , auf dem Grundstück auf ein bisschen Zoomen Sie oben haben, und das Hinzufügen einer geglätteten Linie (unter Verwendung von type=c("p","smooth")in plot.merMododer bewegen zu , ggplotwenn Sie Konfidenzintervall wollen) ist , dass es aussieht gibt es eine kleine , aber signifikante Muster, die Sie Ist möglicherweise in der Lage, das Problem zu beheben, indem Sie eine andere Verknüpfungsfunktion verwenden. Das war's schon ...
Ben Bolker
@ BenBolker Danke. Und können Sie dies und den Link zu freakonomics nicht einfach als Antwort auf die Frage posten? Dann würdest du wenigstens die 150 Punkte bekommen.
Henrik
3
Ich fand diesen CV-Thread, stats.stackexchange.com/questions/63566/… , sehr hilfreich. In diesem Beitrag wird erklärt, wie in R.
Nova
@Henrik Würden Sie mir bitte erklären, wie das Modell true ~ distance*(consequent+direction+dist)^2 + (direction+dist|V1)funktioniert? Wird das Modell give Schätzung der Wechselwirkung zwischen distance*consequent, distance*direction, distance*distund die Steigung von directionund dist daß Variiert mit V1? Was bedeutet das Quadrat in (consequent+direction+dist)^2?
ABC
@Henrik Ich habe deinen Code ausgeführt und er zeigt den 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). Warum ?
ABC

Antworten:

18

Kurze Antwort, da ich keine Zeit zum Besseren habe: Dies ist ein herausforderndes Problem. Binärdaten erfordern fast immer eine Art Binning oder Glättung, um die Anpassungsgüte zu beurteilen. Es war etwas hilfreich, fortify.lmerMod(von lme4, experimentell) in Verbindung mit ggplot2und insbesondere geom_smooth()zum Zeichnen im Wesentlichen der gleichen Residuen-gegen-Fit-Darstellung zu verwenden, die Sie oben haben, aber mit Konfidenzintervallen (ich habe auch die y-Grenzen etwas verkleinert, um auf die ( -5,5) Region). Das deutete auf eine systematische Variation hin, die durch Optimierung der Link-Funktion verbessert werden könnte. (Ich habe auch versucht, Residuen gegen die anderen Prädiktoren zu zeichnen, aber es war nicht allzu nützlich.)

Ich habe versucht, das Modell mit allen 3-Wege-Interaktionen zu versehen, aber es hat sich weder in Bezug auf die Abweichung noch in Bezug auf die Form der geglätteten Restkurve wesentlich verbessert.

Dann habe ich dieses bisschen rohe Gewalt benutzt, um Inverse-Link-Funktionen der Form für Bereich von 0,5 bis 2,0 zu testen: λ(logistic(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
}

Ich fand heraus, dass ein von 0,75 etwas besser war als das ursprüngliche Modell, wenn auch nicht signifikant - ich habe die Daten möglicherweise überinterpretiert.λ

Siehe auch: http://freakonometrics.hypotheses.org/8210

Ben Bolker
quelle
3

Dies ist ein weit verbreitetes Thema in Kursen zu Biostatistik / Epidemiologie, und es gibt aufgrund der Art des Modells keine sehr guten Lösungen dafür. Oft bestand die Lösung darin, eine detaillierte Diagnose unter Verwendung der Residuen zu vermeiden.

Ben schrieb bereits, dass für die Diagnose häufig Binning oder Smoothing erforderlich sind. Das Binning von Residuen ist (oder war) im R-Paket-Arm verfügbar, siehe z . B. diesen Thread . Darüber hinaus wurden einige Arbeiten durchgeführt, die vorhergesagte Wahrscheinlichkeiten verwenden. Eine Möglichkeit ist das Trennungsdiagramm, das zuvor in diesem Thread erläutert wurde . Diese könnten oder könnten nicht direkt in Ihrem Fall helfen, aber könnten möglicherweise die Interpretation helfen.

JTT
quelle
-1

Sie können AIC anstelle von Residuendiagrammen verwenden, um die Modellanpassung zu überprüfen. Mit dem Befehl in R: AIC (model1) erhalten Sie eine Nummer. Dann müssen Sie diese mit einem anderen Modell (z. B. mit mehr Prädiktoren) vergleichen - AIC (model2), das eine andere Nummer ergibt. Vergleichen Sie die beiden Ausgänge, und Sie möchten das Modell mit dem niedrigeren AIC-Wert.

Übrigens sind Dinge wie AIC und Log Likelihood Ratio bereits aufgelistet, wenn Sie die Zusammenfassung Ihres glmer-Modells erhalten, und beide geben Ihnen nützliche Informationen zur Modellanpassung. Sie möchten eine große negative Zahl für das Log-Likelihood-Verhältnis, um die Nullhypothese abzulehnen.

user108972
quelle
3
Dies wäre nützlicher, wenn OP versucht, konkurrierende Modelle zu vergleichen, dies jedoch nicht zu tun scheint und AIC nicht zur Bewertung der absoluten Modellanpassung verwendet werden kann.
Patrick Coulombe