如何运行多个 t 检验或 ANOVA 来控制 R 中的多个变量？

Question

我有 df1：

    Rate       Dogs        MHI_2018  Points      Level AGE65_MORE P_Elderly
1  0.10791173  0.00000000    59338 236.4064         C       8653  15.56267
2  0.06880040  0.00000000    57588 229.4343         C      44571  20.44335
3  0.08644537  0.00000000    50412 200.8446         C      10548  18.23651
4  0.29591635  0.00000000    29267 116.6016         A       1661  16.38390
5  0.05081301  0.00000000    37365 148.8645         B       3995  20.29980
6  0.02625200  0.00000000    45400 180.8765         D      20247  17.71748
7  0.80321285  0.02974862    39917 159.0319         D       6562  19.52105
8  0.07682852  0.00000000    42132 167.8566         D       5980  22.97173
9  0.18118814  0.00000000    47547 189.4303         B       7411  16.78482
10 0.07787555  0.00000000    39907 158.9920         B       2953  22.99665
11 0.15065913  0.00000000    39201 156.1793         C       2751  20.72316
12 0.33362247  0.00000000    46495 185.2390         B       2915  19.45019
13 0.03652168  0.00000000    49055 195.4382         B      10914  19.92988
14 0.27998133  0.00000000    42423 169.0159         A       2481  23.15446
15 0.05407451  0.00000000    40203 160.1713         A       7790  21.06202
16 0.07233796  0.00000000    39057 155.6056         A       2629  19.01765
17 0.08389061  0.00000000    45796 182.4542         B      15446  18.51106
18 0.05220569  0.00000000    34035 135.5976         B       6921  18.06578
19 0.05603418  0.00000000    39491 157.3347         B      12322  17.26133
20 0.15875536  0.00000000    60367 240.5060         C      12400  15.14282

我想测试四个不同Level 组（A、B、C、D）的Rate 的平均值是否显着不同。我知道如果有两组水平，我通常可以运行 t 检验，但是由于有四个组，我想我可以运行 6 个 t 检验，或者我可以运行 ANOVA，ANOVA 运行情况如何并解释？

此外，我想看看变量P_Elderly 是否是一个重要的协变量，它解释了Level 和Rate 之间的一些关系。如果我想稍后添加其他协变量，我该怎么做？

Answer 1

您可以拟合线性模型，其中 Rate 由 Level 解释：

fit0 = lm(Rate ~ Level,data=df)

你可以看看系数：

coefs = coefficients(fit0)
coefs
(Intercept)      LevelB      LevelC      LevelD 
 0.17557754 -0.06655862 -0.06106314  0.12652025

这里以 A 为参考，系数表示它们的均值与 A 的差异有多大。因此我们可以测试 Level B : D 是否为零，这意味着该模型的一个公共截距就足够了：

library(car)
linearHypothesis(fit0,names(coefs)[-1],test="F")
Linear hypothesis test

Hypothesis:
LevelB = 0
LevelC = 0
LevelD = 0

Model 1: restricted model
Model 2: Rate ~ Level

  Res.Df     RSS Df Sum of Sq      F Pr(>F)
1     19 0.59848                           
2     16 0.50688  3  0.091608 0.9639 0.4338

这类似于方差分析，您可以一次性测试所有水平系数的显着性。

anova(fit0)
Analysis of Variance Table

Response: Rate
          Df  Sum Sq  Mean Sq F value Pr(>F)
Level      3 0.09161 0.030536  0.9639 0.4338
Residuals 16 0.50688 0.031680

所以很可能这些方法并没有太大的不同，按照上述方法。您也可以像这样进行成对测试：

library(multcomp)
summary(glht(fit0,linfct = mcp(Level = "Tukey")))

对于您的下一个问题，如何添加协变量，您将拟合另一个模型：

fit_full = lm(Rate ~ Level+P_Elderly,data=df)

并将其与仅具有 Level 的模型进行比较：

anova(f0,fit_full)

Analysis of Variance Table

Model 1: Rate ~ Level
Model 2: Rate ~ Level + P_Elderly
  Res.Df     RSS Df Sum of Sq      F Pr(>F)
1     16 0.50688                           
2     15 0.50150  1 0.0053721 0.1607 0.6942

反对，老者好像没有多大作用..

Answer 2

您可以从方差分析开始，然后使用TukeyHSD 函数获取每次比较的 p 值：

AOV <- aov(Rate~Level, data = df)

Call:
   aov(formula = Rate ~ Level, data = df)

Terms:
                    Level Residuals
Sum of Squares  0.0916076 0.5068768
Deg. of Freedom         3        16

Residual standard error: 0.1779882
Estimated effects may be unbalanced

TukeyHSD(AOV)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = Rate ~ Level, data = df)

$Level
            diff        lwr       upr     p adj
B-A -0.066558621 -0.3783957 0.2452784 0.9272012
C-A -0.061063140 -0.4026635 0.2805372 0.9551663
D-A  0.126520253 -0.2624089 0.5154494 0.7890519
C-B  0.005495482 -0.2848090 0.2958000 0.9999404
D-B  0.193078874 -0.1516699 0.5378277 0.4049948
D-C  0.187583392 -0.1843040 0.5594708 0.4923479

它回答了你的问题吗？

---

可重现的例子

structure(list(Row = 1:20, Rate = c(0.10791173, 0.0688004, 0.08644537, 
0.29591635, 0.05081301, 0.026252, 0.80321285, 0.07682852, 0.18118814, 
0.07787555, 0.15065913, 0.33362247, 0.03652168, 0.27998133, 0.05407451, 
0.07233796, 0.08389061, 0.05220569, 0.05603418, 0.15875536), 
    Dogs = c(0, 0, 0, 0, 0, 0, 0.02974862, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0), MHI_2018 = c(59338L, 57588L, 50412L, 
    29267L, 37365L, 45400L, 39917L, 42132L, 47547L, 39907L, 39201L, 
    46495L, 49055L, 42423L, 40203L, 39057L, 45796L, 34035L, 39491L, 
    60367L), Points = c(236.4064, 229.4343, 200.8446, 116.6016, 
    148.8645, 180.8765, 159.0319, 167.8566, 189.4303, 158.992, 
    156.1793, 185.239, 195.4382, 169.0159, 160.1713, 155.6056, 
    182.4542, 135.5976, 157.3347, 240.506), Level = c("C", "C", 
    "C", "A", "B", "D", "D", "D", "B", "B", "C", "B", "B", "A", 
    "A", "A", "B", "B", "B", "C"), AGE65_MORE = c(8653L, 44571L, 
    10548L, 1661L, 3995L, 20247L, 6562L, 5980L, 7411L, 2953L, 
    2751L, 2915L, 10914L, 2481L, 7790L, 2629L, 15446L, 6921L, 
    12322L, 12400L), P_Elderly = c(15.56267, 20.44335, 18.23651, 
    16.3839, 20.2998, 17.71748, 19.52105, 22.97173, 16.78482, 
    22.99665, 20.72316, 19.45019, 19.92988, 23.15446, 21.06202, 
    19.01765, 18.51106, 18.06578, 17.26133, 15.14282)), row.names = c(NA, 
-20L), class = c("data.table", "data.frame"))