Estou tendo problemas para entender a saída do meu lmer()
modelo. É um modelo simples de uma variável de resultado (Suporte) com interceptações de estado variáveis / efeitos aleatórios de estado:
mlm1 <- lmer(Support ~ (1 | State))
Os resultados de summary(mlm1)
são:
Linear mixed model fit by REML
Formula: Support ~ (1 | State)
AIC BIC logLik deviance REMLdev
12088 12107 -6041 12076 12082
Random effects:
Groups Name Variance Std.Dev.
State (Intercept) 0.0063695 0.079809
Residual 1.1114756 1.054265
Number of obs: 4097, groups: State, 48
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.13218 0.02159 6.123
Entendo que a variação do estado variável intercepta / efeitos aleatórios é 0.0063695
. Mas quando eu extraio o vetor desses efeitos aleatórios do estado e calculo a variação
var(ranef(mlm1)$State)
O resultado é 0.001800869
:, consideravelmente menor que a variação reportada por summary()
.
Tanto quanto eu o entendo, o modelo que eu especifiquei pode ser escrito:
Se isso estiver correto, a variação dos efeitos aleatórios ( ) deve ser σ 2 α . No entanto, estes não são realmente equivalentes no meu ajuste.lmer()
r
mixed-model
random-effects-model
lme4-nlme
nomad545
fonte
fonte
lmer()
? Parece que você postular que é estimado pela variação empírica dos efeitos aleatórios estimados alfa s . A descrição do seu modelo não está clara (perharps y i deve ser y i s ). É um design equilibrado?Respostas:
Esta é uma anova de sentido único clássica. Uma resposta muito curta para sua pergunta é que o componente de variação é composto de dois termos.
Portanto, o termo que você calculou é o primeiro termo no rhs (como efeitos aleatórios têm média zero). O segundo termo depende se REML de ML é usado e a soma dos erros padrão ao quadrado dos seus efeitos aleatórios.
fonte
1/48 * sum((se.ranef(mlm1)$State)^2)
- é0.004557198
. A variação das estimativas pontuais das ERs (obtidas, como acima, usandovar(ranef(mlm1)$State)
) é0.001800869
. A soma é0.006358067
, que é a variação relatada usandosummary()
nolmer()
modelo acima, com pelo menos 4 ou 5 dígitos. Muito obrigado @probabilityarm
pacote R para ase.ranef()
função.