R：计算 lm() 与直接计算的因子的平均值和平均值的标准误差 - 编辑

机器算法验证 r 分类数据 意思是 流明
2022-04-06 15:21:45

在处理具有因子 R 的数据时，可使用 lm() 函数计算每个组的均值。这也给出了估计均值的标准误差。但是这个标准误差与我手工计算得到的不同。

这是一个示例（取自此处预测 R 中两组之间的差异）

首先用 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

截距是第一组 4 缸汽车的平均值。为了通过直接计算获得手段，我使用这个：

  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

直接计算给出了相同的平均值，但两种方法的标准误差不同，我预计会得到相同的标准误差。这里发生了什么？它与 lm() 拟合每个组的平均值和错误项有关吗？

编辑： 在斯文斯回答（下）之后，我可以更简洁明了地提出我的问题。

对于分类数据，我们可以通过使用不带截距的 lm() 来计算不同组的变量均值。

  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

我们可以将其与直接计算均值及其标准误差进行比较：

  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

这两种方法的方法完全相同，但标准误差不同（Sven 也注意到）。我的问题是为什么它们不同而不相同？

（编辑我的问题时，我应该删除原始文本还是像我一样添加我的版本）

4个回答

标准误差的差异是因为在回归中您计算方差的组合估计，而在其他计算中您计算方差的单独估计。

该lm函数不估计因子水平的均值和标准误，而是估计与因子水平相关的对比。

如果未手动指定对比，则在 R 中使用处理对比。这是分类数据的默认设置。

该因子mtcars$cyl具有三个水平（4、6 和 8）。默认情况下，第一个级别 4 用作参考类别。线性模型的截距对应于参考类别中因变量的平均值。但其他影响来自一个因素水平与参考类别的比较。因此， 的估计和标准误差与和cyl6之间的差异有关。效果与和之间的差异有关。cyl == 6cyl == 4cyl8cyl == 8cyl == 4

如果您希望lm函数计算因子水平的均值，则必须排除截距项 ( 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

如您所见，这些估计值与因子水平的均值相同。但请注意，估计的标准误与数据的标准误不同。

顺便说一句：可以使用以下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

除了 Sven Hohenstein 所说的，mtcars数据不平衡。通常一个aov用于 lm 的分类数据（这只是一个包装器lm），它特别说明?aov：

aov 是为平衡设计而设计的，如果没有平衡，结果可能很难解释：请注意，响应中的缺失值可能会失去平衡。

我认为您也可以在模型矩阵的奇怪相关性上看到这一点：

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

因此，从aov(or lm) 获得的标准错误可能是虚假的（如果您与lmeorlmer标准错误进行比较，您可以检查这一点。

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
其它你可能感兴趣的问题
上一篇什么叫多重测试？ 下一篇在 R 中使用 ARMA 模型进行过滤