R: Calculating mean and standard error of mean for factors with lm() vs. direct calculation -edited

When dealing with data with factors R can be used to calculate the means for each group with the lm() function. This also gives the standard errors for the estimated means. But this standard error differs from what I get from a calculation by hand.

Here is an example (taken from here Predicting the difference between two groups in R )

First calculate the mean with 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  The intercept is the mean for the first group, the 4 cylindered cars. To get the means by direct calculation I use this:  with(mtcars, tapply(mpg, cyl, mean)) 4 6 8 26.66364 19.74286 15.10000  To get the standard errors for the means I calculate the sample standard variation and divide by the number of observations in each group:  with(mtcars, tapply(mpg, cyl, sd)/sqrt(summary(mtcars$cyl)) )

4         6         8
1.3597642 0.5493967 0.6842016


The direct calculation gives the same mean but the standard error is different for the 2 approaches, I had expected to get the same standard error. What is going on here? It is related to lm() fitting the mean for each group and an error term?

Edited: After Svens answer (below) I can formulate my question more concise and clearly.

For categorical data we can calculate the means of a variable for different groups is by using lm() without an intercept.

  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  We can compare this with an direct calculation of the means and their standard errors:  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


The means are exactly the same but the standard errors are different for these 2 methods (as Sven also notices). My question is why are they different and not the same?

(when editing my question, should I delete the original text or adding my edition as I did )

The difference in standard errors are because in the regression you compute a combined estimate of the variance, while in the other calculation you compute separate estimates of the variance.

• Thank you for clarifying this. I just found a very good answer for a similar question here, with a nice worked example: stats.stackexchange.com/questions/29479/… – SRJ Feb 25 '13 at 0:14
• Yes, that looks relevant. Well spotted. – Glen_b Feb 25 '13 at 0:15

The lm function does not estimate means and standard errors of the factor levels but of the contrats associated with the factor levels.

If no contrast is specified manually, treatment contrasts are used in R. This is the default for categorical data.

The factor mtcars\$cyl has three levels (4,6, and 8). By default, the first level, 4, is used as reference category. The intercept of the linear model corresponds to the mean of the dependent variable in the reference category. But the other effects result from a comparison of one factor level with the reference category. Hence, the estimate and standard error for cyl6 are related to the difference between cyl == 6 and cyl == 4. The effect cyl8 is related to the difference between cyl == 8 and cyl == 4.

If you want the lm function to calculate the means of the factor levels, you have to exclude the intercept term (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


As you can see, these estimates are identical to the means of the factor levels. But note that the standard errors of the estimates are not identical with the standard errors of the data.

By the way: Data can be aggregated easily with the aggregate function:

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

• Thanks for the response. I already know that coefficients are not the means, as I wrote the intercept is the mean of the first level, the other coefficents are the difference in mean of the other levels to that level. You also notice that with your remark "standard errors of the estimates are not identical with the standard errors of the data." Does that mean that lm() estimates the means and calculate standard errors for those estimates – SRJ Feb 22 '13 at 19:48
• Ooops, I wanted to edit that comment for clarity but did not know I only could edit for 5 minutes, can I delete a comment? I did not realize that I could get mean estimates directly by omitting the intercept, thanks for that tip. If I understand you correctly, the standard errors of the estimated means are not the same as the standard errors calculated directly from the data. Is it a different set of equations used in each case? And what are these equations? I would like to have some more details to u nderstand the difference better – SRJ Feb 22 '13 at 20:01

In addition to what Sven Hohenstein said, the mtcars data is not balanced. Usually one uses aov for lm with categorical data (which is just a wrapper for lm) which specifically says on ?aov:

aov is designed for balanced designs, and the results can be hard to interpret without balance: beware that missing values in the response(s) will likely lose the balance.

I think you can also see this on the weird correlations of the model matrix:

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


Hence, the standard errors obtained from aov (or lm) will likely be bogus (you can check this if you compare with lme or lmer standard errors.

• How would you apply lme here? – SRJ Feb 25 '13 at 0:17
• The correlations of the model matrix values are not weird. Since the constant (intercept) inherently is equal to one, there's no variation between its values. Due to this, you can't compute a correlation coefficient between a variable and the constant. – Sven Hohenstein Feb 25 '13 at 6:39
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

• This is unreadable and lacks any sort of commentary about what it means or does. Can you edit it to give more clarity? – mdewey Aug 27 at 12:57
• It is not possible to understand the answer. It does require some commentary to explain what you are doing. – Michael R. Chernick Aug 27 at 14:04