Agregar efectos aleatorios influye en las estimaciones de coeficientes

Siempre me han enseñado que los efectos aleatorios solo influyen en la varianza (error), y que los efectos fijos solo influyen en la media. Pero he encontrado un ejemplo en el que los efectos aleatorios también influyen en la media: la estimación del 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)

Puede ver que el coeficiente estimado para el xmodelo m1es -0.013780, mientras que para el modelo m3es 0.0011713, ambos significativamente diferentes de cero.

Tenga en cuenta que cuando elimino la línea que simula datos faltantes, los resultados son los mismos (es una matriz completa).

¿Porqué es eso?

PD: tenga en cuenta que no soy un estadístico profesional, por lo que si está a punto de responder con muchas matemáticas, haga también un resumen simple para tontos :-)

"Siempre me han enseñado que los efectos aleatorios solo influyen en la varianza (error), y que los efectos fijos solo influyen en la media".

Como ha descubierto, esto solo es cierto para conjuntos de datos equilibrados, completos (es decir, sin datos faltantes) sin predictores continuos. En otras palabras, para los tipos de datos / modelos discutidos en los textos clásicos de ANOVA. En estas circunstancias ideales, los efectos fijos y los efectos aleatorios pueden estimarse independientemente uno del otro.

Cuando estas condiciones no se cumplen (como muy a menudo no lo hacen en el "mundo real"), los efectos fijos y aleatorios no son independientes. Como comentario interesante, esta es la razón por la cual los modelos mixtos "modernos" se estiman utilizando métodos de optimización iterativos, en lugar de resolverse exactamente con un poco de álgebra matricial como en el caso ANOVA mixto clásico: para estimar los efectos fijos, tenemos que Conozca los efectos aleatorios, pero para estimar los efectos aleatorios, ¡tenemos que conocer los efectos fijos! Más relevante para la presente pregunta, esto también significa que cuando los datos están desequilibrados / incompletos y / o hay predictores continuos en el modelo, el ajuste de la estructura de efectos aleatorios del modelo mixto puede alterar las estimaciones de la parte fija del modelo. , y viceversa.

Editar 2016-07-05. De los comentarios: " ¿Podría elaborar o proporcionar una cita sobre por qué los predictores continuos influirán en las estimaciones de la parte fija del modelo? "

Las estimaciones para la parte fija del modelo dependerán de las estimaciones para la parte aleatoria del modelo, es decir, los componentes de varianza estimados, si (pero no solo si) la varianza de los predictores difiere entre los grupos. Lo que seguramente será cierto si alguno de los predictores es continuo (al menos en los datos del "mundo real", en teoría sería posible que esto no sea cierto, por ejemplo, en un conjunto de datos construido).

En el primer nivel, creo que todos ustedes ignoran la contracción hacia los valores de la población; " las pendientes e intersecciones por sujeto del modelo de efectos mixtos están más cerca de las estimaciones de población que las estimaciones de mínimos cuadrados dentro del sujeto " . [ref. 1] El siguiente enlace probablemente también sea de ayuda ( ¿Cuáles son los descriptivos adecuados para mis modelos mixtos? ), Vea la respuesta de Mike Lawrence).

Además, creo que eres un poco desafortunado en tu ejemplo de juguete porque tienes un diseño perfectamente equilibrado que hace que tengas exactamente la misma estimación en el caso de que no falten valores.

Pruebe el siguiente código que tiene el mismo proceso sin ningún valor perdido ahora:

 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 

Donde ahora, debido a que su diseño no está perfectamente equilibrado, no tiene las mismas estimaciones de coeficientes.

En realidad, si juegas junto con tu patrón de valor perdido de una manera tonta (por ejemplo:) y[ c(1:10, 100 + 1:10, 200 + 1:10, 300 + 1:10, 400 +1:10)] <- NApara que tu diseño aún esté perfectamente equilibrado, obtendrás los mismos coeficientes nuevamente.

 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 

Usted está marginalmente equivocado por el diseño perfecto de su experimento original. Cuando insertó los NA en una distancia no equilibrada, cambió el patrón de cuánta "fuerza" podían prestarse los sujetos individuales entre sí.

En resumen, las diferencias que ve se deben a los efectos de contracción y más específicamente porque distorsionó su diseño perfectamente equilibrado original con valores perdidos no perfectamente equilibrados.

Ref. 1: Douglas Bates lme4: Modelado de efectos mixtos con R , páginas 71-72