R语言 方差分析
Posted 基督徒Isaac
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了R语言 方差分析相关的知识，希望对你有一定的参考价值。
函数总结：
处理数据
pivot_longer()
mutate(factor)
arrange()
count()
summarise()
levels()查看因子变量各水平，table()查看频数，hist(as.numeric(unlist))绘制直方图
建模判断
aov() %>% summary() # 一元因变量（任意自变量、协变量个数）、多元因变量
manova() %>% summary() # 多元方差分析
manova() %>% summary.aov() # 对多元因变量分别做单因素分析
effects::effect() # 去除协变量效应后的组均值
anova()
car::Anova()
lm() # 线性模型要求预测变量是数值型，lm()用一系列与因子水平相对应的数值型对照变量来代替因子
ancova(~) # 不需要假设回归斜率同质性，允许斜率和截距项依据组别而发生变化
sm::sm.ancova(~) # 不需要假设回归斜率同质性的非参数ANCOVA方法
rrcov::Wilks.test(method=“mcd”) # 稳健单因素MANOVA
vegan::adonis() # 非参数MANOVA
初级可视化（均值、方差、分布）
group_by() %>% plot()
boxplot()
ggplot + geom_boxplot
gplots::plotmeans()
高级可视化
gplots::plotmeans(interaction) # 交互效应
interaction.plot() # 交互效应
boxplot() # 交互效应的不同侧面
HH::ancovaplot() # 含协变量
HH::interaction2wt() # 含多因素、交互效应
多重比较
TukeyHSD()
multcomp::glht(linfct = mcp(“Tukey”))
multcomp::glht(linfct = mcp(rbind))
假设检验
正态性
car::qqPlot(lm)
car::qqPlot(
qchisq(ppoints(nrow), df = ncol),
mahalanobis(cbind, colMeans, cov)) %>% # 检验多元正态性
identify(labels = row.names) # 识别异常点，点击返回索引值
方差齐性
bartlett.test
fligner.test
HH::hov
Box’s M Test # 检验方差-协方差矩阵同质性，即各组的协方差矩阵是否相同。该检验对正态性假设很敏感。
离群点
car::outlierTest
mvoutlier::ap.plot() # 多元离群点
aov(~*) # 回归斜率的同质性
参考文献：
https://www.math.pku.edu.cn/teachers/lidf/docs/Rbook/html/_Rbook/stat-aov.html
https://bbs.pinggu.org/thread-5556558-1-1.html
R语言实战（第二版）
案例：单因素方差分析
（数据见https://www.math.pku.edu.cn/teachers/lidf/docs/Rbook/html/_Rbook/stat-aov.html）
# 载入 ----------------------------------------------------------------------
library(tidyverse)
d.manu.w <- read.table("e://d.manu.w.txt",
header = T)
# ctrl + 点击 可以View表格
# 处理 ----------------------------------------------------------------------
# 宽变长\\转化因子变量\\排序
d.manu <- d.manu.w %>%
pivot_longer(A:C,
names_to = "grp",
values_to = "y") %>%
mutate(grp = factor(grp, levels=c("A", "B", "C"))) %>%
arrange(grp)
# pivot_longer创建新数据框，mutate为其创建新的变量，并覆盖同名变量
# 建模 ----------------------------------------------------------------------
knitr::kable(d.manu) # 绘制表格
d.manu %>% select(1) %>% table(grp = .) # 分组结果表面，实验为均衡设计
aov.manu <- aov(y ~ grp, data = d.manu)
aov.manu %>% summary()
# 主效应grp（分组）的F检验的p值为0.004，表明
# 95%检验水平下（α=0.05），分组效应显著，各组之间有显著差异。
# 检验 ----------------------------------------------------------------------
# 我们对于结果的信心依赖于做统计检验时数据满足假设条件的程度
# 单因素方差分析中，我们假设因变量服从正态分布、各组方差相等
car::qqPlot(lm(y ~ grp, data = d.manu),
simulate=TRUE, main="Q-Q Plot", labels=FALSE)
# 注意qqPlot()要求用lm()拟合
# 数据落在95%的置信区间范围内，说明满足正态性假设
bartlett.test(y ~ grp, data = d.manu)
# P值不显著，不拒绝各组方差相等的原假设，通过方差齐性检验
# 结果表明各组方差并没有显著不同
# Bartlett检验：通过求取不同组之间的卡方统计量，判断组间方差是否相等。
# 该方法极度依赖于数据的正态分布假设，否则结果偏差很大。
# 其他检验有Fligner-Killeen检验（fligner.test()函数）和Brown-Forsythe检验（HH::hov()函数）
aov.manu %>% car::outlierTest()
# 方差齐性分析对离群点非常敏感
# 当p>1时将产生NA，从检测结果来看，并没有证据说明数据含有离群点
# 因此根据Q-Q图、Bartlett检验和离群点检验，该数据似乎可以用ANOVA模型拟合得很好。
# 多重比较 --------------------------------------------------------------------
# 为了找到各组两两之间是否有显著差异，可以进行两两的独立两样本t检验，
# 但这样不能利用共同的模型参数，进行多次重复检验也会使得总第一类错误概率变得比较高，发生过度拟合。
#
# 进行多个假设检验（如均值比较）的操作称为“多重比较”，多次检验会使得总第一类错误概率增大。
# 以两两比较为例，当比较次数很多时，即使所有的组之间都没有真实的差异，很多个两两比较也会找到差异最大的一对，使得发现的的显著差异有很大可能性不是真实存在的。
#
# 为此，可以进行一些调整， 使得报告的检验p值能够控制总第一类错误概率。
# multcomp包的glht()函数可以对方差分析结果进行多重比较并控制总错误率，其中一种方法是利用J·W·Tukey的HSD方法，能将所有各对平均值同时比较
d.manu %>% group_by(grp) %>% summarise(mean = y %>% mean(),
sd = y %>% sd())
# C组的工作效率显著低于A组和B组， A、B两组之间差异不明显。
# group_by %>% plot
# geom_boxplot
# boxplot
# gplots::plotmeans
# TukeyHSD %>% plot
# multcomp::glht
d.manu %>% group_by(grp) %>% plot(main = "箱线图") # 分组数据自动绘制箱线图
boxplot(y ~ grp,
data = d.manu,
main = "箱线图") # 分组数据箱线图
ggplot(data = d.manu,
mapping = aes(x = grp, y = y)) +
geom_boxplot() +
ggtitle(label = "箱线图")
gplots::plotmeans(formula = y ~ grp,
data = d.manu,
main = "组均值图")
aov.manu %>% TukeyHSD()
aov.manu %>% TukeyHSD() %>% plot()
# AB组间均值差异不显著，p值不显著，不拒绝组间均值相等的原假设
# 成对比较图形中，AB置信区间包含0，说明差异不显著
library(multcomp)
aov.manu %>% multcomp::glht(model = .,
linfct = mcp(grp = "Tukey")) %>% summary()
aov.manu %>% multcomp::glht(model = .,
linfct = mcp(grp = "Tukey")) %>% plot()
aov.manu %>% multcomp::glht(model = .,
linfct = mcp(grp = "Tukey")) %>%
cld(level=.05) %>% plot() # 默认95%CI
# AB组间均值差异不显著，p值不显著，不拒绝组间均值相等的原假设
# Tukey HSD检验的结果显示，0.05水平下，A和B没有显著差异，C与A、B均有显著差异
协方差分析
# 载入 ----------------------------------------------------------------------
library(tidyverse)
data(litter, package="multcomp") # ctrl + click
litter %>% str()
litter %>% tibble()
litter %>% count(dose) %>% knitr::kable()
litter$dose %>% table(dose = .) %>% knitr::kable()
# 建模 ----------------------------------------------------------------------
aov.lit1 <- aov(weight ~ gesttime + dose,
data = litter)
aov.lit1 %>% summary()
# 排除协变量gesttime影响后， 在0.05水平下， 四个处理组之间有显著差异
# 多重比较 --------------------------------------------------------------------
library(multcomp)
aov.lit1 %>% multcomp::glht(
linfct = mcp(
dose = rbind("服药与不服药比较" = c(3, -1, -1, -1))
)
) %>% summary()
# 假设检验的t统计量（2.581）在p<0.05水平下显著,服药的和不服药组的均值有显著差异
# 检验 ----------------------------------------------------------------------
# ANCOVA与ANOVA相同，都需要正态性和同方差性假设。另外，ANCOVA还假定回归斜率相同。
# car::qqPlot()
# bartlett.test()
# car::outlierTest()
# aov(*)
# ANCOVA模型包含怀孕时间×剂量的交互项时，可对回归斜率的同质性进行检验。
# 交互效应若显著，则意味着时间和幼崽出生体重间的关系依赖于药物剂量的水平。
litter %>% aov(data = .,
formula = weight ~ gesttime*dose) %>% summary()
# ~ x协*x自 包含了x协、x自 以及二者交互项对因变量的回归斜率
# 可以看到交互效应不显著，支持了斜率相等的假设。
# 若假设不成立，变换协变量或因变量，或使用能对每个斜率独立解释的模型，
# 或使用不需要假设回归斜率同质性的非参数ANCOVA方法，sm::sm.ancova()
# 可视化 ---------------------------------------------------------------------
# group_by() %>% plot()
# ggplot + geom_boxplot
# boxplot()
#
# gplots::plotmeans()
# TukeyHSD()
# multcomp::glht %>% cld %>% plot#
# HH::ancovaplot
HH::ancovaplot(weight ~ gesttime + dose, data = litter)
多因素方差分析
# 函数总结 --------------------------------------------------------------------
# 处理数据
pivot_longer()
mutate(factor)
arrange()
count()
summarise()
# levels()查看因子变量各水平，table()查看频数，hist(as.numeric(unlist))绘制直方图
# 建模判断
aov() %>% summary() # 一元因变量（任意自变量、协变量个数）、多元因变量
manova() %>% summary() # 多元方差分析
manova() %>% summary.aov() # 对多元因变量分别做单因素分析
effects::effect() # 去除协变量效应后的组均值
# 其他建模
lm()
# 线性模型要求预测变量是数值型，lm()用一系列与因子水平相对应的数值型对照变量来代替因子
# 对照处理用于无序因子，正交多项式用于有序因子，通过contrasts()函数查看它的编码过程
anova()
car::Anova()
ancova(~*) # 不需要假设回归斜率同质性，允许斜率和截距项依据组别而发生变化
sm::sm.ancova(~*) # 不需要假设回归斜率同质性的非参数ANCOVA方法
rrcov::Wilks.test(method="mcd") # 稳健单因素MANOVA
vegan::adonis() # 非参数MANOVA
# 初级可视化（均值、方差、分布）
group_by() %>% plot()
boxplot()
ggplot + geom_boxplot
gplots::plotmeans()
# 高级可视化
gplots::plotmeans() # 单因素或交互效应
interaction.plot() # 交互效应
boxplot() # 交互效应的不同侧面
HH::ancovaplot() # 含协变量
HH::interaction2wt() # 含多因素、交互效应
# 多重比较
TukeyHSD()
multcomp::glht(linfct = mcp("Tukey"))
multcomp::glht(linfct = mcp(rbind))
# 假设检验
car::qqPlot(lm)
car::qqPlot(
qchisq(ppoints(nrow), df = ncol),
mahalanobis(cbind, colMeans, cov)) %>% # 检验多元正态性
identify(labels = row.names) # 识别异常点，点击返回索引值
bartlett.test
fligner.test
HH::hov
Box’s M Test # 检验方差-协方差矩阵同质性，即各组的协方差矩阵是否相同。
# 该检验对正态性假设很敏感，导致在大部分案例中直接拒绝同质性假设。
car::outlierTest
mvoutlier::ap.plot() # 多元离群点
aov(~*) # 回归斜率的同质性
# 载入、处理 -------------------------------------------------------------------
library(tidyverse)
data(ToothGrowth, package="datasets") # ctrl + click
ToothGrowth %>% str()
ToothGrowth$dose %>% table(dose = .)
d.tooth <- ToothGrowth %>% mutate(
supp = factor(supp, levels=c("OJ", "VC")),
dose = factor(dose, levels=c(0.5, 1.0, 2.0))
) # 创建并覆盖
d.tooth %>% count(supp,dose) %>% knitr::kable() # 均衡设计
# 建模 ----------------------------------------------------------------------
# 由上可知，这个数据结果来自均衡、重复测量的实验设计
aov.to1 <- aov(len ~ dose + supp + dose:supp,
data = d.tooth)
aov.to1 %>% summary()
# 在0.05水平下， 两个主效应和交互作用效应都显著
原书：R in action
# 单因素方差分析 -------------------------------------------------------
library(tidyverse)
library(multcomp)
cholesterol %>% tibble()
cholesterol %>% attach()
trt %>% table(treatment = .) # 各组样本大小
# 每10个患者接受其中一个药物疗法
# 各组均值、标准差
response %>% aggregate(by = trt %>% list(),
FUN = "mean")
response %>% aggregate(by = trt %>% list(),
FUN = "sd")
cholesterol %>% group_by(trt) %>%
summarise(mean = response %>% mean(),
sd = response %>% sd())
# 均值显示drugE降低胆固醇最多，而1time降低胆固醇最少
# 各组的标准差相对恒定，在2.88到3.48间浮动
# 检验组间差异（ANOVA）
fit <- aov(response ~ trt, data = cholesterol)
fit %>% summary()
# ANOVA对治疗方式的F检验非常显著，说明五种疗法的效果不同
# 绘制各组均值及其置信区间
library(gplots)
plotmeans(response ~ trt,
data = cholesterol,
xlab = "治疗方案",
ylab = "治疗效果",
main = "❥\\t95%置信区间的均值\\t❥")
# 带有95%的置信区间的各疗法均值，可以清楚看到它们之间的差异
cholesterol %>% detach()
# 多重均值比较：对各组均值差异的成对检验 -----------------------------------------------------
# 多重比较方法是一个复杂但正迅速发展的领域，
# 可参考Bretz、Hothorn和Westfall的Multiple Comparisons Using R（2010）一书
par(las = 1) # 刻度线上的标签 0平行 2垂直 1向上 3向左
par(mar = c(4,7,2,1)) # 边距大小：底部、左侧、顶部和右侧
fit %>% TukeyHSD()
fit %>% TukeyHSD() %>% plot()
# 1time和2times的均值差异不显著（p=0.138）
# 而1time和4times间的差异非常显著（p<0.001）
# 图形中置信区间包含0的疗法说明差异不显著（p>0.5）#
par(mar = c(4,4,6,.5))
fit %>% multcomp::glht(linfct = mcp(trt = "Tukey"))
fit %>% multcomp::glht(linfct = mcp(trt = "Tukey")) %>%
multcomp::cld(level = .05) %>% plot(col = "lightgrey")
# multcomp包中的glht()函数适用于本例中线性模型，也适用于广义线性模型
# cld()函数的level选项设置了显著水平（0.05，即本例中的95%的置信区间）
# plot(cld(glht))比plot(TukeyHSD)更好理解，还提供了各组得分的分布信息
# 有相同字母的组（用箱线图表示）说明均值差异不显著 #
# 1time和2times差异不显著（有相同的字母a）
# 2times和4times差异也不显著（有相同的字母b）
# 而1time和4times差异显著（它们没有共同的字母）#
# 从结果来看，一天四次5mg剂量比一天一次20mg剂量效果更佳，
# 且与候选药物drugD差异不大，但药物drugE比其他所有药物和疗法都更优。
# 假设检验 --------------------------------------------------------------------
# 我们对于结果的信心依赖于做统计检验时数据满足假设条件的程度
# 单因素方差分析中，我们假设因变量服从正态分布、各组方差相等
# 使用Q-Q图来检验正态性假设 #
car::qqPlot(lm(formula = response ~ trt, data = cholesterol),
simulate = TRUE,
main = "QQ图",
labels = FALSE)
# 注意qqPlot()要求用lm()拟合
# 数据落在95%的置信区间范围内，说明满足正态性假设
# 方差齐性检验：Bartlett检验
bartlett.test(response ~ trt, data = cholesterol)
# Bartlett检验表明五组的方差并没有显著不同（p=0.97）
# 其他检验如Fligner-Killeen检验（fligner.test()函数）
# 和Brown-Forsythe检验（HH包中的hov()函数）结果与Bartlett检验相同
# 方差齐性分析对离群点非常敏感。利用car包的outlierTest()函数来检测离群点
car::outlierTest(fit) # fit <- aov()
# 没有证据说明胆固醇数据中含有离群点（当p>1时将产生NA）。
# 因此根据Q-Q图、Bartlett检验和离群点检验，该数据似乎可用ANOVA模型拟合得很好
# 单因素协方差检验 ----------------------------------------------------------------
# data(litter, package="multcomp")
litter %>% tibble()
litter %>% attach()
# 各组样本大小
dose %>% unlist() %>% unlist() %>% as.numeric() %>%
hist(xaxt = "n")
axis(1, at = c(1,2,3,4), labels = c(0,5,50,500))
dose %>% table("剂量" = .) / litter %>% nrow()
# 可以看到每种剂量下所产的幼崽数并不相同：
# 0剂量时（未用药）产崽20个，500剂量时产崽17个
# 各组体重描述性统计
weight %>% aggregate(by = dose %>% list(),
FUN = "mean")
weight %>% aggregate(by = dose %>% list(),
FUN = "sd")
litter %>% group_by(dose) %>%
summarise(mean = weight %>% mean(),
sd = weight %>% sd())
# 发现未用药组幼崽体重均值最高（32.3）。
# 单因素ANCOVA
fit <- aov(weight ~ gesttime + dose, # 协变量在前 #
data = litter)
fit %>% summary()
# ANCOVA的F检验表明：(a)怀孕时间与幼崽出生体重相关（p=0.00597）；
# (b)控制怀孕时间，药物剂量与出生体重相关：
# 控制怀孕时间，每种药物剂量下幼崽出生体重均值不同（p=0.04988）。
# 计算去除协变量效应后的组均值 ----------------------------------------------------------
# 由于使用了协变量，获取调整的组均值，即去除协变量效应后的组均值。
fit %>% effects::effect("dose", .)
fit %>% effects::effect("dose", .) %>% plot()
plot(c(32.3, 29.3, 29.9, 29.6),
type = "b",
main = "调整前各组均值")
points(c(32.4, 28.9, 30.6, 29.3),
col = 2,
type = "b")
# 本例中调整的均值与未调整的均值类似，但并非所有的情况都是如此
# 自定义对照的多重均值比较 -----------------------------------------------------
# 剂量的F检验虽然表明了不同的处理方式幼崽的体重均值不同，
# 但无法告知我们哪种处理方式与其他方式不同。
# 同样，我们使用multcomp包来对所有均值进行成对比较。
# 另外，multcomp包还可以用来检验用户自定义的均值假设。
# 假定你对未用药条件与其他三种用药条件影响是否不同感兴趣
fit %>% multcomp::glht(.,
linfct = mcp(
dose = rbind(
"no drug vs. drug" = c(3, -1, -1, -1)))) %>%
summary()
# 对照c(3, -1, -1, -1)设定第一组和其他三组的均值进行比较 #
# 假设检验的t统计量（2.581） 在p<0.05水平下显著，
# 因此，可以得出未用药组比其他用药条件下的出生体重高的结论。
# 其他对照可用rbind()函数添加（详见help(glht)）。
# 评估检验的假设条件 ---------------------------------------------------------------
# ANCOVA与ANOVA相同，都需要正态性和同方差性假设；
car::qqPlot(lm(formula = weight ~ gesttime + dose,
data = litter),
simulate = TRUE,
main = "QQ图",
labels = FALSE)
bartlett.test(weight ~ dose, data = litter)
bartlett.test(weight ~ gesttime, data = litter)
bartlett.test(gesttime ~ dose, data = litter)
fit %>% car::outlierTest()
# 另外，ANCOVA还假定回归斜率相同：
# 本例中，假定四个处理组通过gesttime来预测weight的回归斜率都相同
# ANCOVA模型包含gesttime × dose的交互项时，可对回归斜率的同质性进行检验
# 交互效应若显著，则意味着gesttime和weight间的关系依赖于药物剂量的水平
fit2 <- aov(formula = weight ~ gesttime*dose,
data = litter)
fit2 %>% summary() # summary.aov() 效果相同
# 可以看到交互效应不显著，支持了斜率相等的假设 #
# 若假设不成立，可以尝试变换协变量或因变量，
# 或使用能对每个斜率独立解释的模型，
# 或使用不需要假设回归斜率同质性的非参数ANCOVA方法。
# sm包中的sm.ancova()函数为后者提供了一个例子。
?sm::sm.ancova()
sm::sm.ancova(x = weight,
y = gesttime,
group = dose,
data = litter,
model = "none")
# 结果可视化 -------------------------------------------------------------------
library(HH)
ancova(response ~ trt + trt, data = cholesterol)
ancova(weight ~ gesttime + dose, data = litter) # groups = dose 效果一样
# 用怀孕时间来预测出生体重的回归线相互平行，只是截距项不同。
# 随着怀孕时间增加，幼崽出生体重也会增加。
# 另外，还可以看到0剂量组截距项最大，5剂量组截距项最小。
以上是关于R语言 方差分析的主要内容，如果未能解决你的问题，请参考以下文章