Почему линейная регрессия и ANOVA дают различное

22

Я пытался подогнать данные одного временного ряда (без повторов), используя регрессионную модель. Данные выглядят следующим образом:

> xx.2
          value time treat
    1  8.788269    1     0
    2  7.964719    6     0
    3  8.204051   12     0
    4  9.041368   24     0
    5  8.181555   48     0
    6  8.041419   96     0
    7  7.992336  144     0
    8  7.948658    1     1
    9  8.090211    6     1
    10 8.031459   12     1
    11 8.118308   24     1
    12 7.699051   48     1
    13 7.537120   96     1
    14 7.268570  144     1

Из-за отсутствия дубликатов я рассматриваю время как непрерывную переменную. Колонка «Лечить» показывает данные случая и контроля соответственно.

Сначала я подгоняю модель «значение = время * лечить» с помощью «lm» в R:

summary(lm(value~time*treat,data=xx.2))

Call:
lm(formula = value ~ time * treat, data = xx.2)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.50627 -0.12345  0.00296  0.04124  0.63785 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  8.493476   0.156345  54.325 1.08e-13 ***
time        -0.003748   0.002277  -1.646   0.1307    
treat       -0.411271   0.221106  -1.860   0.0925 .  
time:treat  -0.001938   0.003220  -0.602   0.5606    

Ценность времени и удовольствия не имеет значения.

Хотя с anova я получил разные результаты:

 summary(aov(value~time*treat,data=xx.2))
            Df Sum Sq Mean Sq F value Pr(>F)  
time         1 0.7726  0.7726   8.586 0.0150 *
treat        1 0.8852  0.8852   9.837 0.0106 *
time:treat   1 0.0326  0.0326   0.362 0.5606  
Residuals   10 0.8998  0.0900                 

Ценность времени и удовольствия изменилась.

С линейной регрессией, если я прав, это означает, что время и удовольствие не оказывают существенного влияния на ценность, но с ANOVA это означает, что время и удовольствие оказывают существенное влияние на ценность.

Может ли кто-нибудь объяснить мне, почему есть разница в этих двух методах, и какой использовать?

шао
источник
3
Вы можете посмотреть различные виды квадратов. В частности, я считаю, что линейная регрессия возвращает сумму квадратов типа III, тогда как anova возвращает другой вид.
принято нормальным
3
Если вы сохраните результаты lmи aovсможете проверить, что они дают одинаковые результаты; например, сравните их остатки с residualsфункцией или проверьте их коэффициенты ( $coefficientsинтервал в обоих случаях).
whuber

Ответы:

18

Подход для lm () и aov () идентичны, но отчетность отличается. T-тесты - это предельное влияние рассматриваемых переменных, учитывая наличие всех других переменных. F-тесты являются последовательными - поэтому они проверяют важность времени в присутствии ничего, кроме перехвата, обработки в присутствии ничего, кроме перехвата и времени, и взаимодействия в присутствии всего вышеперечисленного.

Предполагая, что вы заинтересованы в значении лечения, я предлагаю вам сравнить две модели, одну с, а другую без, сравнить две, поместив обе модели в anova (), и использовать этот F-тест. Это будет проверять удовольствие и взаимодействие одновременно.

Учтите следующее:

> xx.2 <- as.data.frame(matrix(c(8.788269, 1, 0,
+ 7.964719, 6, 0,
+ 8.204051, 12, 0,
+ 9.041368, 24, 0,
+ 8.181555, 48, 0,
+ 8.041419, 96, 0,
+ 7.992336, 144, 0,
+ 7.948658, 1, 1,
+ 8.090211, 6, 1,
+ 8.031459, 12, 1,
+ 8.118308, 24, 1,
+ 7.699051, 48, 1,
+ 7.537120, 96, 1,
+ 7.268570, 144, 1), byrow=T, ncol=3))
> names(xx.2) <- c("value", "time", "treat")
> 
> mod1 <- lm(value~time*treat, data=xx.2)
> anova(mod1)
Analysis of Variance Table

