R: Calcul de la moyenne et de l'erreur standard de la moyenne pour les facteurs avec lm () vs calcul direct - édité

Lorsqu'il s'agit de données avec des facteurs, R peut être utilisé pour calculer les moyennes pour chaque groupe avec la fonction lm (). Cela donne également les erreurs types pour les moyennes estimées. Mais cette erreur standard diffère de ce que j'obtiens d'un calcul à la main.

Voici un exemple (tiré d'ici Prédire la différence entre deux groupes dans R )

Calculez d'abord la moyenne avec 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

L'ordonnée à l'origine est la moyenne du premier groupe, les voitures à 4 cylindres. Pour obtenir les moyens par calcul direct, j'utilise ceci:

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

         4        6        8 
    26.66364 19.74286 15.10000 

Pour obtenir les erreurs standard pour les moyennes, je calcule la variation standard de l'échantillon et je divise par le nombre d'observations dans chaque groupe:

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

         4         6         8 
   1.3597642 0.5493967 0.6842016 

Le calcul direct donne la même moyenne mais l'erreur standard est différente pour les 2 approches, je m'attendais à obtenir la même erreur standard. Qu'est-ce qui se passe ici? Il est lié à lm () ajustant la moyenne pour chaque groupe et un terme d'erreur?

Modifié: Après la réponse de Svens (ci-dessous), je peux formuler ma question de manière plus concise et plus claire.

Pour les données catégorielles, nous pouvons calculer la moyenne d'une variable pour différents groupes en utilisant lm () sans interception.

  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

On peut comparer cela avec un calcul direct des moyennes et de leurs erreurs types:

  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 

Les moyennes sont exactement les mêmes mais les erreurs standard sont différentes pour ces 2 méthodes (comme Sven le remarque également). Ma question est pourquoi sont-ils différents et pas les mêmes?

(lors de la modification de ma question, dois-je supprimer le texte d'origine ou ajouter mon édition comme je l'ai fait)

La différence dans les erreurs standard est due au fait que dans la régression, vous calculez une estimation combinée de la variance, tandis que dans l'autre calcul, vous calculez des estimations distinctes de la variance.

La lmfonction n'évalue pas les moyennes et les erreurs-types des niveaux de facteurs mais des contrats associés aux niveaux de facteurs.

Si aucun contraste n'est spécifié manuellement, les contrastes de traitement sont utilisés dans R. Il s'agit de la valeur par défaut pour les données catégorielles.

Le facteur mtcars$cyla trois niveaux (4,6 et 8). Par défaut, le premier niveau, 4, est utilisé comme catégorie de référence. L'ordonnée à l'origine du modèle linéaire correspond à la moyenne de la variable dépendante dans la catégorie de référence. Mais les autres effets résultent d'une comparaison d'un niveau de facteur avec la catégorie de référence. Par conséquent, l'estimation et l'erreur-type de cyl6sont liées à la différence entre cyl == 6et cyl == 4. L'effet cyl8est lié à la différence entre cyl == 8et cyl == 4.

Si vous voulez que la lmfonction calcule la moyenne des niveaux de facteur, vous devez exclure le terme d'interception ( 0 + ...):

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

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

Residuals:
    Min      1Q  Median      3Q     Max 
-5.2636 -1.8357  0.0286  1.3893  7.2364 

Coefficients:
                Estimate Std. Error t value 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 

Residual standard error: 3.223 on 29 degrees of freedom
Multiple R-squared: 0.9785, Adjusted R-squared: 0.9763 
F-statistic: 440.9 on 3 and 29 DF,  p-value: < 2.2e-16 

Comme vous pouvez le voir, ces estimations sont identiques aux moyennes des niveaux des facteurs. Mais notez que les erreurs-types des estimations ne sont pas identiques aux erreurs-types des données.

Soit dit en passant: les données peuvent être agrégées facilement avec la aggregatefonction:

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

En plus de ce que Sven Hohenstein a dit, les mtcarsdonnées ne sont pas équilibrées . Habituellement, on utilise aovpour lm avec des données catégorielles (qui sont juste un wrapper pour lm) qui dit spécifiquement sur ?aov:

aov est conçu pour des conceptions équilibrées et les résultats peuvent être difficiles à interpréter sans équilibre: prenez garde que les valeurs manquantes dans la ou les réponses perdront probablement l'équilibre.

Je pense que vous pouvez également le voir sur les corrélations étranges de la matrice du modèle:

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
Warning message:
In cor(mf) : the standard deviation is zero

Par conséquent, les erreurs standard obtenues à partir de aov(ou lm) seront probablement fausses (vous pouvez le vérifier si vous comparez avec lmeou avec lmerdes erreurs standard.

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