单因素方差分析

目录

前言

1. 为什么不能两两比较？
2. 方差分析（ANOVA）原理
3. 实例讲解
   3.1 提出问题
   3.2 画图观察
   3.3 计算各误差平方和
   3.4 计算F检验值
   3.5 R语言代码
4. 判定系数
5. 事后检验
6. 参考资料
   后记

前言

我们知道，在比较两个分组之间有没有差异时，我们会首选T test进行分析。如果样本量太小或者数据分布不满足正态性时，我们会选择[Wilcoxon检验]Wilcoxon检验 - 简书 (jianshu.com)。

但是，在我们课题中，我们的实验组可能不止2组，例如：用A药组+用B药组+用C药组+用D药组+……

在这种情况下，我们该怎么办呢？

1. 为什么不能两两比较？

最简单来说，我们可能会想着把所有的组分成两两比较的组。例如：对于A组+B组+C组+D组来说，我们可能会用T test去分别比较：

A+B

A+C

A+D

B+C

B+D

C+D
但是这里会存在一些问题：

检验过程繁琐：比较次数多，这里因为只有4组，所以一共需要比较C42=6次T test进行检验。当分组更多时，检验次数则呈现指数化的增长。

检验的灵敏度降低：每次比较只用了部分数据，没有考虑整体数据变化，使得计算结果的准确性降低，从而降低了检验的灵敏度。

结果可靠性低：在我们认同的a=0.05的前提下，也就是每次检验正确的可能性是0.95。那么6次检验同时正常的概率=(0.95)6=0.735。那么这个时候检验的a=1-0.735=0.265，远远大于我们普遍接受的0.05

所以我们需要一种检验来同时比较多组均值之间的关系，这就引出了我们的方差分析检验

2. 方差分析（ANOVA）原理

1.方差分析的简介

它的基本思想是将测量数据的总变异（即总方差）按照变异来源分为处理（组间）效应和误差（组内）效应，并作出其数量估计，从而确定实验处理对研究结果影响力的大小。

按因素划分，可分为单因素方差分析、二因素方差分析和多因素方差分析。

方差分析的步骤为：总平方和分解→总自由度分解→F检验。若F检验显著，则可以进行多重比较，从而发现哪些处理具有两两间的差异。

2.方差分析的基本原理

在一次实验中，可以得到一系列不同的观测值。造成观测值不同的原因可能是①由于处理因素不同引起的，即处理效应；也可能是②由于实验过程中偶然性因素的干扰和测量误差所致，即误差效应。
反应测量数据变异性的指标有多个，在方差分析中选用方差来度量资料的变异程度。要正确认识观测值的变异是由于处理效应还是误差效应引起的，我们可以分别计算出处理效应的方差以及误差效应的方差，在一定显著水平下进行比较，如果二者相差不大，说明实验处理对观测值的影响不大；如果差异较大则说明实验处理对于观测的影响较大。

