How to run a multiple t-tests or ANOVA that controls for multiple variables in R?

Question

I have 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

I would like to test if the means of Rate of four different Level groups (A,B,C,D) significantly differ. I know that I can usually run a t-test if there are two groups in level, but since there are four groups, I was thinking I could either run 6 t-tests, or I could run an ANOVA, how is an ANOVA run and interpreted?

Additionally, I would like to see if the variable P_Elderly is a significant covariate that explains some of the relationship between Level and Rate. If I have additional covariates that I would like to add later, how would I do that?

score 2 · Accepted Answer · answered Mar 31 '20 at 17:15

You can fit a linear model, where the Rate is explained by Level:

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

You can look at the coefficients:

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

Here A is set as the reference, and the coefficients indicate how much their means differ from that of A. So we can test whether Level B : D are zero, meaning one common intercept is enough for this model:

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

This is similar to an anova, where you test the significance of all your level coefficients at one go.

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

So most likely the means are not very different, going by the above. You can also do a pairwise test like this:

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

For your next question, how to add a covariate, you will fit another model:

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

And compare it against the model with only 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

Against, doesn't seem like Elderly has much of an effect..

This is great, if elderly DID have an effect, what would be the finding in the final anova? — Evan, Mar 31 '20 at 18:14
If it has an effect, you will have predicted values that fit the observed better , leading to lesser residual sum of squares (RSS bottom row), so a larger F and a lower P value — StupidWolf, Mar 31 '20 at 18:18
Great, is there anyway to judge if the effect itself is significant. What if the RSS increases but only by a small small amount? — Evan, Mar 31 '20 at 18:22
you can either do it with anova(f0,fit_full), so the difference between the two models is only P_Elderly, if you get a significant p-value, it means P_Elderly has a significant effect. You can also do summary(fit_full) from above, it gives you a p-value of each coefficient, include P_Elderly — StupidWolf, Mar 31 '20 at 18:45

score 1 · Answer 2 · answered Mar 31 '20 at 16:31

You can start by an anova test and then use TukeyHSD function to get the p values for each comparisons:

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

Does it answer your question ?

Reproducible example

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"))

This is great! Can you explain how you interpret the outputs? How do I control for a covariate? — Evan, Mar 31 '20 at 16:40
What do you mean by covariate ? To interpret the output, you can see the pvalues of each combination of levels on the last column. Then, you can interpret it as you will do with a t. test — dc37, Mar 31 '20 at 16:43
How can I test to see if P_Elderly is a mediating or confounding variable? — Evan, Mar 31 '20 at 16:44

How to run a multiple t-tests or ANOVA that controls for multiple variables in R?

2 Answers2