Entendo que usamos modelos de efeitos aleatórios (ou efeitos mistos) quando acreditamos que alguns parâmetros do modelo variam aleatoriamente em algum fator de agrupamento. Desejo ajustar um modelo em que a resposta tenha sido normalizada e centralizada (não perfeitamente, mas bastante próxima) em um fator de agrupamento, mas uma variável independente x
não tenha sido ajustada de forma alguma. Isso me levou ao teste a seguir (usando dados fabricados ) para garantir que encontraria o efeito que estava procurando se realmente estivesse lá. Eu executei um modelo de efeitos mistos com uma interceptação aleatória (entre os grupos definidos por f
) e um segundo modelo de efeito fixo com o fator f como um preditor de efeito fixo. Eu usei o pacote R lmer
para o modelo de efeito misto e a função baselm()
para o modelo de efeito fixo. A seguir estão os dados e os resultados.
Observe que y
, independentemente do grupo, varia em torno de 0. E isso x
varia de acordo com o y
grupo, mas varia muito mais entre grupos do quey
> data
y x f
1 -0.5 2 1
2 0.0 3 1
3 0.5 4 1
4 -0.6 -4 2
5 0.0 -3 2
6 0.6 -2 2
7 -0.2 13 3
8 0.1 14 3
9 0.4 15 3
10 -0.5 -15 4
11 -0.1 -14 4
12 0.4 -13 4
Se você estiver interessado em trabalhar com os dados, aqui está a dput()
saída:
data<-structure(list(y = c(-0.5, 0, 0.5, -0.6, 0, 0.6, -0.2, 0.1, 0.4,
-0.5, -0.1, 0.4), x = c(2, 3, 4, -4, -3, -2, 13, 14, 15, -15,
-14, -13), f = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L,
4L, 4L, 4L), .Label = c("1", "2", "3", "4"), class = "factor")),
.Names = c("y","x","f"), row.names = c(NA, -12L), class = "data.frame")
Ajustando o modelo de efeitos mistos:
> summary(lmer(y~ x + (1|f),data=data))
Linear mixed model fit by REML
Formula: y ~ x + (1 | f)
Data: data
AIC BIC logLik deviance REMLdev
28.59 30.53 -10.3 11 20.59
Random effects:
Groups Name Variance Std.Dev.
f (Intercept) 0.00000 0.00000
Residual 0.17567 0.41913
Number of obs: 12, groups: f, 4
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.008333 0.120992 0.069
x 0.008643 0.011912 0.726
Correlation of Fixed Effects:
(Intr)
x 0.000
Observo que o componente de variação de interceptação é estimado em 0 e, o que x
é mais importante para mim, não é um preditor significativo de y
.
Em seguida, ajustei o modelo de efeito fixo f
como um preditor em vez de um fator de agrupamento para uma interceptação aleatória:
> summary(lm(y~ x + f,data=data))
Call:
lm(formula = y ~ x + f, data = data)
Residuals:
Min 1Q Median 3Q Max
-0.16250 -0.03438 0.00000 0.03125 0.16250
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.38750 0.14099 -9.841 2.38e-05 ***
x 0.46250 0.04128 11.205 1.01e-05 ***
f2 2.77500 0.26538 10.457 1.59e-05 ***
f3 -4.98750 0.46396 -10.750 1.33e-05 ***
f4 7.79583 0.70817 11.008 1.13e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1168 on 7 degrees of freedom
Multiple R-squared: 0.9484, Adjusted R-squared: 0.9189
F-statistic: 32.16 on 4 and 7 DF, p-value: 0.0001348
Agora percebo que, como esperado, x
é um preditor significativo de y
.
O que estou procurando é intuição sobre essa diferença. De que maneira meu pensamento está errado aqui? Por que espero incorretamente encontrar um parâmetro significativo para x
esses dois modelos, mas apenas o vejo no modelo de efeito fixo?
x
variável não seja significativa. Suspeito que seja o mesmo resultado (coeficientes e SE) que você teria conseguido executarlm(y~x,data=data)
. Não há mais tempo para diagnosticar, mas gostaria de salientar isso.Respostas:
Há várias coisas acontecendo aqui. Essas são questões interessantes, mas será necessário bastante tempo / espaço para explicar tudo.
Primeiro de tudo, isso tudo se torna muito mais fácil de entender se traçarmos os dados . Aqui está um gráfico de dispersão em que os pontos de dados são coloridos por grupo. Além disso, temos uma linha de regressão específica para cada grupo, bem como uma linha de regressão simples (ignorando grupos) em negrito:
O modelo de efeito fixo
lm()
O modelo misto
Aqui estão os coeficientes para o modelo de regressão simples (a linha em negrito tracejada no gráfico):
Como você pode ver, os coeficientes aqui são idênticos aos que obtivemos no modelo misto. Isso é exatamente o que esperávamos encontrar, pois, como você já observou, temos uma estimativa de 0 para as interceptações aleatórias, fazendo a relação mencionada anteriormente / correlação intra-classe 0. Portanto, as estimativas do modelo misto neste caso são apenas as estimativas de regressão linear simples e, como podemos ver no gráfico, a inclinação aqui é muito menos pronunciada do que as inclinações dentro do cluster.
Isso nos leva a uma questão conceitual final ...
Por que a variação das interceptações aleatórias é estimada em 0?
A resposta a esta pergunta pode se tornar um pouco técnica e difícil, mas tentarei mantê-la o mais simples e não técnica possível (pelo bem de nós!). Mas talvez ainda seja um pouco demorado.
O modelo misto que estamos considerando não está usando o método de correlação intra-classe para representar a dependência nos dados. Em vez disso, descreve a dependência em termos de componentes de variação . Tudo bem, desde que a correlação intra-classe seja positiva. Nesses casos, a correlação intra-classe pode ser facilmente escrita em termos de componentes de variância, especificamente como a razão mencionada anteriormente da variância de interceptação aleatória para a variância total. (Veja a página wiki sobre correlação intra-classepara obter mais informações sobre isso.) Mas, infelizmente, os modelos de componentes de variância têm dificuldade em lidar com situações em que temos uma correlação intra-classe negativa. Afinal, escrever a correlação intra-classe em termos dos componentes da variação envolve escrevê-la como uma proporção da variação, e as proporções não podem ser negativas.
Então o que nós podemos fazer?
Para fazer isso, levamos nossax xb x xw x
Finalmente, ainda temos uma estimativa de 0 para a variação das interceptações aleatórias, pelas razões que elaborei na seção anterior. Não tenho muita certeza do que podemos fazer com isso, pelo menos, sem mudar para outro software que não seja
lmer()
, e também não sei até que ponto isso ainda afetará adversamente nossas estimativas neste modelo misto final. Talvez outro usuário possa entrar em contato com algumas idéias sobre esse problema.Referências
fonte
lme
restringe por padrão seja> = 0? Veja esta pergunta e sua resposta selecionada , ou seja, ajustando uma correlação de simetria composta através de umgls
ajuste ou configuraçãocorrelation = corCompSymm(form = ~1|f)
emlme
Após considerável contemplação, acredito ter descoberto minha própria resposta. Acredito que um economista definiria minha variável independente como endógena e, assim, estaria correlacionada com as variáveis independentes e as dependentes. Nesse caso, essas variáveis são omitidas ou não observadas . No entanto, observo os agrupamentos entre os quais a variável omitida deve variar.
Acredito que o economista sugeriria um modelo de efeito fixo . Ou seja, um modelo que inclua um manequim para cada nível de agrupamento (ou uma especificação equivalente que condicione o modelo de modo que muitos manequins de agrupamento não sejam necessários) neste caso. Com um modelo de efeito fixo, a esperança é que todas as variáveis não observadas e invariantes no tempo possam ser controladas condicionando as variações de grupo (ou de indivíduo). De fato, o segundo modelo em minha pergunta é precisamente um modelo de efeito fixo e, como tal, fornece a estimativa que espero.
Congratulo-me com comentários que iluminarão ainda mais essa circunstância.
fonte