Response: value
           Df  Sum Sq Mean Sq F value  Pr(>F)  
time        1 0.77259 0.77259  8.5858 0.01504 *
treat       1 0.88520 0.88520  9.8372 0.01057 *
time:treat  1 0.03260 0.03260  0.3623 0.56064  
Residuals  10 0.89985 0.08998                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1 
> mod2 <- aov(value~time*treat, data=xx.2)
> anova(mod2)
Analysis of Variance Table

Response: value
           Df  Sum Sq Mean Sq F value  Pr(>F)  
time        1 0.77259 0.77259  8.5858 0.01504 *
treat       1 0.88520 0.88520  9.8372 0.01057 *
time:treat  1 0.03260 0.03260  0.3623 0.56064  
Residuals  10 0.89985 0.08998                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1 
> summary(mod2)
            Df Sum Sq Mean Sq F value Pr(>F)  
time         1 0.7726  0.7726   8.586 0.0150 *
treat        1 0.8852  0.8852   9.837 0.0106 *
time:treat   1 0.0326  0.0326   0.362 0.5606  
Residuals   10 0.8998  0.0900                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1 
> summary(mod1)

Call:
lm(formula = value ~ time * treat, data = xx.2)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.50627 -0.12345  0.00296  0.04124  0.63785 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  8.493476   0.156345  54.325 1.08e-13 ***
time        -0.003748   0.002277  -1.646   0.1307    
treat       -0.411271   0.221106  -1.860   0.0925 .  
time:treat  -0.001938   0.003220  -0.602   0.5606    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1 

Residual standard error: 0.3 on 10 degrees of freedom
Multiple R-squared: 0.6526,     Adjusted R-squared: 0.5484 
F-statistic: 6.262 on 3 and 10 DF,  p-value: 0.01154 
Питер Эллис
источник
Спасибо за подробное объяснение, оно напоминает мне ANCOVA (анализ ковариации). Первым шагом ANCOVA является проверка взаимодействия между категориальным фактором и ковариатой, чтобы увидеть, имеют ли они одинаковый наклон для обоих условий. Это очень похоже на то, что я сделал здесь. В ANCOVA он дает одинаковое значение для взаимодействия в t-тесте и F-тесте, поскольку взаимодействие является последним членом в aov.
Шао
17

Ответ Питера Эллиса превосходен, но есть еще одно замечание. -test статистика (и его -значение) является испытанием ли . -test на распечатке , является ли значительно снижает добавленный переменная остаточную сумму квадратов.p β = 0 FTпβзнак равно0Fanova()

-test заказать независимый, в то время как -test нет. Отсюда и предложение Питера попробовать переменные в разных порядках. Также возможно, что переменные, значимые в одном тесте, могут не быть значимыми в другом (и наоборот).FtF

Я чувствую (и другие участники могут поправить меня), что когда вы пытаетесь предсказать явления (как в системном приложении), вы больше всего заинтересованы в уменьшении дисперсии с наименьшим количеством предикторов и, следовательно, хотите получить anova()результаты. Однако, если вы пытаетесь установить предельное влияние на , вас больше всего заинтересует значение вашей конкретной , и все остальные переменные будут просто контролировать альтернативные объяснения, которые ваши коллеги-рецензенты попытаются найти.y βИксYβ

gregmacfarlane
источник
2

Приведенные выше два ответа великолепны, но я подумал, что добавлю немного больше. Другой кусок информации можно почерпнуть отсюда .

Когда вы сообщаете о lm()результатах с помощью термина взаимодействия, вы говорите что-то вроде: «лечения 1 отличается от лечения 0 (бета! = 0, р = 0,0925), когда время установлено на базовое значение 1 ». Принимая во внимание, что anova()результаты ( как упомянуто ранее ) игнорируют любые другие переменные и касаются только различий.

Вы можете доказать это, удалив член взаимодействия и используя простую модель только с двумя основными эффектами ( m1 ):

