《R语言与统计分析》-方差分析

####单因子方差分析####
aov(formula, data = NULL, projections = FALSE, qr = TRUE,
    contrasts = NULL, ...)
#在单因子方差分析中，formula为x~A
#aov(formula)等价于anova(lm(formula))

####均值的多重比较####
####多重t检验-p.adjust####
p.adjust(p, method = p.adjust.methods, n = length(p))

p.adjust.methods
# c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY",
#   "fdr", "none")

####多重t检验-pairwise.t.test####
#当比较次数较多时，Bonferroni方法较好
pairwise.t.test(x, g, p.adjust.method = p.adjust.methods,
                pool.sd = !paired, paired = FALSE,
                alternative = c("two.sided", "less", "greater"),
                ...)

####同时置信区间：Tukey法-TukeyHSD####
TukeyHSD(x, which, ordered = FALSE, conf.level = 0.95, ...)

#例:
sales<-data.frame(
  X=c(23,19,21,13,24,25,28,27,20,18,19,15,
      22,25,26,23,24,23,26,27),
  A=factor(rep(1:5,c(4,4,4,4,4)))
)

summary(aov(X~A,sales))

TukeyHSD(aov(X~A,sales))

####方差齐性检验####
#要进行方差分析，有三个先行条件：
#1.可加性
#2.独立正态性
#3.方差齐性
####Bartlett检验####
bartlett.test(x, g, ...)
bartlett.test(formula, data, subset, na.action, ...)

####Levene检验####
library(car)
levenetest(x,group)

#单因子方差分析还可使用函数
oneway.test(formula, data, subset, na.action, var.equal = FALSE)
#各水平数据的方差相等时（ var.equal = T），等同于aov()
#方差不相等时，（ var.equal = F），使用Welch的近似方法进行方差分析

#各水平的分布未知时，利用
kruskal.test(x, g, ...)
#进行方差分析

####双因子方差分析####
####无交互作用的方差分析####
aov(formula, data = NULL, projections = FALSE, qr = TRUE,
    contrasts = NULL, ...)
#formula:x~A+B

#例：
#A1,A2,A3,A4
#B1,B2,B3,B4,B5,B6
juice<-data.frame(
  X=c(0.05,0.46,0.12,0.16,0.84,1.3,0.08,0.38,0.4,
      0.1,0.92,1.57,0.11,0.43,0.05,0.1,0.94,1.1,
      0.11,.44,.08,.03,.93,1.15),
  A=gl(4,6),
  B=gl(6,1,24)
)
#gl(n, k, length = n*k, labels = seq_len(n), ordered = FALSE)
#n是水平数，k是每一水平上的重复次数，length是总观测值数，
#ordered为各水平是否先排序
juice.aov<-aov(X~A+B,data=juice)
summary(juice.aov)
# Df Sum Sq Mean Sq F value Pr(>F)    
# A            3  0.057  0.0190   1.629  0.225    
# B            5  4.902  0.9804  83.976  2e-10 ***
#   Residuals   15  0.175  0.0117                   
# ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#对两个因素分别进行齐性检验
bartlett.test(X~A,juice)#对因素A
# Bartlett test of homogeneity of variances
# 
# data:  X by A
# Bartlett's K-squared = 0.26802, df = 3, p-value = 0.9659

bartlett.test(X~B,juice)#对因素B
# Bartlett test of homogeneity of variances
# 
# data:  X by B
# Bartlett's K-squared = 17.422, df = 5, p-value = 0.003766

####有交互作用的方差分析####
aov(formula, data = NULL, projections = FALSE, qr = TRUE,
    contrasts = NULL, ...)
#formula:x~A+B+A:B

#例：
#A,B,C,D
#I，II，III
rats<-data.frame(
  Time=c(.31,.45,.46,.43,.82,1.1,.88,.72,.43,.45,.63,.76,
      .45,.71,.66,.62,.38,.29,.4,.23,.92,.61,.49,1.24,
      .44,.35,.31,.4,.56,1.02,.71,.38,.22,.21,.18,.23,
      .3,.37,.38,.29,.23,.25,.24,.22,.3,.36,.31,.33),
  Toxicant=gl(3,16,48,labels=c("I","II","III")),
  Cure=gl(4,4,48,labels = c("A","B","C","D"))
)
op<-par(mfrow=c(1,2))
plot(Time~Toxicant+Cure,rats)

with(rats,interaction.plot(Toxicant,Cure,Time,trace.label = "Cure"))
with(rats,interaction.plot(Cure,Toxicant,Time,trace.label = "Toxicant"))
#两图中的曲线没有明显的相交情况，初步认为两个因素没有交互作用

rats.aov<-aov(Time~Toxicant*Cure,rats)
summary(rats.aov)

library(car)
levenetest(rats$Time,rats$Toxicant)
levenetest(rats$Time,rats$Cure)
bartlett.test(Time~Toxicant,rats)
bartlett.test(Time~Cure,rats)

####协方差分析####
library(HH)
ancova(formula, data.in = NULL, ...,x, groups, ...)

#formula是协方差公式，data.in是数据框，x为协方差分析中的协变量
#例：
feed<-rep(c("A","B","C"),each=8)
Weight_Initial<-c(15,13,11,12,12,16,14,17,
                  17,16,18,18,21,22,19,18,
                  22,24,20,23,25,27,30,32)
Weight_Increment<-c(85,83,65,76,80,91,84,90,
                    97,90,100,95,103,106,99,94,
                    89,91,83,95,100,102,105,110)
data_feed<-data.frame(feed,Weight_Initial,Weight_Increment)

#检验初始体重与饲料组的交互作用是否显著：
ancova(Weight_Increment~Weight_Initial*feed,
       data.in= data_feed)

ancova(Weight_Increment~Weight_Initial+feed,
       data.in= data_feed)