给女朋友写的生统资料_Part14

作者: 城管大队哈队长 | 来源:发表于2019-06-10 19:26 被阅读0次

之前我们提到了如果做多次的假设检验，就要考虑多重比较矫正的问题了。那有没有只用做一次检验就可以搞定的方法呢。其实是有的，就是方差分析（ANOVA）了。不过由于ANOVA假设的是大家都是一样的，所以统计假设的意义不是那么的明确，还需要做事后的检验。

单因素ANOVA

首先把方差分析的假设性检验放出来（考试应该要写这个）：
$H0:\mu_1=\mu_2=……=\mu_k \\ H1：\mu_1,\mu_2……\mu_k不全相等$
同时设定 $\alpha$ 值为0.05.

别忘了设 $\alpha$ 值为0.05

但其实H0和H1还有另一种写法，就是因素A对结果没有影响：
$H0:\alpha_1=\alpha_2……=\alpha_k=0 \\ H1：\alpha_1，\alpha_2……不全为0$
这里面 $\alpha_k$ 代表的是处理的各个分组或者说水平。

ANOVA用的是 F 检验，即比较两种方差的大小。一种是组间方差，一种是组内方差。为什么这样就能检验各组均值是否是一样的呢。想象一下，我们总的变异应该来源于两部分，一部分是因素（比如某些处理，某些药物）导致的变异（组间变异），另一部分是随机抽样导致的抽样误差（组内变异）。如果大家都是来源于一个总体的，那么因素造成的变异就应该跟抽样误差的变异是差不多的。具体的推导建议大家去看书。

我这里把F检验的步骤列一下（考试应该也要写这个）。

我们假设总共有 k 组，每组有 m 个重复，总共为 n 个数。首先我们要把总平方和分解成组内平方和和组间平方和
$SS_{total} = SS_{between}+SS_{within}$
总平方和为
$SS_{total} = \sum(x_{ij}-\bar{\bar{X}})^2$
组内平方和为
$SS_{between}=\sum_j \sum_i (\bar{X_j}-\bar{\bar{X}})^2$

当然组间平方和比较让人看得懂的写法应该写成下面的式子，不过老师PPT上是按上面这么写的，反正也不会让你算。。。。。：

$SS_{between}=\sum_{j=1}^k m(\bar{X_j}-\bar{\bar{X}})^2$

组间平方和为
$SS_{within}=\sum_j \sum_i(x_{ij}-{\bar{X_j}})^2$

看得懂的版本

$SS_{within}=\sum_{j=1}^k\sum_{i=1}^{m}(x_{ij}-{\bar{X_j}})^2$

在计算完了组间和组内平方和之后，就可以计算 F 统计量了。
$F =\frac{MS_{between}}{MS_{within}}=\frac{SS_{between}/(k-1)}{SS_{within}/(n-k)}$
k-1 我们称为组间自由度 $df_{between}$ ，n-k 我们称为组内自由度 $df_{within}$ 。

当然，R里面只需要用到 aov 这个函数就可以了。不过在用 aov 这个函数之前，即做ANOVA分析之前。还需要做正态性检验，以及方差齐性检验。（就像之前讲过的 t-test 一样）。

关于这个流程，用一张PPT上的图表示：

14_1.png

首先检验正态性，用到的检验函数是 shapiro.test
即如果每组样本不是正态性的样本，那么就只能做非参数的ANOVA了。非参数用的函数为：kruskal.test
如果每组样本都是正态性的样本，那么接下来就要检验方差齐性。检验的函数是用的是 bartlett.test
如果方差不是齐性的，那么就要用 oneway.test()。在 oneway.test 那里设置 var.equal = FALSE（不过FASLE是默认的。。。。。。）
如果方差是齐性的，才可以最后用 aov 函数。

举个例子，用的是我们最一开始练习R操作时候的 test1数据集

# 读入数据，并把 seed 那列变成因子
test1 <- read.table("rawdata/test1.txt",header = T)
test1$seed <- factor(test1$seed)

# 正态性检验，发现都是正态的
# 假设检验为  H0：种子1的产量满足正态分布；H1：种子1的产量不满足正态分布。
# 种子2,3假设同理

shapiro.test(subset(test1,seed == 1)$yield)
shapiro.test(subset(test1,seed == 2)$yield)
shapiro.test(subset(test1,seed == 3)$yield)

# 方差齐性检验，发现方差其实是齐性
# H0：三组数据方差没有差别；H1：方差有差别

bartlett.test(yield ~ seed, data = test1)

# 都满足之后就可以用 aov 了，用summary检查结果
> test1_aov <- aov(yield ~ seed, data = test1)
> summary(test1_aov)
            Df Sum Sq Mean Sq F value   Pr(>F)    
seed         2   4364  2182.2   10.04 0.000586 ***
Residuals   26   5649   217.3                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

可以看到aov输出的结果，把 $SS_{between}$ 、 $SS_{within}$ 、F等表示出来了。

我们之前提到过，方差分析并不能知道哪组跟哪组有差异。但如果我们想知道的话，就要用事后检验了。R 里面的函数是 TukeyHSD()

> TukeyHSD(test1_aov)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = yield ~ seed, data = test1)

$seed
         diff        lwr       upr     p adj
2-1 -25.32727 -41.330375 -9.324171 0.0015630
3-1  -0.10000 -17.473293 17.273293 0.9998872
3-2  25.22727   8.208575 42.245971 0.0029482

稍微要注意的一点是，aov 函数要求输入的数据都是我们之前提到过的长数据，如果碰到了宽数据，记得要转换一下。