> m1 = lm(value~time+treat,data=dat)
> summary(m1)

Call:
lm(formula = value ~ time + treat, data = dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.54627 -0.10533 -0.04574  0.11975  0.61528 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  8.539293   0.132545  64.426 1.56e-15 ***
time        -0.004717   0.001562  -3.019  0.01168 *  
treat       -0.502906   0.155626  -3.232  0.00799 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2911 on 11 degrees of freedom
Multiple R-squared:   0.64, Adjusted R-squared:  0.5746 
F-statistic: 9.778 on 2 and 11 DF,  p-value: 0.003627

> anova(m1)
Analysis of Variance Table

Response: value
          Df  Sum Sq Mean Sq F value   Pr(>F)   
time       1 0.77259 0.77259  9.1142 0.011677 * 
treat      1 0.88520 0.88520 10.4426 0.007994 **
Residuals 11 0.93245 0.08477                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

В этом случае мы видим, что сообщаемые значения p одинаковы; это потому, что в случае этой более простой модели,

Константино
источник
Этот ответ, к сожалению, выглядит незаконченным. Еще +1 за ссылку и за упоминание того, что эффект связан с разными схемами кодирования.
говорит амеба, восстанови Монику
2
Следует также добавить это summary(lm)и anova(lm)не всегда будет давать идентичный результат, если отсутствует термин взаимодействия. Так уж получилось, что в этих данных timeи treatони ортогональны и поэтому суммы квадратов типа I (последовательные) и III (предельные) дают одинаковые результаты.
говорит амеба: восстанови Монику
2
  • Разница связана с типом парных сравнений каскадных моделей.
  • Кроме того, функция aov () имеет проблему с тем, как она выбирает степени свободы. Кажется, что смешаны два понятия: 1) сумма квадратов из пошаговых сравнений, 2) степени свободы от общей картины.

ПРОБЛЕМА РЕПРОДУКЦИИ

> data <- list(value = c (8.788269,7.964719,8.204051,9.041368,8.181555,8.0414149,7.992336,7.948658,8.090211,8.031459,8.118308,7.699051,7.537120,7.268570), time = c(1,6,12,24,48,96,144,1,6,12,24,48,96,144), treat = c(0,0,0,0,0,0,0,1,1,1,1,1,1,1) )
> summary( lm(value ~ treat*time, data=data) )
> summary( aov(value ~ 1 + treat + time + I(treat*time),data=data) )

НЕКОТОРЫЕ МОДЕЛИ, ИСПОЛЬЗУЕМЫЕ В ОБЪЯСНЕНИИ

#all linear models used in the explanation below
> model_0                      <- lm(value ~ 1, data)
> model_time                   <- lm(value ~ 1 + time, data)
> model_treat                  <- lm(value ~ 1 + treat, data)
> model_interaction            <- lm(value ~ 1 + I(treat*time), data)
> model_treat_time             <- lm(value ~ 1 + treat + time, data)
> model_treat_interaction      <- lm(value ~ 1 + treat + I(treat*time), data)
> model_time_interaction       <- lm(value ~ 1 + time + I(treat*time), data)
> model_treat_time_interaction <- lm(value ~ 1 + time + treat + I(treat*time), data)

КАК LM T_TEST работает и имеет отношение к F-TEST

# the t-test with the estimator and it's variance, mean square error, is
# related to the F test of pairwise comparison of models by dropping 1
# model parameter

> anova(model_treat_time_interaction, model_time_interaction)

Analysis of Variance Table

Model 1: value ~ 1 + time + treat + I(treat * time)
Model 2: value ~ 1 + time + I(treat * time)
  Res.Df     RSS Df Sum of Sq      F  Pr(>F)  
1     10 0.89985                              
2     11 1.21118 -1  -0.31133 3.4598 0.09251 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

> anova(model_treat_time_interaction, model_treat_interaction)

Analysis of Variance Table

