Как запустить несколько t-тестов или ANOVA, который контролирует несколько переменных в R? - PullRequest
Новая
       113

Как запустить несколько t-тестов или ANOVA, который контролирует несколько переменных в R?

2 голосов
/ 31 марта 2020

У меня есть 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

Я хотел бы проверить, значительно ли отличаются значения Rate для четырех разных Level групп (A, B, C, D). Я знаю, что обычно я могу запустить t-тест, если на уровне есть две группы, но так как есть четыре группы, я подумал, что я мог бы или выполнить 6 t-тестов, или я мог бы запустить ANOVA, как выполняется ANOVA и интерпретировать?

Кроме того, я хотел бы посмотреть, является ли переменная P_Elderly значимым ковариатом, который объясняет некоторые отношения между Level и Rate. Если у меня есть дополнительные ковариаты, которые я хотел бы добавить позже, как бы я это сделал?

Ответы [ 2 ]

2 голосов
/ 31 марта 2020

Вы можете выбрать линейную модель, где скорость объясняется уровнем: 

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

Вы можете посмотреть на коэффициенты: 

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

Здесь A задается в качестве эталона , а коэффициенты показывают, насколько их средние значения отличаются от средних для A. Таким образом, мы можем проверить, равен ли уровень 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, где вы проверяете значимость всех ваших коэффициентов уровня в одном 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

Так что, скорее всего, средства не очень разные, если судить по вышеприведенному. Вы также можете выполнить попарный тест следующим образом: 

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

Для вашего следующего вопроса, как добавить ковариату, вам подойдет другая модель: 

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

И сравните ее с модель только с уровнем: 

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

Против, не похоже, что пожилые люди имеют большой эффект ..

1 голос
/ 31 марта 2020

Вы можете начать с теста anova, а затем использовать функцию 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"))
...