多因素ANOVA

设有两个因素A和B，两因素的ANOVA的总平方和分解为
$SS_{total}=SS_{A}+SS_{B}+SS_{AB}+SS_E$
即总平方和可以分解为A因素造成的变异，B因素造成的变异，AB互作造成的变异以及抽样误差。具体的分解计算看PPT算了，太过于复杂，我觉得考试应该不用写。

两因素方差分析的假设为：

零假设有三个，分别是因素A（有l个水平）没有效果，因素B（有m个水平）没有效果，A和B互作因素没有效果
$H0_A：\alpha_1=\alpha_2=……=\alpha_l=0 \\ H0_B: \beta_1=\beta_2=……=\alpha_m=0 \\ H0_{AB}: ab_{ij}=0 \\ i=1,2……l;j=1,2……,m$
备则假设就是不全为0。

两因素ANOVA的 F统计量有三个。见PPT图

14_2.png

14_3.png

再次举个例子，用的是我们之前讲过的 test4 的数据。A和B为两种因素，A有三个水平，B有2个水平。

关于 aov 函数的用法，可以去看看《R语言实战》第二版的201页。这里放出我感觉比较有用的一张截图。

14_4.png

# 读取和准备数据（这部分先前已经讲过了）

test4 <- read.table("rawdata/test4.txt",header = T)

test4_long <- gather(test4, key = temperature, value = weight, A1, A2, A3)
test4_long$temperature <- factor(test4_long$temperature)
feed <- c(rep("B1",10),rep("B2",10))
test4_data <- cbind(feed,test4_long)

> head(test4_data)
  feed temperature weight
1   B1          A1  282.1
2   B1          A1  264.2
3   B1          A1  274.2
4   B1          A1  276.4
5   B1          A1  283.7
6   B1          A1  288.0

# 正态性检验
# 检验体重数据对于因素A和因素B是否是正态的
shapiro.test(subset(test4_data, feed == "B1")$weight)
shapiro.test(subset(test4_data, feed == "B2")$weight)
shapiro.test(subset(test4_data, temperature == "A1")$weight)
shapiro.test(subset(test4_data, temperature == "A2")$weight)
shapiro.test(subset(test4_data, temperature == "A3")$weight)

# 方差齐性检验
bartlett.test(weight ~ feed, test4_data)
bartlett.test(weight ~ temperature, test4_data)

# ANOVA
> test4_aov <- aov(weight ~ feed * temperature, data = test4_data)
> summary(test4_aov)
                 Df Sum Sq Mean Sq F value   Pr(>F)    
feed              1    127     127   1.589    0.213    
temperature       2   9080    4540  56.809 5.22e-14 ***
feed:temperature  2     17       9   0.108    0.897    
Residuals        54   4316      80                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

可以看到有3个p-value。就可以看到温度对作用是显著的。

多因素也可以做事后检验。可以看到，事后检验也是分因素A，因素B，因素AB互作的。

> test4_pos <- TukeyHSD(test4_aov)
> test4_pos
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = weight ~ feed * temperature, data = test4_data)

$feed
      diff       lwr      upr     p adj
B2-B1 2.91 -1.717782 7.537782 0.2128406

$temperature
       diff       lwr       upr     p adj
A2-A1 25.39 18.576905 32.203095 0.0000000
A3-A1 26.75 19.936905 33.563095 0.0000000
A3-A2  1.36 -5.453095  8.173095 0.8805321

$`feed:temperature`
             diff        lwr       upr     p adj
B2:A1-B1:A1  1.87  -9.942078 13.682078 0.9970585
B1:A2-B1:A1 25.09  13.277922 36.902078 0.0000009
B2:A2-B1:A1 27.56  15.747922 39.372078 0.0000001
B1:A3-B1:A1 25.49  13.677922 37.302078 0.0000006
B2:A3-B1:A1 29.88  18.067922 41.692078 0.0000000
B1:A2-B2:A1 23.22  11.407922 35.032078 0.0000050
B2:A2-B2:A1 25.69  13.877922 37.502078 0.0000005
B1:A3-B2:A1 23.62  11.807922 35.432078 0.0000035
B2:A3-B2:A1 28.01  16.197922 39.822078 0.0000001
B2:A2-B1:A2  2.47  -9.342078 14.282078 0.9892562
B1:A3-B1:A2  0.40 -11.412078 12.212078 0.9999985
B2:A3-B1:A2  4.79  -7.022078 16.602078 0.8358408
B1:A3-B2:A2 -2.07 -13.882078  9.742078 0.9952518
B2:A3-B2:A2  2.32  -9.492078 14.132078 0.9919360
B2:A3-B1:A3  4.39  -7.422078 16.202078 0.8800277

缺失数据

对于缺失值，我感觉我们考试如果遇到的话，应该只要考虑直接去掉就可以了。举个例子

# 数据模拟
> dat <- data.frame(V1 = 1:3,V2 = c(2,3,NA))
> dat
  V1 V2
1  1  2
2  2  3
3  3 NA

# 检查缺失值
> is.na(dat)
        V1    V2
[1,] FALSE FALSE
[2,] FALSE FALSE
[3,] FALSE  TRUE
> table(is.na(dat))

FALSE  TRUE 
    5     1 
    
    
# 凡是有缺失值的那一行就去掉
> na.omit(dat)
  V1 V2
1  1  2
2  2  3

给女朋友写的生统资料_Part14

单因素ANOVA

多因素ANOVA

缺失数据

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

统计