方差分析R语音3
一、基础原理
常见研究设计的表达式
| 设计 | 表达式 |
|---|---|
| 单因素ANOVA | y ~ A |
| 含单个协变量的单因素ANCOVA | y ~ x + A |
| 双因素ANOVA | y ~ A * B |
| 含两个协变量的双因素ANCOVA | y ~ x1 + x2 + A*B |
| 随机化区组 | y ~ B + A（B是区组因子） |
| 单因素组内ANOVA | y ~ A + Error(Subject/A) |
| 含单个组内因子（W）和单个组间因子（B）的重复测量ANOVA | y ~ B * W + Error(Subject/W) |
表达式中效应的顺序在两种情况下会造成影响：
- (a)因子不止一个，并且是非平衡设计；
- (b)存在协变量。出现任意一种情况时，等式右边的变量都与其他每个变量相关。
此时，我们无法清晰地划分它们对因变量的影响。例如，对于双因素方差分析，若不同处理方式中的观测数不同，那么模型 y ~ AB 与模型 y ~ BA 的结果不同。
样本大小越不平衡，效应项的顺序对结果的影响越大。一般来说，越基础性的效应越需要放在表达式前面。具体来讲，首先是协变量，然后是主效应，接着是双因素的交互项，再接着是三因素的交互项，以此类推。对于主效应，越基础性的变量越应放在表达式前面，因此性别要放在处理方式之前。有一个基本的准则：若研究设计不是正交的（也就是说，因子和/或协变量相关），一定要谨慎设置效应的顺序。
二、案例分析
单因素方差分析
案例1（文献1）
案例描述
比较分类因子定义的两个或多个组别中的因变量均值。
以multcomp 包中的 cholesterol 数据集为例（取自Westfall、Tobia、Rom、Hochberg，1999），50个患者均接受降低胆固醇药物治疗（ trt ）五种疗法中的一种疗法。其中三种治疗条件使用药物相同，分别是20mg一天一次（1time）、10mg一天两次（2times）和5mg一天四次（4times）。剩下的两种方式（drugD和drugE）代表候选药物。哪种药物疗法降低胆固醇（响应变量）最多呢？
基础代码
（详细代码见附录D.1）
# 载入数据
#install.packages("multcomp")
library(multcomp)
attach(cholesterol)
#方差分析
fit <- aov(response ~ trt)
summary(fit)
分析结果
Df Sum Sq Mean Sq F value Pr(>F)
trt 4 1351.4 337.8 32.43 9.82e-13 ***
Residuals 45 468.8 10.4
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
结果解读
从输出结果可以看到，每10个患者接受其中一个药物疗法➊。均值显示drugE降低胆固醇最
多，而1time降低胆固醇最少➋，各组的标准差相对恒定，在2.88到3.48间浮动➌。ANOVA对治
疗方式（ trt ）的F检验非常显著（p<0.0001），说明五种疗法的效果不同➍。
虽然ANOVA对各疗法的F检验表明五种药物疗法效果不同，但是并没有告诉你哪种疗法与其
他疗法不同。多重比较可以解决这个问题。例如， TukeyHSD() 函数提供了对各组均值差异的成
对检验
单因素协方差分析
案例1（文献1）
案例描述
单因素协方差分析（ANCOVA）扩展了单因素方差分析（ANOVA），包含一个或多个定量的
协变量。下面的例子来自于 multcomp 包中的 litter 数据集（见Westfall et al.，1999）。
怀孕小鼠被分为四个小组，每个小组接受不同剂量（0、5、50或500）的药物处理。产下幼崽的体重均值
为因变量，怀孕时间为协变量。
基础代码
（详细代码见附录D.2）
data(litter, package="multcomp")
attach(litter)
fit <- aov(weight ~ gesttime + dose)
summary(fit)
分析结果
Df Sum Sq Mean Sq F value Pr(>F)
gesttime 1 134.3 134.30 8.049 0.00597 **
dose 3 137.1 45.71 2.739 0.04988 *
Residuals 69 1151.3 16.69
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
结果解读
双因素方差分析
案例1（文献1）
案例描述
在双因素方差分析中，受试者被分配到两因子的交叉类别组中。以基础安装中的ToothGrowth 数据集为例，随机分配60只豚鼠，分别采用两种喂食方法（橙汁或维生素C），各喂食方法中抗坏血酸含量有三种水平（0.5mg、1mg或2mg），每种处理方式组合都被分配10只豚鼠。牙齿长度为因变量。
基础代码
（详细代码见附录D.3）
attach(ToothGrowth)
dose <- factor(dose)
fit <- aov(len ~ supp*dose)
summary(fit)
分析结果
Df Sum Sq Mean Sq F value Pr(>F)
supp 1 205.4 205.4 15.572 0.000231 ***
dose 2 2426.4 1213.2 92.000 < 2e-16 ***
supp:dose 2 108.3 54.2 4.107 0.021860 *
Residuals 54 712.1 13.2
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
结果解读
dose 变量被转换为因子变量，这样 aov() 函数就会将它当做一个分组变量，而不是一个数值型协变量。用 summary() 函数得到方差分析表，可以看到主效应（supp和dose）和交互效应都非常显著。
重复测量方差分析
重复测量方差分析，即受试者被测量不止一次。
案例1（文献1）
案例描述
示例来源于生理生态学领域，研究方向是生命系统的生理和生化过程如何响应环境因素的变异（此为应对全球变暖的一个非常重要的研究领域）。基础安装包中的 CO2 数据集包含了北方和南方牧草类植物Echinochloa crus-galli （Potvin、Lechowicz、Tardif，1990）的寒冷容忍度研究结果，在某浓度二氧化碳的环境中，对寒带植物与
非寒带植物的光合作用率进行了比较。
研究所用植物一半来自于加拿大的魁北克省（Quebec），另一半来自美国的密西西比州（Mississippi）。
首先，我们关注寒带植物。因变量是二氧化碳吸收量（ uptake ），单位为ml/L，自变量是植物类型 Type （魁北克VS.密西西比）和七种水平（95~1000 umol/m^2 sec）的二氧化碳浓度（ conc ）。另外， Type 是组间因子， conc 是组内因子。 Type 已经被存储为一个因子变量，但你还需要先将conc 转换为因子变量。
基础代码
（详细代码见附录D.4）
CO2$conc <- factor(CO2$conc)
w1b1 <- subset(CO2, Treatment=='chilled')
fit <- aov(uptake ~ conc*Type + Error(Plant/(conc)), w1b1)
summary(fit)
分析结果
Error: Plant
Df Sum Sq Mean Sq F value Pr(>F)
Type 1 2667.2 2667.2 60.41 0.00148 **
Residuals 4 176.6 44.1
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Error: Plant:conc
Df Sum Sq Mean Sq F value Pr(>F)
conc 6 1472.4 245.40 52.52 1.26e-12 ***
conc:Type 6 428.8 71.47 15.30 3.75e-07 ***
Residuals 24 112.1 4.67
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
结果解读
方差分析表表明在0.01的水平下，主效应类型和浓度以及交叉效应类型×浓度都非常显著，
通过 interaction.plot() 函数展示了交互效应。
多元方差分析
当因变量（结果变量）不止一个时，可用多元方差分析（MANOVA）对它们同时进行分析
案例描述.
以 MASS 包中的 UScereal 数据集为例（Venables，Ripley（1999）），我们将研究美国谷物中的卡路
里、脂肪和糖含量是否会因为储存架位置的不同而发生变化；其中1代表底层货架，2代表中层货
架，3代表顶层货架。卡路里、脂肪和糖含量是因变量，货架是三水平（1、2、3）的自变量。
基础代码
详细代码见D.5
library(MASS)
attach(UScereal)
shelf <- factor(shelf)
y <- cbind(calories, fat, sugars)
fit <- manova(y ~ shelf)
summary(fit)
分析结果
Df Pillai approx F num Df den Df Pr(>F)
shelf 2 0.4021 5.1167 6 122 0.0001015 ***
Residuals 62
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
结果解读
由于多元检验是显著的，可以使用 summary.aov() 函数对每一个变量做单因素方差分析➊。
从上述结果可以看到，三组中每种营养成分的测量值都是不同的。另外，还可以用均值比较步骤
（比如 TukeyHSD ）来判断对于每个因变量，哪种货架与其他货架都是不同的
案例描述
基础代码
结果解读
分析结果
三、案例扩展
附录
附录A 参考文献
A.1 Robert ，R语言实战[M] ,2016
附录D
D.1 单因素方差分析
案例1
# 加载 multcomp 包，该包提供了用于多重比较的函数和方法
library(multcomp)
# 将 cholesterol 数据集添加到搜索路径，这样可以直接使用数据集中的变量名
attach(cholesterol)
# 使用 table() 函数创建 trt 变量的频数表，展示 trt 各个水平的出现次数
table(trt)
# aggregate() 函数用于对数据进行分组计算。
# 这里是按照 trt 变量对 response 变量进行分组，然后计算每组的均值
aggregate(response, by = list(trt), FUN = mean)
# 同样按照 trt 变量对 response 变量进行分组，计算每组的标准差
aggregate(response, by = list(trt), FUN = sd)
# 使用 aov() 函数进行单因素方差分析，研究 trt 对 response 的影响
fit <- aov(response ~ trt)
# 查看方差分析的结果摘要，包括自由度、平方和、均方、F 值和 p 值等信息
summary(fit)
# 安装 gplots 包（如果未安装），gplots 包提供了一些绘图函数
# install.packages("gplots")
# 加载 gplots 包
library(gplots)
# 使用 plotmeans() 函数绘制带有 95% 置信区间的均值图
# xlab 为 x 轴标签，ylab 为 y 轴标签，main 为图形的标题
plotmeans(response ~ trt, xlab = "Treatment", ylab = "Response",
main = "Mean Plot\nwith 95% CI")
# 使用 TukeyHSD() 函数进行 Tukey 多重比较，分析不同处理组之间的差异
TukeyHSD(fit)
# 设置图形中坐标轴标签的方向为垂直，方便查看较长的标签
par(las = 2)
# 设置图形的边距，mar 参数依次为下、左、上、右的边距大小
par(mar = c(5, 8, 4, 2))
# 绘制 Tukey 多重比较的结果图，直观展示不同处理组之间的差异
plot(TukeyHSD(fit))
# 再次加载 multcomp 包（可能前面代码运行后会话状态有变化）
library(multcomp)
# 重新设置图形的边距
par(mar = c(5, 4, 6, 2))
# 使用 glht() 函数进行广义线性假设检验，这里使用 Tukey 方法对 trt 变量进行多重比较
tuk <- glht(fit, linfct = mcp(trt = "Tukey"))
# 绘制使用 compact letter display (CLD) 方法表示的多重比较结果图
# col 参数设置图形的颜色为浅灰色
plot(cld(tuk, level = .05), col = "lightgrey")
# 正态性检验
# 加载 car 包，该包提供了一些用于回归分析和诊断的工具
library(car)
# 使用 qqPlot() 函数绘制 Q-Q 图，用于检验响应变量在不同处理组下是否服从正态分布
# simulate = TRUE 表示进行模拟以得到置信区间，main 为图形标题
qqPlot(lm(response ~ trt, data = cholesterol),
simulate = TRUE, main = "Q-Q Plot", labels = row.names(cholesterol))
# 方差齐性检验
# 使用 bartlett.test() 函数进行 Bartlett 检验，检验不同处理组的方差是否齐性
bartlett.test(response ~ trt, data = cholesterol)
# 离群点检测
# 再次加载 car 包（确保其在当前会话中可用）
library(car)
# 使用 outlierTest() 函数检测方差分析模型中的离群点
outlierTest(fit)
# 显示 cholesterol 数据集的前五行，方便查看数据的基本结构和内容
head(cholesterol, 5)
# 保存 cholesterol 数据集到指定路径
# write.csv() 函数将数据框保存为 CSV 文件
# row.names = FALSE 表示不保存行名到文件中
write.csv(cholesterol,
file = "C:\\Users\\Administrator\\Desktop\\code\\cholesterol_data.csv",
row.names = FALSE)
D.2 单因素协方差分析
案例1
# 清除当前工作空间中的所有变量，避免已有变量对后续操作产生干扰
rm(list = ls())
# 检查工作空间是否为空，输出为空则表明工作空间已清空
ls()
# 从 multcomp 包中加载 litter 数据集
data(litter, package="multcomp")
# 将 litter 数据集添加到搜索路径，这样可以直接使用数据集中的变量名，无需使用 litter$变量名 的形式
attach(litter)
# 创建 dose 变量的频数表，展示 dose 各个水平的出现次数
table(dose)
# 显示 litter 数据集的前五行，用于快速查看数据的基本结构和内容
head(litter,5)
# 使用 aggregate 函数对 weight 变量按照 dose 变量进行分组，并计算每组的均值
aggregate(weight, by=list(dose), FUN=mean)
# 构建一个方差分析模型，研究 gesttime 和 dose 对 weight 的影响
# 这里是一个双因素的方差分析模型，gesttime 和 dose 作为自变量，weight 作为因变量
fit <- aov(weight ~ gesttime + dose)
# 查看方差分析模型 fit 的结果摘要，包括自由度、平方和、均方、F 值和 p 值等信息
summary(fit)
# 安装 effects 包（如果未安装），该包用于可视化和分析线性模型的效应
# install.packages("effects")
# 加载 effects 包
library(effects)
# 计算并输出模型 fit 中 dose 变量的边际效应，用于评估 dose 对 weight 的影响
effect("dose", fit)
# 加载 multcomp 包，该包提供了用于多重比较和假设检验的功能
library(multcomp)
# 定义一个对比矩阵，这里定义了一个名为 "no drug vs. drug" 的对比
# 对比矩阵表示了不同组之间的比较关系，这里是将第一个水平与其他三个水平进行对比
contrast <- rbind("no drug vs. drug" = c(3, -1, -1, -1))
# 使用 glht 函数对模型 fit 进行广义线性假设检验
# linfct 参数指定了线性假设，这里使用 mcp 函数将对比矩阵应用到 dose 变量上
# 最后输出检验结果的摘要
summary(glht(fit, linfct=mcp(dose=contrast)))
# 再次加载 multcomp 包（可能前面代码运行后会话状态有变化）
library(multcomp)
# 构建一个新的方差分析模型 fit2，考虑 gesttime 和 dose 的交互作用
# 使用 * 表示同时包含主效应和交互效应
fit2 <- aov(weight ~ gesttime*dose, data=litter)
# 查看新方差分析模型 fit2 的结果摘要
summary(fit2)
# 安装 HH 包（如果未安装），该包提供了一些高级的统计分析和可视化功能
# install.packages("HH")
# 加载 HH 包
library(HH)
# 进行协方差分析（ANCOVA），研究在控制 gesttime 变量的情况下，dose 对 weight 的影响
# 这里使用 ancova 函数进行分析，并指定数据来源为 litter 数据集
ancova(weight ~ gesttime + dose, data=litter)
D.3 双因素方差分析
案例1
# 加载 ToothGrowth 数据集，并将其添加到搜索路径中，
# 这样可以直接使用数据集中的变量而无需使用 $ 符号来引用
attach(ToothGrowth)
# 创建一个列联表，展示变量 supp（补充剂类型）和 dose（剂量）的各种组合的频数
table(supp, dose)
# 查看 ToothGrowth 数据集的前 5 行数据，以便快速了解数据的结构和内容
head(ToothGrowth, 5)
# 按照 supp（补充剂类型）和 dose（剂量）进行分组，
# 计算每组中 len（牙齿长度）的平均值
aggregate(len, by = list(supp, dose), FUN = mean)
# 按照 supp（补充剂类型）和 dose（剂量）进行分组，
# 计算每组中 len（牙齿长度）的标准差
aggregate(len, by = list(supp, dose), FUN = sd)
# 将 dose 变量转换为因子类型，以便在后续的方差分析中正确处理
dose <- factor(dose)
# 进行双因素方差分析（ANOVA），模型为 len（牙齿长度） ~ supp（补充剂类型）* dose（剂量），
# 这里的 * 表示包含主效应和交互效应
fit <- aov(len ~ supp * dose)
# 输出方差分析模型的摘要信息，包括各因素的 F 值、P 值等，用于评估模型的显著性
summary(fit)
# 绘制交互作用图，展示剂量（dose）和补充剂类型（supp）对牙齿长度（len）的交互影响
# type = "b" 表示绘制点和线的组合图
# col = c("red", "blue") 设置不同补充剂类型对应的线条和点的颜色
# pch = c(16, 18) 设置不同补充剂类型对应的点的形状
# main = "Interaction between Dose and Supplement Type" 设置图形的主标题
interaction.plot(dose, supp, len, type = "b",
col = c("red", "blue"), pch = c(16, 18),
main = "Interaction between Dose and Supplement Type")
D.4 重复测量方差分析
案例1
# 将 CO2 数据集中的 conc 变量转换为因子类型
# 这样做是为了在后续的方差分析和绘图中，将浓度作为分类变量处理
CO2$conc <- factor(CO2$conc)
# 从 CO2 数据集中筛选出 Treatment 变量值为 'chilled' 的数据子集
# 存储在 w1b1 中，后续分析将基于这个子集进行
w1b1 <- subset(CO2, Treatment == 'chilled')
# 进行嵌套设计的双因素方差分析
# 因变量是 uptake（二氧化碳吸收量）
# 自变量包括 conc（浓度）和 Type（植物类型），并考虑它们之间的交互作用
# Error(Plant/(conc)) 表示嵌套结构，即每个 Plant 下有不同的 conc 水平
# 整体意思是在 chilled 处理下，分析浓度和植物类型对二氧化碳吸收量的影响
fit <- aov(uptake ~ conc*Type + Error(Plant/(conc)), w1b1)
# 输出方差分析模型的摘要信息
# 包括各因素的 F 值、P 值等，用于评估模型中各个效应的显著性
summary(fit)
# 设置图形中坐标轴标签的方向
# las = 2 表示坐标轴标签垂直于坐标轴，这样可以避免标签过长时相互重叠
par(las = 2)
# 设置图形的边距
# mar = c(10, 4, 4, 2) 分别表示图形下、左、上、右边距的大小
# 这里将下边距设置为 10，以确保 x 轴标签有足够的空间显示
par(mar = c(10, 4, 4, 2))
# 使用 with 函数，在 w1b1 数据子集的环境中执行绘图操作
# 绘制交互作用图，展示 conc（浓度）和 Type（植物类型）对 uptake（二氧化碳吸收量）的交互影响
# type = "b" 表示绘制点和线的组合图
# col = c("red", "blue") 设置不同植物类型对应的线条和点的颜色
# pch = c(16, 18) 设置不同植物类型对应的点的形状
# main = "Interaction Plot for Plant Type and Concentration" 设置图形的主标题
with(w1b1, interaction.plot(conc, Type, uptake,
type = "b", col = c("red", "blue"), pch = c(16, 18),
main = "Interaction Plot for Plant Type and Concentration"))
# 绘制箱线图，展示不同 Type（植物类型）和 conc（浓度）组合下 uptake（二氧化碳吸收量）的分布情况
# data = w1b1 表示使用 w1b1 数据子集
# col = (c("gold", "green")) 设置不同箱线的颜色
# main = "Chilled Quebec and Mississippi Plants" 设置图形的主标题
# ylab = "Carbon dioxide uptake rate (umol/m^2 sec)" 设置 y 轴的标签
boxplot(uptake ~ Type*conc, data = w1b1, col = (c("gold", "green")),
main = "Chilled Quebec and Mississippi Plants",
ylab = "Carbon dioxide uptake rate (umol/m^2 sec)")
多元方差分析
案例1
library(MASS)
attach(UScereal)
shelf <- factor(shelf)
y <- cbind(calories, fat, sugars)
aggregate(y, by=list(shelf), FUN=mean)
cov(y)
fit <- manova(y ~ shelf)
summary(fit)
浙公网安备 33010602011771号