• 全体样本的总误差平方和（总变异），也成为SST：[图片上传失败

• 每个组内的误差平方和（组内变异），也称为SSE：[图片上传失败

• 每个组间的误差平方和（组间变异），也称为SSA：[图片上传失败

• 他们的之间的数学关系：SST=SSE+SSA

i：i组，例子中的i可以理解为a,b,c,d

nj：第j组的总人数

j：每个组中的单个样本

：全体样本数据的总均值

：每组样本的均值

：第 i 组中的第 j 个样本

image.png

我们把变异分为组内变异和组间变异两部分：

image.png

仔细观察下面的例子，左图各组样本均数不相同，而右图则样本均数一致。

那我们再看下其两类组间变异：左图中组间差异很大，右图中组间差异很小。

我们再看下其两类组内变异：，左图中组内差异较小，而右图中组内差异较大。

image.png

如果组间差异（B）远大于组内差异（A），就意味着各组样本均数不一致呢？

方差分析就是基于这样的思路：

以组内差异（A）为参考基准，考察组间差异（B）的大小。

如果组间差异（B）远大于组内差异（A），则认为组间存在区别。

而组内差异，我们认为是因为（完全）随机而产生的。以这样一个完全随机的尺度作为标杆，也甚是巧妙。

我们也引出了F值，即B比A的值。

image.png

基于F分布，就很容易看出，当组间差异越大，那么F=B/A值也越大，横坐标越向右，越容易进入我们拒绝原假设（H0，各组均数相同）的区域。

也就是说，当组间差异越大，p值越小，越有理由拒绝原假设H0，也就是两组之间的差异越有显著性。

3.实例讲解

3.1 提出问题

假设我们收集了不同药物降低体重的数据：

image.png

那么这4种药物之间的效果有没有差异呢？

3.2 画图观察

我们首先画个散点图看看，具体画图代码如下：

library(ggplot2) 2dat <- data.frame(a=c(57,66,49,40,34,53,44), 3 b=c(68,39,29,45,56,51,NA), 4 c=c(31,49,21,34,40,NA,NA), 5 d=c(44,51,65,77,58,NA,NA)) 6gdat <- melt(dat,na.rm = T) 7colnames(gdat) 8gdat$variable <- as.character(gdat$variable) 9mdat <- data.frame(x=c('a','b','c','d'),10 y=sapply(dat, mean,na.rm=T))11ggplot()+12 geom_point(gdat, mapping = aes(y=value, x=variable))+13 geom_line(mdat,mapping = aes(y=y, x=x),color='red',group=1)

结果如下：

image.png

从图中我们可以大致看出这4组药物的效果之间存在一定差异。

3.3 计算各误差平方和

SST的计算：直接计算即可，计算得到4164.609

SSA的计算：直接计算即可，计算得到1456.609

SSE的计算：直接计算即可，计算得到2708

这些指标的计算代码如下：

1rm(list = ls()) 2options(stringsAsFactors = F) 3library(reshape2) 4dat <- data.frame(a=c(57,66,49,40,34,53,44), 5                  b=c(68,39,29,45,56,51,NA), 6                  c=c(31,49,21,34,40,NA,NA), 7                  d=c(44,51,65,77,58,NA,NA)) 8gdat <- melt(dat,na.rm = T) 9# 定义函数diff sum10ds <- function(n,m){11  sum((n-m)**2,na.rm=T)12}13#SST 14ds(gdat$value,mean(gdat$value))#4164.60915#SSA 16sum(((sapply(dat, mean, na.rm=T)-mean(gdat$value))**2)*table(gdat$variable))#1456.60917#SSE 18sum(c(ds(dat$a,mean(dat$a,na.rm=T)),19      ds(dat$b,mean(dat$b,na.rm=T)),20      ds(dat$c,mean(dat$c,na.rm=T)),21      ds(dat$d,mean(dat$d,na.rm=T))))#2708

这里再次印证了前面的结论：SST=SSE+SSA

3.4 计算F检验值

1.因为误差平方和会受到样本数的影响，所以为了消除样本数的影响，我们将各个误差平方和除以（样本数-1），其实也就是方差

2.综合前面所有数据，得到的结果值统计如下：

image.png

(*) ：此处的方差不可通过MSA+MSE计算而来

3.如果不同组药物之间的效果一致，或者相差不大，那么MSA和MSE之间的比值与1接近。所以我们需要对F-test=MSA/MSE进行检验，得到P值。

4.F-test值服从分子自由度为k-1、分母自由度为n-k的F分布，即：MSA/MSE=F-test~F(i-1,n-i)

5.查表得知P<0.05

3.5 R语言代码

rm(list = ls()) 2options(stringsAsFactors = F) 3library(reshape2) 4dat <- data.frame(a=c(57,66,49,40,34,53,44), 5                  b=c(68,39,29,45,56,51,NA), 6                  c=c(31,49,21,34,40,NA,NA), 7                  d=c(44,51,65,77,58,NA,NA)) 8gdat <- melt(dat,na.rm = T) 9fit <- aov(value ~ variable,data = gdat)10summary(fit)11#               Df Sum Sq Mean Sq F value Pr(>F)  12#   variable     3   1457   485.5   3.407 0.0388 *13#   Residuals   19   2708   142.5                 14# ---15#   Signif. codes:  0 '***’ 0.001 '**’ 0.01 '*’ 0.05 '.’ 0.1 ' ’ 1

可以看到，这个结果和我们手动计算的结果是一致的！F=3.407，对应的P=0.0388

4. 判定系数

1.判定系数（Coefficient of Determination，记为R2）在统计学中用于度量因变量的变异中可由自变量解释部分所占的比例，以此来判断统计模型的解释力。

2.决定系数是相关系数的二次幂，例如在这个例子中，我们只有一个自变量——不同的治疗方案。通过计算不同方案的决定系数，我们可以得知不同治疗方案和效果的相关系数R。

3.决定系数意义：拟合优度越大，自变量对因变量的解释程度越高，自变量引起的变动占总变动的百分比高。观察点在回归直线附近越密集。

4.在方差分析中，我们可以通过计算全部变异（SST）中，由分组不同所解释变异（SSA）的占比（SSA/SST）来计算判定系数：

image.png

所以可以计算得知相关系数R：

image.png

上述数据说明了：

R2=0.3498：不同的药物解释了数据34.98%的变异。

R=0.5914：不同的药物选择与疗效之间的相关性为0.5914。

5. 事后检验

1.有了方差分析的结果，我们只能说明j个总体均值不全相等。若想进一步了解哪些两个总体均值不等，需进行多个样本均值间的两两比较或称多重比较（multiple comparison），也称为post hoc（事后）检验。

2.事后检验的方法非常多，这里不展开说明，仅列出常见的方法：

SNK-q检验

LSD-t检验：同过检验两个总体均值是否相等的T检验，其中总体方差用MSE来代替而得到的。

image.png

Dunnett检验

Tukey HSD检验

Duncan检验

Scheffe检验
1.下面我用Tukey HSD检验对我们的结果进行检验，并告诉大家如何解读自己的结果。这里只展示代码，具体原理日后有机会再补充。

1rm(list = ls()) 2options(stringsAsFactors = F) 3library(reshape2) 4dat <- data.frame(a=c(57,66,49,40,34,53,44), 5 b=c(68,39,29,45,56,51,NA), 6 c=c(31,49,21,34,40,NA,NA), 7 d=c(44,51,65,77,58,NA,NA)) 8gdat <- melt(dat,na.rm = T) 9fit <- aov(value ~ variable,data = gdat)10TukeyHSD(fit)11# Tukey multiple comparisons of means12# 95% family-wise confidence level13#14# Fit: aov(formula = value ~ variable, data = gdat)15#16# $variable17# diff lwr upr p adj18# b-a -1 -19.676096 17.676096 0.998739019# c-a -14 -33.656024 5.656024 0.221814420# d-a 10 -9.656024 29.656024 0.496679021# c-b -13 -33.327070 7.327070 0.304637822# d-b 11 -9.327070 31.327070 0.444784823# d-c 24 2.769068 45.230932 0.0234078

可以看到，在事后检验中，仅仅只有d-c之间存在显著差异，而其他检验之间没有显著差异。

我们可以把结果进行绘制图形：

1plot(TukeyHSD(fit))