Comparando modelos de efeitos mistos e efeitos fixos (testando a significância dos efeitos aleatórios)

10

Dadas três variáveis, ye x, que são positivas contínuas e z, que são categóricas, tenho dois modelos candidatos dados por:

fit.me <- lmer( y ~ 1 + x + ( 1 + x | factor(z) ) )

e

fit.fe <- lm( y ~ 1 + x )

Espero comparar esses modelos para determinar qual modelo é mais apropriado. Parece-me que, em certo sentido, fit.feestá aninhado por dentro fit.me. Normalmente, quando esse cenário geral ocorre, um teste qui-quadrado pode ser realizado. Em R, podemos executar este teste com o seguinte comando,

anova(fit.fe,fit.me)

Quando ambos os modelos contêm de efeitos aleatórios (gerados por lmerdo lme4pacote), o anova()comando funciona bem. Devido aos parâmetros de contorno, normalmente é aconselhável testar a estatística resultante do qui-quadrado via simulação, no entanto, ainda podemos usar a estatística no procedimento de simulação.

Quando os dois modelos contêm apenas efeitos fixos, essa abordagem - e o anova()comando associado - funcionam bem.

No entanto, quando um modelo contém efeitos aleatórios e o modelo reduzido contém apenas efeitos fixos, como no cenário acima, o anova()comando não funciona.

Mais especificamente, eu recebo o seguinte erro:

 > anova(fit.fe, fit.me)
 Error: $ operator not defined for this S4 class

Há algo de errado em usar a abordagem Chi-Square de cima (com simulação)? Ou isso é simplesmente um problema de anova()não saber lidar com modelos lineares gerados por diferentes funções?

Em outras palavras, seria apropriado gerar manualmente a estatística do qui-quadrado derivada dos modelos? Em caso afirmativo, quais são os graus de liberdade adequados para comparar esses modelos? Pela minha conta:

F=((SSEredvocêced-SSEfvocêeueu)/(p-k))((SSEfvocêeueu)/(n-p-1 1))Fp-k,n-p-1 1

Estamos estimando dois parâmetros no modelo de efeitos fixos (inclinação e interceptação) e mais dois parâmetros (parâmetros de variação para a inclinação aleatória e interceptação aleatória) no modelo de efeitos mistos. Tipicamente, o parâmetro de intercepção não é contado nos graus de liberdade computação, de modo que implica que e ; tendo dito que não tenho certeza se os parâmetros de variação para os parâmetros de efeitos aleatórios devem ser incluídos nos cálculos dos graus de liberdade; as estimativas de variância para parâmetros de efeito fixo não são consideradas , mas acredito que seja porque as estimativas de parâmetro para efeitos fixos são consideradas constantes desconhecidas enquanto são consideradas variáveis ​​aleatórias desconhecidask=1 1p=k+2=3para efeitos mistos. Eu gostaria de receber alguma ajuda sobre esse assunto.

Finalmente, alguém tem uma Rsolução ( baseada em) mais apropriada para comparar esses modelos?

user9171
fonte
4
Se você substituir lm()por gls()do nlmepacote e lmer()por lme()(novamente do nlmepacote), tudo funcionará bem. Mas observe que você obterá um teste conservador (valores- p muito grandes ), pois os parâmetros para o modelo mais simples estão no limite do espaço de parâmetros. E realmente a escolha de incluir os efeitos aleatórios deve basear-se na teoria (por exemplo, no plano de amostragem), não em um teste estatístico.
Karl Ove Hufthammer 14/03
11
O que você quer fazer com os modelos? Um modelo pode ser melhor para alguns propósitos e o outro modelo para outros fins. Todos os modelos estão errados; portanto, a questão não é qual modelo está certo, mas qual é mais útil para o seu problema específico.
Kodiologist
11
@ Kodiologist Basicamente, quero garantir que as estimativas de parâmetros para os efeitos fixos sejam confiáveis. Os erros padrão podem não ser confiáveis ​​se as observações forem consideradas independentes. Além disso, seria bom fazer uma declaração sobre como a variável é o efeito aleatório, mas acho que isso não é tão essencial.
user9171
2
@ user9171 Uma boa maneira de verificar a estabilidade (confiabilidade) nas estimativas de parâmetros de um modelo é usar a inicialização. Faça o gráfico das distribuições de inicialização para cada parâmetro que os dois modelos compartilham, com um gráfico por parâmetro e modelo. Distribuições mais apertadas implicam maior estabilidade. Você provavelmente descobrirá que o modelo mais simples gera estimativas mais estáveis, porque menos parâmetros permitem uma estimativa mais precisa de cada parâmetro.
Kodiologist

Respostas:

6

Tecnicamente, você pode fazê-lo funcionar apenas alternando a ordem dos parâmetros:

> anova(fit.me, fit.fe) 

Vai funcionar muito bem. Se você passar um objeto gerado pela lmerprimeira vez, anova.merModele será chamado em vez de anova.lm(que não sabe como lidar com lmerobjetos). Vejo:

?anova.merMod

Embora a escolha de um modelo misto ou fixo seja uma opção de modelagem que precise levar em consideração o projeto experimental, não um problema de seleção de modelo. Consulte https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#testing-significance-of-random-effects do BenBolker para obter mais detalhes:

Considere não testar a importância dos efeitos aleatórios.

witek
fonte
+1. Tomei a liberdade de inserir um link para o FAQ do @BenBolker, que contém algumas discussões e referências adicionais.
ameba