Adicionar efeito aleatório influencia as estimativas do coeficiente

Sempre fui ensinado que efeitos aleatórios influenciam apenas a variação (erro) e que efeitos fixos influenciam apenas a média. Mas encontrei um exemplo em que efeitos aleatórios também influenciam a média - a estimativa do coeficiente:

require(nlme)
set.seed(128)
n <- 100
k <- 5
cat <- as.factor(rep(1:k, each = n))
cat_i <- 1:k # intercept per kategorie
x <- rep(1:n, k)
sigma <- 0.2
alpha <- 0.001
y <- cat_i[cat] + alpha * x + rnorm(n*k, 0, sigma)
plot(x, y)

# simulate missing data
y[c(1:(n/2), (n*k-n/2):(n*k))] <- NA

m1 <- lm(y ~ x)
summary(m1)

m2 <- lm(y ~ cat + x)
summary(m2)

m3 <- lme(y ~ x, random = ~ 1|cat, na.action = na.omit)
summary(m3)

Você pode ver que o coeficiente estimado para o xmodelo m1é -0,013780, enquanto o modelo m3é 0,0011713 - ambos significativamente diferentes de zero.

Observe que, quando removo a linha que simula os dados ausentes, os resultados são os mesmos (é matriz completa).

Por que é que?

PS: observe que eu não sou um estatístico profissional; portanto, se você está prestes a responder com muita matemática, faça também um resumo simples para manequins :-)

"Eu sempre fui ensinado que efeitos aleatórios influenciam apenas a variação (erro), e que efeitos fixos influenciam apenas a média".

Como você descobriu, isso só é verdade para conjuntos de dados completos e equilibrados (ou seja, sem dados ausentes) sem preditores contínuos. Em outras palavras, para os tipos de dados / modelos discutidos nos textos clássicos da ANOVA. Sob essas circunstâncias ideais, os efeitos fixos e aleatórios podem ser estimados independentemente um do outro.

Quando essas condições não se mantêm (como muitas vezes não acontecem no "mundo real"), os efeitos fixos e aleatórios não são independentes. Como um aparte interessante, é por isso que os modelos mistos "modernos" são estimados usando métodos de otimização iterativos, em vez de serem exatamente resolvidos com um pouco de álgebra matricial, como no caso clássico da ANOVA mista: para estimar os efeitos fixos, precisamos conheça os efeitos aleatórios, mas para estimar os efeitos aleatórios, precisamos conhecer os efeitos fixos! Mais relevante para a presente questão, isso também significa que, quando os dados estão desequilibrados / incompletos e / ou existem preditores contínuos no modelo, o ajuste da estrutura de efeitos aleatórios do modelo misto pode alterar as estimativas da parte fixa do modelo , e vice versa.

Editar 2016-07-05. Dos comentários: " Você poderia elaborar ou fornecer uma citação sobre por que preditores contínuos influenciarão as estimativas da parte fixa do modelo? "

As estimativas para a parte fixa do modelo dependerão das estimativas para a parte aleatória do modelo - isto é, os componentes de variação estimados - se (mas não somente se) a variação dos preditores diferir entre os clusters. O que quase certamente será verdadeiro se algum dos preditores for contínuo (pelo menos nos dados do "mundo real" - em teoria, seria possível que isso não fosse verdade, por exemplo, em um conjunto de dados construído).

No primeiro nível, acho que tudo o que você está ignorando é o encolhimento em relação aos valores da população; " as inclinações e interceptações por sujeito do modelo de efeitos mistos estão mais próximas das estimativas populacionais do que as estimativas de mínimos quadrados dentro do sujeito " . [ref. 1] O link a seguir provavelmente também será útil ( Quais são os descritivos apropriados para os meus modelos mistos? ), Veja a resposta de Mike Lawrence).

Além disso, acho que você é um pouco azarado no seu exemplo de brinquedo, porque possui um design perfeitamente equilibrado que faz com que você tenha exatamente a mesma estimativa no caso de não haver valores ausentes.

Experimente o código a seguir, que tem o mesmo processo sem nenhum valor ausente agora:

 cat <- as.factor(sample(1:5, n*k, replace=T) ) #This should be a bit unbalanced.
 cat_i <- 1:k # intercept per kategorie
 x <- rep(1:n, k)
 sigma <- 0.2
 alpha <- 0.001
 y <- cat_i[cat] + alpha * x + rnorm(n*k, 0, sigma) 

 m1 <- lm(y ~ x)  
 m3 <- lme(y ~ x, random = ~ 1|cat, na.action = na.omit) 

 round(digits= 7,fixef(m3)) ==  round(digits=7, coef(m1)) #Not this time lad.
 #(Intercept)           x 
 #      FALSE       FALSE 

Onde agora, como seu design não é perfeitamente equilibrado, você não tem as mesmas estimativas de coeficiente.

Na verdade, se você jogar junto com seu padrão de valor ausente de uma maneira boba (por exemplo:), y[ c(1:10, 100 + 1:10, 200 + 1:10, 300 + 1:10, 400 +1:10)] <- NApara que seu design ainda esteja perfeitamente equilibrado, você obterá os mesmos coeficientes novamente.

 require(nlme)
 set.seed(128)
 n <- 100
 k <- 5
 cat <- as.factor(rep(1:k, each = n))
 cat_i <- 1:k # intercept per kategorie
 x <- rep(1:n, k)
 sigma <- 0.2
 alpha <- 0.001
 y <- cat_i[cat] + alpha * x + rnorm(n*k, 0, sigma)
 plot(x, y)

 # simulate missing data in a perfectly balanced way
 y[ c(1:10, 100 + 1:10, 200 + 1:10, 300 + 1:10, 400 +1:10)] <- NA

 m1 <- lm(y ~ x)  
 m3 <- lme(y ~ x, random = ~ 1|cat, na.action = na.omit) 

 round(digits=7,fixef(m3)) ==  round(digits=7, coef(m1)) #Look what happend now...
 #(Intercept)           x 
 #       TRUE        TRUE 

Você está marginalmente desorientado pelo design perfeito do seu experimento original. Quando você inseriu os NAs de uma maneira desequilibrada, alterou o padrão de quanta "força" os sujeitos individuais poderiam emprestar um do outro.

Em resumo, as diferenças que você vê são devidas a efeitos de retração e, mais especificamente, porque você distorceu seu design perfeitamente equilibrado original com valores ausentes não perfeitamente equilibrados.

Ref 1: Douglas Bates lme4: modelagem de efeitos mistos com R , páginas 71-72