Model 1: value ~ 1 + time + treat + I(treat * time)
Model 2: value ~ 1 + treat + I(treat * time)
  Res.Df     RSS Df Sum of Sq      F Pr(>F)
1     10 0.89985                           
2     11 1.14374 -1   -0.2439 2.7104 0.1307

> anova(model_treat_time_interaction, model_treat_time)

Analysis of Variance Table

Model 1: value ~ 1 + time + treat + I(treat * time)
Model 2: value ~ 1 + treat + time
  Res.Df     RSS Df Sum of Sq      F Pr(>F)
1     10 0.89985                           
2     11 0.93245 -1 -0.032599 0.3623 0.5606

> # which is the same as
> drop1(model_treat_time_interaction, scope  = ~time+treat+I(treat*time), test="F")

Single term deletions

Model:
value ~ 1 + time + treat + I(treat * time)
                Df Sum of Sq     RSS     AIC F value  Pr(>F)  
<none>                       0.89985 -30.424                  
time             1  0.243896 1.14374 -29.067  2.7104 0.13072  
treat            1  0.311333 1.21118 -28.264  3.4598 0.09251 .
I(treat * time)  1  0.032599 0.93245 -31.926  0.3623 0.56064  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

КАК AOV РАБОТАЕТ И ВЫБИРАЕТ DF В F-ТЕСТАХ

> #the aov function makes stepwise additions/drops
> 
> #first the time, then treat, then the interaction
> anova(model_0, model_time)

Analysis of Variance Table

Model 1: value ~ 1
Model 2: value ~ 1 + time
  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1     13 2.5902                              
2     12 1.8176  1    0.7726 5.1006 0.04333 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

> anova(model_time, model_treat_time)

Analysis of Variance Table

Model 1: value ~ 1 + time
Model 2: value ~ 1 + treat + time
  Res.Df     RSS Df Sum of Sq      F   Pr(>F)   
1     12 1.81764                                
2     11 0.93245  1    0.8852 10.443 0.007994 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

> anova(model_treat_time, model_treat_time_interaction)

Analysis of Variance Table

Model 1: value ~ 1 + treat + time
Model 2: value ~ 1 + time + treat + I(treat * time)
  Res.Df     RSS Df Sum of Sq      F Pr(>F)
1     11 0.93245                           
2     10 0.89985  1  0.032599 0.3623 0.5606

> 
> # note that the sum of squares for within model variation is the same
> # but the F values and p-values are not the same because the aov 
> # function somehow chooses to use the degrees of freedom in the 
> # complete model in all stepwise changes
>

ВАЖНАЯ ЗАМЕТКА

> # Although the p and F values do not exactly match, it is this effect
> # of order and selection of cascading or not in model comparisons. 
> # An important note to make is that the comparisons are made by 
> # stepwise additions and changing the order of variables has an 
> # influence on the outcome!
>
> # Additional note changing the order of 'treat' and 'time' has no 
> # effect because they are not correlated

> summary( aov(value ~ 1 + treat + time +I(treat*time), data=data) )

        Df Sum Sq Mean Sq F value Pr(>F)  
treat            1 0.8852  0.8852   9.837 0.0106 *
time             1 0.7726  0.7726   8.586 0.0150 *
I(treat * time)  1 0.0326  0.0326   0.362 0.5606  
Residuals       10 0.8998  0.0900                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

> summary( aov(value ~ 1 + I(treat*time) + treat + time, data=data) )

                Df Sum Sq Mean Sq F value  Pr(>F)   
I(treat * time)  1 1.3144  1.3144  14.606 0.00336 **
treat            1 0.1321  0.1321   1.469 0.25343   
time             1 0.2439  0.2439   2.710 0.13072   
Residuals       10 0.8998  0.0900                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

> # This is an often forgotten quirck 
> # best is to use manual comparisons such that you know
> # and understand your hypotheses
> # (which is often forgotten in the click and
> #     point anova modelling tools)
> #
> # anova(model1, model2) 
> #     or use 
> # stepAIC from the MASS library
Секст Эмпирик
источник