Comment évaluer l'ajustement d'un GLMM binomial équipé de lme4 (> 1.0)?

19

J'ai un GLMM avec une distribution binomiale et une fonction de lien logit et j'ai le sentiment qu'un aspect important des données n'est pas bien représenté dans le modèle.

Pour tester cela, je voudrais savoir si les données sont bien décrites par une fonction linéaire sur l'échelle logit. Par conséquent, je voudrais savoir si les résidus se comportent bien. Cependant, je ne peux pas savoir à quel résidu tracer à tracer et comment interpréter le tracé.

Notez que j'utilise la nouvelle version de lme4 ( la version de développement de GitHub ):

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

Ma question est la suivante: comment inspecter et interpréter les résidus d'un modèle mixte linéaire généralisé binomial avec une fonction de lien logit?

Les données suivantes ne représentent que 17% de mes données réelles, mais le montage prend déjà environ 30 secondes sur ma machine, je le laisse donc comme ceci:

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)

Le tracé le plus simple ( ?plot.merMod) produit les éléments suivants:

plot(m1)

entrez la description de l'image ici

Est-ce que cela me dit déjà quelque chose?

Henrik
la source
1
Je vais peut-être trouver le temps de revenir et de faire un essai, mais je pense que la réponse générale est qu'il est difficile de faire beaucoup avec les résidus des modèles binaires. Ma principale découverte si loin de zoomer sur un peu sur le terrain que vous avez ci - dessus, et en ajoutant une ligne lissée ( en utilisant type=c("p","smooth")dans plot.merMod, ou de passer à ggplotsi vous voulez des intervalles de confiance) est qu'il semble qu'il y ait un modèle petit mais significatif, que vous pourrait être en mesure de corriger en adoptant une fonction de lien différente. C'est tout pour l'instant ...
Ben Bolker
@BenBolker Merci. Et ne pouvez-vous pas simplement poster ceci et le lien vers la freakonomics en réponse à la question? Ensuite, au moins, vous obtiendriez les 150 points.
Henrik
3
J'ai trouvé ce fil de CV, stats.stackexchange.com/questions/63566/… , très utile. Le billet explique comment créer un tracé des résidus en bac dans R.
Nova
@Henrik Pourriez-vous s'il vous plaît m'expliquer comment fonctionne le modèle true ~ distance*(consequent+direction+dist)^2 + (direction+dist|V1)? Est-ce que l'estimation de donner modèle d'interaction entre distance*consequent, distance*direction, distance*distet la pente directionet dist qui varie avec V1? Que signifie le carré (consequent+direction+dist)^2?
ABC du
@Henrik J'ai exécuté votre code et il montre le 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). Pourquoi ?
ABC du

Réponses:

18

Réponse courte car je n'ai pas le temps de faire mieux: c'est un problème difficile; les données binaires nécessitent presque toujours une sorte de regroupement ou de lissage pour évaluer la qualité de l'ajustement. Il a été quelque peu utile d'utiliser fortify.lmerMod(à partir de lme4, expérimental) conjointement avec ggplot2et en particulier geom_smooth()pour dessiner essentiellement le même graphique résiduel vs ajusté que vous avez ci-dessus, mais avec des intervalles de confiance (j'ai également resserré un peu les limites y pour zoomer sur le ( -5,5) région). Cela suggère une variation systématique qui pourrait être améliorée en modifiant la fonction de liaison. (J'ai également essayé de tracer des résidus par rapport aux autres prédicteurs, mais ce n'était pas trop utile.)

J'ai essayé d'ajuster le modèle avec toutes les interactions à 3 voies, mais ce n'était pas vraiment une amélioration de la déviance ou de la forme de la courbe résiduelle lissée.

Ensuite, j'ai utilisé ce peu de force brute pour essayer les fonctions de liaison inverse de la forme , pour allant de 0,5 à 2,0: λ(la logistique(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
}

J'ai trouvé qu'un de 0,75 était légèrement meilleur que le modèle d'origine, mais pas de manière significative - j'ai peut-être surinterprété les données.λ

Voir aussi: http://freakonometrics.hypotheses.org/8210

Ben Bolker
la source
3

C'est un thème très courant dans les cours de biostatistique / épidémiologie, et il n'y a pas de très bonnes solutions pour cela, essentiellement en raison de la nature du modèle. Souvent, la solution a été d'éviter les diagnostics détaillés en utilisant les résidus.

Ben a déjà écrit que les diagnostics nécessitent souvent un binning ou un lissage. Le regroupement des résidus est (ou était) disponible dans le bras du package R, voir par exemple ce fil . De plus, certains travaux sont effectués et utilisent des probabilités prédites; une possibilité est le tracé de séparation qui a été discuté plus tôt dans ce fil . Ceux-ci pourraient ou pourraient ne pas aider directement dans votre cas, mais pourraient aider à l'interprétation.

JTT
la source
-1

Vous pouvez utiliser l'AIC au lieu des tracés résiduels pour vérifier l'ajustement du modèle. Commandez dans R: AIC (modèle1), il vous donnera un nombre ... alors vous devez le comparer avec un autre modèle (avec plus de prédicteurs, par exemple) - AIC (modèle2), qui donnera un autre nombre. Comparez les deux sorties, et vous voudrez le modèle avec la valeur AIC inférieure.

Soit dit en passant, des éléments comme l'AIC et le rapport de vraisemblance du journal sont déjà répertoriés lorsque vous obtenez le résumé de votre modèle glmer, et les deux vous donneront des informations utiles sur l'ajustement du modèle. Vous voulez qu'un grand nombre négatif pour le rapport de vraisemblance logarithmique rejette l'hypothèse nulle.

user108972
la source
3
Cela serait plus utile si OP essayait de comparer des modèles concurrents, mais il ne semble pas que ce soit ce qu'ils essaient de faire, et AIC ne peut pas être utilisé pour évaluer l'ajustement absolu du modèle.
Patrick Coulombe
-3

Le tracé ajusté par rapport aux résidus ne doit présenter aucun motif (clair). Le graphique montre que le modèle ne fonctionne pas bien avec les données. Voir http://www.r-bloggers.com/model-validation-interpreting-residual-plots/

lankoski
la source
1
Je pense que vous vous trompez, un certain modèle est prédit par la fonction de lien: stats.stackexchange.com/q/25068/442
Henrik
1
Peut-être utile: freakonometrics.hypotheses.org/8210
Ben Bolker