R語言--史上最全方差分析

單因素

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

重復測量