rm(list = ls())
install.packages(\'multcomp\')
library(multcomp)
data(\'cholesterol\')
head(cholesterol)
data <- cholesterol
#正态性分析####
shapiro.test(cholesterol$response) #data#观测值
#不满足正态性
#1.当偏态分布时，且偏态分布不严重时，仍可以使用方差分析，应当报告偏度、峰度
fBasics::basicStats(x)
#2.使用非参数检验Kruskal-Waliis检验
kruskal.test()

#方差齐性检验####
bartlett.test(response~trt,data=cholesterol) #观测值~分组

library(car)
leveneTest(response ~ trt,data=cholesterol, center=mean)#观测值~分组
                                            #center=mean为原始levene,median更稳健

#不满足方差齐性####
#1.应先处理异常值
#2.使用Welch\'s anova方法，两两比较使用Game-Howell方法(尤其各组例数不同)
oneway.test(x~g,data=,var.equal = F)
#Game-Howell方法:http://aoki2.si.gunma-u.ac.jp/R/src/tukey.R

#单因素方差分析

anova <- with(cholesterol,aov(response ~ trt))
summary(anova)

OR

anova <- oneway.test(response ~ trt , cholesterol, var.equal = T)
anova
####描述统计
by(cholesterol$response,cholesterol[,"trt"],summary) 
aggregate(cholesterol$response,by=list(cholesterol$trt),FUN=sd)
fBasics::basicStats(x)

#自定义函数
mystats <- function(x,na.omit=FALSE){   #一个函数返回x的多个描述性统计量
  if (na.omit)
    x <- x[!is.na(x)]
  m <- mean(x)
  n <- length(x)
  s <- sd(x)
  return(c(n=n,mean=m,stdev=s))
}
dstats <- function(x)sapply(x, mystats)

options(digits = 3)
by(data[2],data$trt,dstats)  
#data$response 和 data[,2]数据均为横向显示，而data[2]则为纵向
by(mtcars[1],mtcars$am,dstats)

#画图,建议用graphpad prism画
install.packages(\'gplots\')
library(gplots)
attach(cholesterol)
plotmeans(response ~ trt)
detach(cholesterol)

##单因素方差分析与单因素线性回归的关系
anolm <-lm(response ~trt ,data = cholesterol)  #trt为有4个哑变量，把1time当做reference
summary(anolm) #F值相等

rm(list = ls())
data("ToothGrowth")
head(ToothGrowth)

table(ToothGrowth$supp,ToothGrowth$dose)  #查看各组例数

#各组均数 标准差
aggregate(ToothGrowth$len,by=list(ToothGrowth$supp,ToothGrowth$dose),mean)
aggregate(ToothGrowth$len,by=list(ToothGrowth$supp,ToothGrowth$dose),sd)

#双因素方差分析
twoanova <- aov(len ~ supp*dose,data = ToothGrowth)#+只有主效应| ：只有交互效应|*都有
summary(twoanova)

#画图
interaction.plot(ToothGrowth$dose,ToothGrowth$supp,ToothGrowth$len,
                 type=\'o\',
                 col=c(\'red\',\'blue\'))

rm(list = ls())

data("airquality")
head(airquality)

#两两比较
#两两比较方法很多，具体的选择看个人
#最后的个人小结(仅代表个人观点): 如果各组例数相等,建议首选Tukey法;
#如果例数不等,建议首选Scheffe法(如果比较组数不多,如3组,Bonferroni法也可以作为首选);
#如果要分别比较每个试验组与对照组,建议采用Dunnett法;

model <- aov(Temp~as.factor(Month),data=airquality)
summary(model)
##Tukey

TukeyHSD(model)
plot(TukeyHSD(model))#除了横跨虚线的均有统计学意义


##LSD医学常用
install.packages("agricolae")
library(agricolae)

out <- LSD.test(y = model, trt = "as.factor(Month)", p.adj = "none" ) 
                                            #多重比较，不矫正P值,==t检验
                                            #p.adj="bonferroni"
out$groups # 查看每个组的label,group不一样即有差别，

plot(out) # 可视化展示


##如果例数不等,建议首选Scheffe法
library(agricolae)

out <-scheffe.test (y = model,trt = "as.factor(Month)")

out$group 

plot(out) 


##如果要分别比较每个试验组与对照组,建议采用Dunnett法
library(multcomp)

out <- glht(model1, 
            linfct = mcp(\'as.factor(Month)\' = \'Dunnett\'), #mcp可以是不同的方法
            alternative = \'two.side\')

summary(out)

rm(list = ls())
library(multcomp)
data(litter)
head(litter)
###协方差分析指的是某一变量A影响了自变量对响应变量的效应，变量A被称作协变量
#前提要求：协变量为连续性变量 | 响应变量正态且方差齐|协变量对响应变量的影响呈线性关系


#正态和方差齐
table(litter$dose)
aggregate(litter$weight,by=list(litter$dose),FUN=mean)
aggregate(litter$weight,by=list(litter$dose),FUN=sd)
shapiro.test(litter$weight)
bartlett.test(litter$weight~litter$dose)
library(car)
leveneTest(weight ~ dose, data=litter)
qqPlot(lm(weight ~ gesttime + dose,data=litter),simulate = T,labels=F)
#协变量对响应变量的线性关系可以通过交互项的显著性判断
interaction <- aov(weight ~ dose * gesttime, data = litter)
summary(interaction)#交互项不显著


#协方差分析  
#必须为 响应变量 ~ 协方差 + 控制因素
ancova <- aov(weight ~ gesttime +dose, data = litter)
summary(ancova)#gesttime显著，说明对响应变量产生了影响
#查看调整后的均值
install.packages(\'effects\')
library(effects)
effect(\'dose\',ancova)

#两两比较
K <- rbind(contrMat(table(litter$dose), "Tukey")) 
            #根据‘Tukey’自动生成一个矩阵

#Tukey为两两比较，有6种情况；Dunnet则是实验组与对照组比较，为3种（3行）

          
mc1 <- glht(ancova, linfct = mcp(dose = K),#根据K矩阵计算，K可自定义
            alternative = "less")
summary(mc1)

#自定义两两比较
K <- rbind(\'5-500\'=c(0,1,0,-1),
           \'0-500\'=c(1,0,0,-1))

mc2 <- glht(ancova, linfct = mcp(dose = K), alternative = "less")
summary(mc2)

###根据不同方法调整P
summary(mc1, test = univariate())
summary(mc2, test = adjusted("bonferroni"))


##协变量对响应变量的线性关系可以通过交互项的显著性判断---之可视化判断####
data1 <- subset(litter,dose==0)
data2 <- subset(litter,dose==5)
data3 <- subset(litter,dose==50)
data4 <- subset(litter,dose==500)
par(mfrow=c(2,2))
plot(data1$gesttime,data1$weight,main = dose0)
abline(lm(weight~gesttime,data=data1))

plot(data2$gesttime,data2$weight,main = dose5)
abline(lm(weight~gesttime,data=data2))

plot(data3$gesttime,data3$weight,main = dose50)
abline(lm(weight~gesttime,data=data3))

plot(data4$gesttime,data4$weight,main = dose500)
abline(lm(weight~gesttime,data=data4))

#日照时间对农作物的 糖分、维生素、脂肪含量
#课外阅读量对 学生语文、数学、历史成绩的影响
#小学生每日VC摄入与其身高、胸围、BMI的影响


#多元方差分析应用条件和要求：
#多元方差分析指的是响应变量不止一个，且存在高度线性相关
#对某个概念需要多个维度去进行描述，否则难以反应对象真实情况
#要求自变量为分类变量

rm(list = ls())

library(MASS)
data("UScereal")
head(UScereal)
#货架对各层营养含量的影响
y <- with(UScereal, cbind(calories,fat,sugars))
mav <- manova(y ~ shelf, data=UScereal)
#两控制变量，3响应变量
#mav <- manova(y ~ shelf + control +shelf*control, data=UScereal)

summary(mav)
##更多前提假设参考
#https://www.sohu.com/a/213290294_489312