北河以北

用 R 做方差分析(ANOVA)

文/宋春林

方差分析(Analysis of variance, ANOVA)是统计学中的经典方法,1923年由英国统计学家R. A. Fisher提出,其实质是关于观测值变异原因的数量分析,主要用途是研究外界因素或试验条件对观测结果影响的显著性。方差分析方法就是根据平方和的加和性原理,将总平方和分解为不同影响因素产生的平方和以及代表重复观测值之间的随机变化的误差平方和,并利用假设检验(F-检验)手段分析不同影响因素对总平方和贡献的显著性程度,进而判断这些因素是否对观测结果产生明显影响。

方差分析的前提假设有:采样随机性,总体中每个个体被采集的机会均等;个体独立性,每个样本的采集不受上一个样本的影响;分布正态性,无需多言,数理统计中最重要的假定。对于非正态分布的总体,可以通过数据变换或改用非参数检验方法;方差同质性,要求各个样本(同一因素不同水平的一组重复观测值)的方差大小没有显著差异;方差加和性,即一组观测值的总平方和可以分解成不同影响因素造成的平方和之和。

对于完全随机设计试验且处理数大于2时可以用单因素方差分析(等于2 时用t检验)。单因素方差分析(One-way ANOVA)应用举例:

随机在A、B、C三地抽取家蝇,测量翅膀长度,一共五十个样本,数字特征分别为μ1、μ2和μ3。问题:三地家蝇翅膀长度是否有差异?H0假设:μ1=μ2=μ3,即三地家蝇翅膀长度无显著差异;Ha假设:μ1,μ2,μ3不完全相等,即三地家蝇翅膀长度至少有一个与其他样地有显著差异,α=0.05。R代码和结果:

 
site1 <- c(45, 44, 43, 47, 48, 44, 46, 44, 40, 45, 42, 40, 43, 46, 47, 45, 46, 45, 43, 44) 
site2 <- c(45, 48, 47, 43, 46, 47, 48, 46, 43, 49, 46, 43, 47, 46, 47, 46, 45, 46, 44, 45, 46, 44, 43, 42, 45) 
site3 <- c(47, 48, 45, 46, 46, 44, 45, 48, 49, 50, 49, 48, 47, 44, 45, 46, 45, 43, 44, 45, 46, 43, 42) fly.survey <- data.frame(length = c(site1, site2, site3), site = factor(c(rep("1", 20), rep("2", 25), rep("3", 23)))) 
options(digits = 3) 
tapply(fly.survey$length, fly.survey$site, mean) 
tapply(fly.survey$length, fly.survey$site, var) #
boxplot(length ~ site, data = fly.survey, xlab = "Sites", ylab = "Length") 
bartlett.test(length ~ site, data = fly.survey)
fit <- aov(length ~ site, data = fly.survey) summary(fit) 
plot(fit) 

单因子方差分析结果显示F value = 3.24 ,Pr(>F) = 0.045,因此拒绝H0假设,即认为三地家蝇翅膀长度在统计学上有显著差异。

两因素方差分析(Two-way ANOVA)可用于随机区组实验设计,用来分析两个因素的不同水平对结果是否有显著影响,以及两因素之间是否存在交互效应。应用举例: 有一个关于检验毒品强弱的试验,给48只老鼠注射I, II, III三种毒药(因素A),同时有A, B, C, D 4种治疗方案(因素B),这样的试验在每一种因素组合下都重复四次测试老鼠的存活时间(年)。试分析毒药和治疗方案以及它们的交互作用对老鼠存活时间有无显著影响。

H01假设:α1=α2=α3=0,即因素A毒药种类对老鼠存活时间没有影响;Ha1假设:α1,α2和α3不全相等,即因素A毒药种类对老鼠存活时间有影响。

H02假设:β1=β2=β3=β4=0,即因素B治疗方案对老鼠存活时间没有影响;Ha2假设:β1,β2,β3和β4不全相等,即因素B治疗方案对老鼠存活时间有影响。

H03假设:δ11=δ12=δ13=δ14=δ21=δ22=δ23=δ24=δ31=δ32=δ33=δ34=0,即A因素和B因素没有交互作用;Ha3假设:δ11,δ12,…,δ34不全相等,即A因素和B因素有交互作用。α=0.05。

R代码和结果:

 
rats <- data.frame(time = c(0.39, 0.45, 0.46, 0.43, 0.72, 0.8, 0.78, 0.92, 0.53, 0.45, 0.63, 0.56, 0.65, 0.71, 0.66, 0.62, 0.36, 0.29, 0.23, 0.33, 0.42, 0.61, 0.79, 1.04, 0.94, 0.85, 0.81, 0.6, 0.76, 0.82, 0.71, 0.65, 0.22, 0.31, 0.28, 0.33, 0.65, 0.67, 0.69, 0.78, 0.43, 0.45, 0.54, 0.48, 0.38, 0.36, 0.41, 0.33), toxicant = gl(3, 16, 48, labels = c("I", "II", "III")), cure = gl(4, 4, 48, labels = c("A", "B", "C", "D"))) 
op <- par(mfrow = c(1, 2)) plot(time ~ toxicant + cure, data = rats) 
bartlett.test(time ~ toxicant, data = rats) ## ## Bartlett test of homogeneity of variances ## ## data: time by toxicant 
bartlett.test(time ~ cure, data = rats) 
rats.aov <- aov(time ~ toxicant * cure, data = rats) summary(rats.aov)

由p值可以判断,拒绝H01和H02假设,即因素A毒素和因素B治疗方案对存活时间有显著影响;同时也拒绝H03假设,即认为因素A和因素B交互作用对老鼠存活时间有显著影响。

参考资料:

[1]陶澍编著.应用数理统计方法[M].北京市:中国环境科学出版社.1994.

[2]明道绪主编.生物统计附试验设计[M].北京市:中国农业出版社.2008.

[3]汤银才主编.R语言与统计分析[M].北京市:高等教育出版社.2008.

发表于
分类 学无止境  标签 R  统计学  方差分析