Como avaliar o ajuste de um GLMM binomial equipado com lme4 (> 1.0)?

19

Eu tenho um GLMM com uma distribuição binomial e uma função de link de logit e tenho a sensação de que um aspecto importante dos dados não está bem representado no modelo.

Para testar isso, eu gostaria de saber se os dados estão bem descritos por uma função linear na escala de logit. Por isso, gostaria de saber se os resíduos são bem comportados. No entanto, não consigo descobrir em que plotagem de resíduos traçar e como interpretar a trama.

Observe que estou usando a nova versão do lme4 ( a versão de desenvolvimento do GitHub ):

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

Minha pergunta é: como inspecionar e interpretar os resíduos de um modelo misto linear generalizado binomial com uma função de link logit?

Os dados a seguir representam apenas 17% dos meus dados reais, mas o ajuste já leva cerca de 30 segundos na minha máquina, então deixo assim:

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)

A plotagem mais simples ( ?plot.merMod) produz o seguinte:

plot(m1)

insira a descrição da imagem aqui

Isso já me diz alguma coisa?

Henrik
fonte
1
Talvez eu encontre tempo para voltar e resolver isso, mas acho que a resposta geral é que é difícil fazer muito com os resíduos de modelos binários. Meu principal descoberta tão longe de zoom em um pouco sobre o enredo que você tem acima, e adicionando uma linha alisada (usando type=c("p","smooth")em plot.merMod, ou movendo-se para ggplotse você quiser intervalos de confiança) é que parece que há um pequeno, mas significativo padrão, que você pode ser corrigido adotando uma função de link diferente. É isso até agora ...
Ben Bolker
@BenBolker Obrigado. E você não pode simplesmente postar isso e o link para a freakonomics como resposta à pergunta? Então, pelo menos, você obteria os 150 pontos.
Henrik
3
Eu achei este tópico do CV, stats.stackexchange.com/questions/63566/… , muito útil. O post explica como criar uma trama resíduos binned em R.
Nova
@ Henrik Você poderia me explicar como o modelo true ~ distance*(consequent+direction+dist)^2 + (direction+dist|V1)funciona? Será que a estimativa give modelo de interação entre distance*consequent, distance*direction, distance*diste inclinação de directione dist que varia com V1? O que o quadrado (consequent+direction+dist)^2denota?
ABC
@ Henrik Corri seu código e mostra o 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). Por quê ?
ABC

Respostas:

18

Resposta curta, já que não tenho tempo para melhorar: esse é um problema desafiador; dados binários quase sempre requerem algum tipo de classificação ou suavização para avaliar a qualidade do ajuste. Foi um pouco útil usar fortify.lmerMod(de lme4, experimental) em conjunto ggplot2e particularmente geom_smooth()desenhar essencialmente o mesmo gráfico residual versus ajustado que você tem acima, mas com intervalos de confiança (eu também reduzi um pouco os limites de y para ampliar o ( -5,5) região). Isso sugeriu alguma variação sistemática que poderia ser aprimorada ajustando a função de link. (Também tentei traçar resíduos contra os outros preditores, mas não foi muito útil.)

Tentei ajustar o modelo com todas as interações de três vias, mas não houve muita melhoria no desvio ou no formato da curva residual suavizada.

Então usei esse pouco de força bruta para tentar funções de link inverso da forma , para variando de 0,5 a 2,0: λ(logístico(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
}

Descobri que um de 0,75 era um pouco melhor que o modelo original, embora não significativamente - talvez eu tenha interpretado demais os dados.λ

Veja também: http://freakonometrics.hypotheses.org/8210

Ben Bolker
fonte
3

Esse é um tema muito comum nos cursos de bioestatística / epidemiologia, e não há soluções muito boas para isso, basicamente devido à natureza do modelo. Freqüentemente, a solução foi evitar diagnósticos detalhados usando os resíduos.

Ben já escreveu que os diagnósticos geralmente exigem binning ou suavização. O binning de resíduos está (ou estava) disponível no braço da embalagem R, consulte, por exemplo, esta rosca . Além disso, existem alguns trabalhos realizados que usam probabilidades previstas; uma possibilidade é o gráfico de separação que foi discutido anteriormente neste tópico . Eles podem ou não ajudar diretamente no seu caso, mas podem ajudar na interpretação.

JTT
fonte
-1

Você pode usar o AIC em vez de gráficos residuais para verificar o ajuste do modelo. Comando em R: AIC (modelo1), ele fornecerá um número ... então você precisará comparar isso com outro modelo (com mais preditores, por exemplo) - AIC (modelo2), que produzirá outro número. Compare as duas saídas e você desejará o modelo com o valor mais baixo da AIC.

A propósito, itens como AIC e razão de verossimilhança de log já estão listados quando você obtém o resumo do seu modelo glmer, e ambos fornecem informações úteis sobre o ajuste do modelo. Você deseja que um número negativo grande para a razão de verossimilhança de log rejeite a hipótese nula.

user108972
fonte
3
Isso seria mais útil se o OP estivesse tentando comparar modelos concorrentes, mas não parece que é isso que eles estão tentando fazer, e o AIC não pode ser usado para avaliar o ajuste absoluto do modelo.
Patrick Coulombe