方差分析

借助实例，理解R代码，形成自己的代码库，方便自己使用。
部分参考：https://blog.csdn.net/myl1992/article/details/45865837

---

---

前言

R语言是重要的开源的软件，常用于统计学分析，回归分析和时间序列分析就属于此类。

---

一、单因素方差分析

数据导入

一般地，在R中的数据不是因子型，所以进行数据转化是有必要的。

# 将数据转化为因子
factor(x) 
# factor代表因子水平数，factor_num代表，label代表标签
gl(factor, factor_num, labels = c("s", "m"))

各个因子水平

library("ggpubr")
ggline(data, x="factor", y="weight",
       order=c("1","2","3"))

正态性检验

group <- split(y, x)#将数据按水平分割成list
# 注意前面后面一个为factor，不可颠倒
unlist(lapply(group,function(x){shapiro.test(x)$p.value}))#做每个水平下的正假定假定

观察Q-Q图

library(car)
qqPlot(group[[1]])

方差齐性检验

#方差齐性检验,大于0.05不能拒绝方差齐性
leveneTest(weight~factor, data=data)

#观测离群点,Bonferroni p<0.05则没有离群点
outlierTest(lm(weight~factor, data = data))

方差分析

# 上述步骤想省去，就去干吧
library(car)
library(carData)

# 拟合函数
fit <- aov(weight~factor, data=data)
summary(fit)

# 缺少其他项，采用自定义函数补充,失拟检验
anova.tab <- function(fit){
  tab <- summary(fit)
  k<-length(tab[[1]])-2
  
  # 计算总和
  tdf = sum(tab[[1]][,1])
  SE = sum(tab[[1]][,2])
  tab[[1]]["Total",] <- c(tdf,SE,SE/tdf,rep(NA,2))
  
  # 失拟检验
  F.value = tab[[1]]["Residuals  ",3]/(SE/tdf)
  tab[[1]]["Residuals  ",c(4,5)] <- c(F.value,1-pf(F.value,tab[[1]]["Residuals  ",1],tdf))
  return(tab)
}

anova.tab(fit)

方差不齐

oneway.test(weight~factor, var.equal = F)

二、复杂的方差分析

双因素方差分析

data <- ToothGrowth
data$dose <- factor(data$dose, levels=c(0.5,1,2),labels=c("D0.5","D1","D2"))

fit <- aov(len~supp*dose,data=data)
summary(fit)

带有协变量的方差分析

library(car)
library(carData)

# 带有协变量
fit <- aov((weight~x+factor, data=data))#x为协变量
summary(fit)

# 去除协变量效应后的组均值
effect("x", aov(weight~x+factor, data = data))

# 观察正态性和同方差性假设
fit <- aov(weight~x*factor, data = data)
summary(fit)