> 方差分析

## 数据


```r
head(warpbreaks)
```

```
##   breaks wool tension
## 1     26    A       L
## 2     30    A       L
## 3     54    A       L
## 4     25    A       L
## 5     70    A       L
## 6     52    A       L
```

两种简单的方差分析运行方法：


```r
tens.aov <- aov(breaks ~ tension, data = warpbreaks)
print(tens.aov)
```

```
## Call:
##    aov(formula = breaks ~ tension, data = warpbreaks)
## 
## Terms:
##                  tension Residuals
## Sum of Squares  2034.259  7198.556
## Deg. of Freedom        2        51
## 
## Residual standard error: 11.88058
## Estimated effects may be unbalanced
```

```r
summary(tens.aov)
```

```
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## tension      2   2034  1017.1   7.206 0.00175 **
## Residuals   51   7199   141.1                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

或者，


```r
tens.lm <- lm(breaks ~ tension, data = warpbreaks)
anova (tens.lm)
```

```
## Analysis of Variance Table
## 
## Response: breaks
##           Df Sum Sq Mean Sq F value   Pr(>F)   
## tension    2 2034.3 1017.13  7.2061 0.001753 **
## Residuals 51 7198.6  141.15                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

和常规线性模型一样，检验因子之间的交互作用也很简单：


```r
summary(aov(breaks~wool*tension, warpbreaks))
```

```
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## wool          1    451   450.7   3.765 0.058213 .  
## tension       2   2034  1017.1   8.498 0.000693 ***
## wool:tension  2   1003   501.4   4.189 0.021044 *  
## Residuals    48   5745   119.7                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

可视化此类交互作用的最简单方法是使用一个虽不美观但实用的交互作用图：


```r
# with(warpbreaks,interaction.plot(tension,wool,breaks))
```

（它需要 3 个参数：第一个是 x 轴的因子，然后是不同线条的追踪因子，最后是响应变量）。当然，你也可以用条形图（比如用 ggplot 让它更美观）。

在方差分析中，你还可以选择协方差分析（ANCOVA）—— 使用协变量来吸收方差，或者检验变量 X1 是否在 X2 的基础上还有额外效应，就像在回归中纳入控制变量 / 固定效应一样，所以 R  代码非常相似。如果要将 tension 作为协变量（在这个数据集里没什么意义，不过只是举个例子）：


```r
summary(aov(breaks ~ tension + wool, warpbreaks))
```

```
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## tension      2   2034  1017.1   7.537 0.00138 **
## wool         1    451   450.7   3.339 0.07361 . 
## Residuals   50   6748   135.0                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

关键是协变量必须在代码中列在前面，因为显然 R 是逐步读取方差分析项的。

另外两个主要内容是线性对比和重复测量分析…… 遗憾的是，我对这些主题记得不多，也没有现成的代码  :-/ 但怡丹和敏都很擅长统计和 R，她们可能能帮上忙！

哦，另外一种（更简单的）检验组间真实差异所在的方法（如果你发现因子整体有显著效应）是使用 Tukey HSD 检验：


```r
TukeyHSD(tens.aov)
```

```
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = breaks ~ tension, data = warpbreaks)
## 
## $tension
##           diff       lwr        upr     p adj
## M-L -10.000000 -19.55982 -0.4401756 0.0384598
## H-L -14.722222 -24.28205 -5.1623978 0.0014315
## H-M  -4.722222 -14.28205  4.8376022 0.4630831
```

这会输出因子中所有组间的 pairwise 比较结果。