R语言实战-第九章 方差分析

#第九章方差分析
#需要的packages:car gplots HH rrcov multicomp effects MASS mvotlier
#单因素方差分析
#数据集来源multcomp包的cholesterol数据集
library(multcomp)
attach(cholesterol)
table(trt)
aggregate(response,by=list(treatment=trt),FUN=mean)
aggregate(response,by=list(treatment=trt),FUN=sd)
fit <- aov(response~trt)
summary(fit)
library(gplots)
gplots::plotmeans(response~trt,xlab = "treatment",ylab="response")
#plotmeans可以绘制带有置信区间的组均值图形

#多重比较 明确两两之间的显著性差异
#方法1
TukeyHSD(fit)
par(las=1)# =2时表示旋转轴标签
par(mar=c(5,8,4,2))
plot(TukeyHSD(fit))#图形中置信区间若包含0，表示差异不显著

#方法2
library(multcomp)
par(mar=c(5,4,6,2))
tuk <- multcomp::glht(fit,linfct = mcp(trt="Tukey"))
plot(cld(tuk,level=0.05),col="lightgrey")#level设置显著水平

#评估检验的假设条件
#正态假设性
library(car)
car::qqPlot(lm(response~trt,data = cholesterol),
            simulate=TRUE,main="QQ plot",labels=FALSE)
#数据均落在置信区间内，满足正态性假设

#方差齐性检验
#方法1 bartlett.test
bartlett.test(response~trt,data = cholesterol)
#p值表明五组数据方差没有显著不同，满足方差齐性

#方法2 fligner.test
fligner.test(response~trt,data = cholesterol)

#方法3 HH::hov

#利群点的检测 方差齐性对离群点非常敏感
library(car)
outlierTest(fit)

#综上 根据QQ图 bartlett检验 离群点检验证明数据可以用ANOVA模型模拟的很好
data(litter,package="multcomp")
attach(litter)
aggregate(weight,by=list(dose),mean)
fit <- aov(weight~gesttime+dose)
summary(fit)
#怀孕时间与幼崽出生体重相关，控制怀孕时间，药物剂量与出生体重相关

#如何去除协变量效应后的组均值
library(effects)
effects::effect("dose",fit)

#单因素协方差分析的多重比较
library(multcomp)
contrast <- rbind("no drug vs.drug"=c(3,-1,-1,-1))
summary(glht(fit,linfct=mcp(dose=contrast)))

#评估检验单因素协方差分析检验的假设条件
#正态性 方差齐性 回归斜率相同
library(multcomp)
fit2 <- aov(weight~gesttime*dose,data = litter)
summary(fit2)
#若交互效应显著 则意味着时间和出生体重间的关系依赖于药物剂量
#结果可视化HH::ancova
library(HH)
HH::ancova(weight~gesttime+dose,data=litter)

#双因素方差分析
attach(ToothGrowth)
table(supp,dose)
aggregate(len,by=list(supp,dose),mean)
aggregate(len,by=list(supp,dose),sd)
dose <- as.factor(dose)
fit <- aov(len~supp*dose)
summary(fit)
#主效应和交互响应都十分显著

#可视化双因素方差分析
#method 1
interaction.plot(dose,supp,len,type = "b",
                 col = c("red","blue"),pch = c(16,18),
                 main = "picture")

#method 2
library(gplots)
plotmeans(len~interaction(supp,dose,sep = " "),
          connect=list(c(1,3,5),c(2,4,6)),#设置线条连接方式
          col = c("red","blue"),
          main = "picture 2",
          xlab = "Treatment and dose combination")

#method 3
library(HH)
interaction2wt(len~supp*dose)
aggregate(len,by=list(supp),quantile) #左下角图
aggregate(len,by=list(dose),quantile) #有上角图

#重复测量方差分析及其可视化
#使用数据为基础安装包中的CO2数据集
CO2
str(CO2)
#关注寒带植物 因变量是CO2吸收量 自变量是植物类型(组间因子) 和 CO2浓度(组内因子)
attach(CO2)
conc <- as.factor(conc)
w1b1 <- subset(CO2,Treatment=="chilled")
fit <- aov(uptake~conc*Type+Error(Plant/(conc)))
summary(fit)

#可视化 1
par(las=2)
par(mar=c(10,4,4,2))
with(w1b1,interaction.plot(conc,Type,uptake,
                           type = 'b',col = c("red","blue"),pch = c(16,18),
                           main = "picture 3"))
#可视化 2
boxplot(uptake~Type*conc,data = w1b1,col=c("red","gold"))

#多元方差分析
#使用MASS包中的UScereal数据集
#卡路里 脂肪和糖含量是因变量；货架三水平是自变量
library(MASS)
attach(UScereal)
shelf <- as.factor(shelf)
y <- cbind(calories,fat,sugars)
aggregate(y,by=list(shelf=shelf),FUN=mean)
aggregate(y,by=list(shelf=shelf),FUN=sd)
cov(y)#方差和协变量
fit <- manova(y~shelf)
summary(fit) #是否存在显著差异
summary.aov(fit) #哪些因变量存在显著差异

#评估假设检验
#检验多元正态性
#方法1 QQ图
center <- colMeans(y)
n <- nrow(y)
p <- ncol(y)
cov <- cov(y)
d <- mahalanobis(y,center,cov)
coord <- qqplot(qchisq(ppoints(n),df=p),d)
abline(a=0,b=1)
identify(coord$x,coord$y,labels = row.names(UScereal))
#如果全部点落在斜率为1，截距为0的直线上，表明数据服从多元正态分布
#从图中可以发现，右上角两个点异常，违反了多元正态分布，可以删除这两个点再重新分析

#方法 2 box's M检验
#Box's M 检验验证方差-协方差矩阵同质性(即各组的协方差矩阵相同)
#p 值大于 0.05 即说明各组的协方差矩阵相同
merge <- cbind(shelf,y)
dim(y)
dim(merge)
head(merge)
merge <- as.data.frame(merge)
merge$shelf <- as.factor(merge$shelf)
library(biotools)
boxM(merge[ ,c("calories","fat","sugars")], merge[ ,"shelf"])
#p值小于0.05 各组的协方差矩阵不相同 不符合多元正态性

#方法3 mvoutlier::ap.plot函数检验多元离群点
library(mvoutlier)
outliers <- aq.plot(y)
outliers

#如果多元正态性或方差-协方差均值假设均不满足，可使用稳健单因素 MANOVA检验
#通过 rrcov 包中的 Wilks.test() 函数实现，详情可使用 ?Wilks.test 查看帮助
library(rrcov)
Wilks.test(y,shelf, method = "mcd")
#稳健检验对离群点和违反MANOVA假设的情况不敏感，并再次验证了自变量对3个因变量显著差异

#用回归做ANOVA(模糊 再理解一下)
#用单因素ANOVA为例
library(multcomp)
levels(cholesterol$trt)
fit.aov <- aov(response~trt,data = cholesterol)
summary(fit.aov)
#回归方法
fit.lm <- lm(response~trt,data = cholesterol)
summary(fit.lm)