R: Calculando la media y el error estándar de la media para factores con lm() vs. cálculo directo -editado

Question

Cuando se trata de datos con factores, R se puede usar para calcular las medias de cada grupo con la función lm(). Esto también proporciona los errores estándar para las medias estimadas. Pero este error estándar difiere de lo que obtengo de un cálculo manual.

Aquí hay un ejemplo (tomado de aquí Predicting the difference between two groups in R)

Primero calcula la media con lm():

mtcars$cyl <- factor(mtcars$cyl)
mylm <- lm(mpg ~ cyl, data = mtcars)
summary(mylm)$coef

           Estimate Std. Error   t value     Pr(>|t|)
(Intercept) 26.663636  0.9718008 27.437347 2.688358e-22
cyl6        -6.920779  1.5583482 -4.441099 1.194696e-04
cyl8       -11.563636  1.2986235 -8.904534 8.568209e-10

La intersección es la media para el primer grupo, los autos de 4 cilindros. Para obtener las medias por cálculo directo uso esto:

with(mtcars, tapply(mpg, cyl, mean))

         4        6        8 
    26.66364 19.74286 15.10000 

Para obtener los errores estándar para las medias calculo la desviación estándar de la muestra y divido por el número de observaciones en cada grupo:

with(mtcars, tapply(mpg, cyl, sd)/sqrt(summary(mtcars$cyl)))

         4         6         8 
   1.3597642 0.5493967 0.6842016 

El cálculo directo da la misma media pero el error estándar es diferente para los 2 enfoques, esperaba obtener el mismo error estándar. ¿Qué está pasando aquí? ¿Está relacionado con lm() ajustando la media para cada grupo y un término de error?

Editado: Después de la respuesta de Sven (abajo) puedo formular mi pregunta de forma más concisa y clara.

Para datos categóricos podemos calcular las medias de una variable para diferentes grupos utilizando lm() sin una intersección.

mtcars$cyl <- factor(mtcars$cyl)
mylm <- lm(mpg ~ cyl, data = mtcars)
summary(mylm)$coef

     Estimate Std. Error
cyl4 26.66364  0.9718008
cyl6 19.74286  1.2182168
cyl8 15.10000  0.8614094

Podemos comparar esto con un cálculo directo de las medias y sus errores estándar:

with(mtcars, tapply(mpg, cyl, mean))

         4        6        8 
    26.66364 19.74286 15.10000 

with(mtcars, tapply(mpg, cyl, sd)/sqrt(summary(mtcars$cyl)))

         4         6         8 
   1.3597642 0.5493967 0.6842016 

Las medias son exactamente las mismas pero los errores estándar son diferentes para estos 2 métodos (como también señala Sven). Mi pregunta es ¿por qué son diferentes y no iguales?

(al editar mi pregunta, ¿debo borrar el texto original o agregar mi edición como lo hice?)

Answer 1

La diferencia en los errores estándar se debe a que en la regresión se calcula una estimación combinada de la varianza, mientras que en el otro cálculo se calculan estimaciones separadas de la varianza.

Answer 2

La función lm no estima las medias y errores estándar de los niveles del factor, sino de los contrastes asociados con los niveles del factor.

Si no se especifica manualmente ningún contraste, se utilizan los contrastes de tratamiento en R. Este es el valor predeterminado para los datos categóricos.

El factor mtcars$cyl tiene tres niveles (4, 6 y 8). Por defecto, el primer nivel, 4, se utiliza como categoría de referencia. La intersección del modelo lineal corresponde a la media de la variable dependiente en la categoría de referencia. Pero los otros efectos resultan de una comparación de un nivel del factor con la categoría de referencia. Por lo tanto, la estimación y el error estándar para cyl6 están relacionados con la diferencia entre cyl == 6 y cyl == 4. El efecto cyl8 está relacionado con la diferencia entre cyl == 8 y cyl == 4.

Si deseas que la función lm calcule las medias de los niveles del factor, debes excluir el término de intersección (0 + ...):

summary(lm(mpg ~ 0 + as.factor(cyl), mtcars))

Llamada:
lm(formula = mpg ~ 0 + as.factor(cyl), data = mtcars)

Residuos:
    Min      1Q  Mediana      3Q     Máx 
-5.2636 -1.8357  0.0286  1.3893  7.2364 

Coeficientes:
                Estimado Error estándar valor t Pr(>|t|)    
as.factor(cyl)4  26.6636     0.9718   27.44  < 2e-16 ***
as.factor(cyl)6  19.7429     1.2182   16.21 4.49e-16 ***
as.factor(cyl)8  15.1000     0.8614   17.53  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Error estándar residual: 3.223 on 29 degrees of freedom
R-cuadrado múltiple: 0.9785, R-cuadrado ajustado: 0.9763 
Estadístico F: 440.9 on 3 and 29 DF,  valor p: < 2.2e-16 

Como puedes ver, estas estimaciones son idénticas a las medias de los niveles del factor. Pero ten en cuenta que los errores estándar de las estimaciones no son idénticos a los errores estándar de los datos.

Por cierto: Los datos se pueden agregar fácilmente con la función aggregate:

aggregate(mpg ~ cyl, mtcars, function(x) c(M = mean(x), SE = sd(x)/sqrt(length(x))))

  cyl      mpg.M     mpg.SE
1   4 26.6636364  1.3597642
2   6 19.7428571  0.5493967
3   8 15.1000000  0.6842016

Answer 3

Además de lo que dijo Sven Hohenstein, los datos mtcars no están balanceados. Normalmente se usa aov para lm con datos categóricos (que es simplemente un envoltorio para lm) que específicamente dice en ?aov:

aov está diseñado para diseños balanceados, y los resultados pueden ser difíciles de interpretar sin equilibrio: ten en cuenta que los valores faltantes en la(s) respuesta(s) probablemente perderán el equilibrio.

Creo que también puedes ver esto en las correlaciones extrañas de la matriz del modelo:

mf <- model.matrix(mpg ~ cyl, data = mtcars)
cor(mf)
            (Intercept)       cyl6       cyl8
(Intercept)           1         NA         NA
cyl6                 NA  1.0000000 -0.4666667
cyl8                 NA -0.4666667  1.0000000
Mensaje de advertencia:
In cor(mf) : la desviación estándar es cero

Por lo tanto, es probable que los errores estándar obtenidos de aov (o lm) sean incorrectos (puedes verificar esto si comparas con los errores estándar de lme o lmer).

Answer 4

Y = matrix(0,5,6)
Y[1,] = c(1250, 980, 1800, 2040, 1000, 1180)
Y[2,] = c(1700, 3080,1700,2820,5760,3480)
Y[3,] = c(2050,3560,2800,1600,4200,2650)
Y[4,] = c(4690,4370,4800,9070,3770,5250)
Y[5,] = c(7150,3480,5010,4810,8740,7260)

n = ncol(Y)
R = rowMeans(Y)
M = mean(R)

s = mean(apply(Y,1,var))

v = var(R)  -s/n

#z = n/(n+(E(s2)/var(m)))
Q = 6/(6+(s/v))
t = Q*R[1] + (1-Z)*M