R之方差分析与秩和Kruskal-Wallis

一、单因素方差分析与协方差分析理论部分

1、方差分析基本思想

MS组内=MS组间

如何判断是否有统计学意义呢？

将计算得到的F值与F分布的界值相比较

2、方差分析用途

检验两个或多个样本均数间的差异有无统计学意义；注意：两个样本均数的比较可以采用t检验或 F检验，两个以上样本均数的比较只能用F检验
回归方程的线性假设检验；
检验两个或多个因素间有无交互作用

3、方差分析应用条件

各个样本是相互独立的随机样本，变量是连续性变量
各个样本来自正态总体
各个处理组的总体方差方差相等，即方差齐

4、不满足方差分析应用条件时的处理方法

进行变量变换，以达到方差齐或正态的要求
采用非参数法（秩和检验）
使用近似F检验

5、完全随机设计的方差分析

6、配对设计与配伍（区组）设计

7、配对设计与配伍（区组）设计举例

8、多个样本均数的两两比较

由于涉及的对比组数大于2，就不能应用前面介绍的t检验。若仍用前述前述的t检验方法，对每两个对比组作比较，会使犯第I类错误(拒绝了实际上成立的H0所犯的错误)的概率α增大，即可能把本来无差别的两个总体均数判为有差别

9、因素和水平

降压药物对收缩压SBP的影响，因素是不同疗法与不同时间

二、单因素方差分析与协方差分析之R语言实现

ANOVA和回归方法都是独立发展，但从函数形式上，它们都是广义线性模型的特例。回归中的lm()函数也能分析ANOVA模型。两个函数的结果是等同的。

aov()函数的语法为：

 aov(formula, data=dataframe)

表9-4列举了表达式中可以使用的特殊符号。表9-4中的y是因变量，字母A、 B、 C代表因子。

注意交互作用的表达

1、常见研究设计的表达式

表9-5列举了一些常见的研究设计表达式。在表9-5中，小写字母表示定量变量，大写字母表示组别因子， Subject是对被试者独有的标识变量

备注：Error(Subject/A)指误差项

1.1 表达式中各项的顺序

表达式中效应的顺序在两种情况下会造成影响： (a)因子不止一个，并且是非平衡设计； (b)存在协变量。出现任意一种情况时，等式右边的变量都与其他每个变量相关。此时，我们无法清晰地划分它们对因变量的影响。例如，对于双因素方差分析，若不同处理方式中的观测数不同，那么模型y ~ AB与模型y ~ BA的结果不同

样本大小越不平衡，效应项的顺序对结果的影响越大。一般来说，越基础性的效应越需要放在表达式前面。具体来讲，首先是协变量，然后是主效应，接着是双因素的交互项，再接着是三因素的交互项，以此类推。对于主效应，越基础性的变量越应放在表达式前面，因此性别要放在处理方式之前。有一个基本的准则：若研究设计不是正交的（也就是说，因子和/或协变量相关），一定要谨慎设置效应的顺序。

1.2 顺序很重要！

含因子A、 B和因变量y的双因素不平衡因子设计，有三种效应： A和B的主效应， A和B的交互效应。假设你正使用如下表达式对数据进行建模： Y ~ A + B + A:B 有三种类型的方法可以分解等式右边各效应对y所解释的方差。

类型Ⅰ （序贯型）效应根据表达式中先出现的效应做调整。 A不做调整， B根据A
调整， A:B交互项根据A和B调整
类型Ⅱ（分层型）效应根据同水平或低水平的效应做调整。 A根据B调整， B依据
A调整， A:B交互项同时根据A和B调整
类型Ⅲ（边界型）每个效应根据模型其他各效应做相应调整。 A根据B和A:B做调
整， A:B交互项根据A和B调整
R默认调用类型I方法，其他软件（比如SAS和SPSS）默认调用类型Ⅲ方法，按照R语言，第一个因素不回被调整，因而注意顺序。
请注意car包中的Anova()函数（不要与标准anova()函数混淆）提供了使用类型，Ⅱ或类型Ⅲ方法的选项，而aov()函数使用的是类型I方法。若想使结果与其他软件（如SAS和SPSS）提供的结果保持一致，可以使用Anova()函数，细节可参考help(Anova,package=“car”)。

2、单因素方差分析

单因素方差分析用于比较分类因子定义的两个或多个组别中的因变量均值。

以multcomp 包中的 cholesterol 数据集为例（取自Westfall、 Tobia、 Rom、Hochberg， 1999）， 50个患者均接受降低胆固醇药物治疗（trt）五种疗法中的一种疗法。其中三种治疗条件使用药物相同，分别是20mg一天一次（1time）、 10mg一天两次（2times）和5mg一天四次（4times）。剩下的两种方式（drugD和drugE）代表候选药物。哪种药物疗法降低胆固醇（响应变量）最多呢？分析过程见如下代码：

备注：trt一定要是因子变量，否则要重新定义

3、多重比较

虽然ANOVA对各疗法的F检验表明五种药物疗法效果不同，但是并没有告诉你哪种疗法与其他疗法不同。多重比较可以解决这个问题。例如， TukeyHSD()函数提供了对各组均值差异的成对检验（见如下代码）：

还可以画图：

TukeyHSD(fit)
par(las=1) #可以试一下1有什么效果
par(mar=c(5,8,4,2))
plot(TukeyHSD(fit))
par(opar)library(multcomp)
par(mar=c(5,4,6,2))
tuk <- glht(fit, linfct=mcp(trt="Tukey"))
plot(cld(tuk, level=.05),col="lightgrey")
par(opar)

横轴表示均数差，与0相交（无效线）则没有差别

也可以使用Dunnett检验法：

library(multcomp)
par(mar=c(5,4,6,2))
tuk <- glht(fit, linfct=mcp(trt="Dunnett"))

4、评估方差分析的假设条件

单因素方差分析中，我们假设因变量服从正态分布，各组方差相等。可使用Q-Q图来检验正态性假设：

library(car)
qqPlot(lm(response ~ trt, data=cholesterol),simulate=TRUE, main="Q-Q
Plot", labels=FALSE)

看点是否在虚线以内，一般样本量至少40个以上图性比较准确

注意qqPlot()要求用lm()拟合，落在95%的置信区间内，说明满足正态性假设。

R提供了一些可用来做方差齐性检验的函数。例如，可以通过如下代码来做：

bartlett检验：
bartlett.test(response ~ trt, data=cholesterol)

方差齐性分析对离群点非常敏感。可利用car包中的outlierTest()函数来检测离群点：

library(car)
outlierTest(fit)

5、单因素协方差分析

单因素协方差分析（ANCOVA）扩展了单因素方差分析（ANOVA），包含一个或多个定量协变量。下面的例子来自于multcomp 包中的 litter 数据集（见Westfall et al.， 1999）。怀孕小鼠被分为四个小组，每个小组接受不同剂量（0、 5、 50或500）的药物处理。产下幼崽的体重均值为因变量，怀孕时间为协变量。分析代码如下：

备注：gesttime是数值型变量（协变量），dose是因子变量**（分组效应）**

5.1 校正均值

由于使用了协变量，你可能想要获取调整的组均值，即去除协变量效应后的组均值。可使用effects包中的effects()函数来计算调整的均值：

library(effects)
effect("dose", fit)

5.2 对用户定义对照的多重比较

ibrary(multcomp)
contrast <- rbind("no drug vs. drug" = c(3, -1, -1, -1))
summary(glht(fit, linfct=mcp(dose=contrast))) #fit是前方定义过的，后面只是一下比较标准

对照c(3, -1, -1, -1)设定第一组和其他三组的均值进行比较。假设检验的t统计量（2.581）在p<0.05水平下显著，因此，可以得出未用药组比其他用药条件下的出生体重高的结论。其他对照可用rbind()函数添加

5.3 评估检验的假设条件

ANCOVA与ANOVA相同，都需要正态性和同方差性假设，可以用上节中相同的步骤检验这些假设条件。另外， ANCOVA还假定回归斜率相同。本例中，假定四个处理组通过怀孕时间来预测出生体重的回归斜率都相同。 ANCOVA模型包含怀孕时间×剂量的交互项时，可对回归斜率的同质性进行检验。交互效应若显著，则意味着时间和幼崽出生体重间的关系依赖于药物剂量的水平。代码如下：

library(multcomp)
fit2 <- aov(weight ~ gesttime*dose, data=litter) #还考虑了协变量与分组变量的交互作用
summary(fit2)

三、双因素方差分析

1、双因素方差分析

先来看一个例子：

双因素方差分析中，受试者被分配到两因子交叉类别组中。以基础安装包中的ToothGrowth 数据集为例，随机分配60只豚鼠，分别采用两种喂食方法（橙汁或维生素C），各喂食方法中抗坏血酸含量有三种水平（0.5mg、 1mg或2mg），每种处理方式组合都被分配10只豚鼠。牙齿长度为因变量

代码如下：

attach(ToothGrowth)
table(supp, dose) #table语句的预处理表明该设计是均衡设计（各设计单元中样本大小都相同）# aggregate语句处理可获得各单元的均值和标准差
aggregate(len, by=list(supp, dose), FUN=mean)
aggregate(len, by=list(supp, dose), FUN=sd) #dose变量被转换为因子变量，这样aov()函数就会将它当做一个分组变量，而不是一个数值型协变量
dose <- factor(dose)
fit <- aov(len ~ supp*dose) # fit1<- aov(len~supp+dose+supp:dose)等价 #用summary()函数得到方差分析表，可以看到主效应（supp和dose）和交互效应
summary(fit)

有多种方式可以对双因素方差分析结果进行可视化处理，此处可用interaction.plot()函数来展示双因素方差分析的交互效应。

interaction.plot(dose, supp, len, type="b",col=c("red","blue"), pch=c(16,18),main = "Interaction between Dose and Supplement Type") library(gplots)
plotmeans(len ~ interaction(supp, dose, sep=" "),connect=list(c(1, 3, 5),c(2, 4, 6)), col=c("red","darkgreen"),main = "Interaction Plot with 95% CIs", xlab="Treatment and Dose Combination")
library(HH)
interaction2wt(len~supp*dose)

备注：图不能完全判断是否有交互效应，还是看前方的p值

2、重复测量的方差分析

一般来说，以下几种情况需要考虑重复测量的方法分析：

研究主要目的之一是考察某指标在不同时间的变化情况。
研究个体间变异很大，应用普通研究设计的方差分析时，方差分析表中的误差项值将很大，即计算F值时的分母很大，对反应变量有作用的因素常难以识别。
有的研究中研究对象很难征募到足够多的数量，此时可考虑对所征募到的对象在不同条件下的反应进行测量。

2.1 重复测量的方差分析原理

**基本思想：**仍然应用方差分析的基本思想，将反应变量的变异分解成以下四个部分：研究对象内的变异（即测量时间点或测量条件下的效应）、研究对象间的变异（即处理因素效应）、上述两者的交互作用、随机误差。

因素：

受试者内因素----用于区分重复测量次数的变量
受试者间因素----在重复测量时保持恒定的因素

**分析目的：**一是分析受试者间因素的作用；二是考察随着测量次数的增加，测量指标是如何发生变化的，以及分组因素的作用是否会随时间发生，即是否和时间存在交互作用。

2.2 重复测量的方差分析应用条件

反应变量之间存在相关关系。
反应变量的均数向量服从多元正态分布。
对于自变量的各取值水平组合而言，反应变量的方差、协方差矩阵相等。

2.3 案例1

因变量是二氧化碳吸收量（uptake），单位为ml/L，自变量是植物类型Type（魁北克VS.密西西比）和七种水平（95~1000 umol/m^2 sec）的二氧化碳浓度（conc）。另外， Type是组间因子， conc是组内因子。 Type已经被存储为一个因子变量，但你还需要先将conc转换为因子变量。

CO2$conc <- factor(CO2$conc)
w1b1 <- subset(CO2, Treatment=='chilled’) #取子集，treatment为chilled的全都取出来
fit <- aov(uptake ~ (conc*Type) + Error(Plant/(conc)), w1b1)
summary(fit)

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="Interaction Plot for Plant Type and Concentration"))
boxplot(uptake ~ Type*conc, data=w1b1, col=(c("gold","green")),main="Chilled Quebec and Mississippi Plants", ylab="Carbon dioxide uptake rate (umol/m^2 sec)")
par(opar)

2.4 案例2

Example8_12  <- read.table ("example8_12.csv", header=TRUE, sep=",")
attach(Example8_12)
type  <-factor(type, order=FALSE)
time  <-factor(time, order=FALSE)
subject  <-factor(subject, order=FALSE)
fit <- aov(rate ~type*time +Error(subject/time))
summary(fit)
detach(Example8_12)

2.5 案例3

Example8_13  <- read.table ("example8_13.csv", header=TRUE, sep=",")
attach(Example8_13)
a  <-factor(a, order=FALSE)
b  <-factor(b, order=FALSE)
s  <-factor(s, order=FALSE)
time  <-factor(time, order=FALSE)
fit <- aov(y ~a*b*time +Error(s/time))
summary(fit)
detach(Example8_13)

四、多元方差分析

将某肝炎病人随机分成两组，一组施以新的疗法，另一组仍用传统的治疗方法，考察用二种方法治疗后对反映病人肝功能指标（如SGPT、 AST 、 ALT、 HCV等）的影响。

1、基本思想

与一个反应变量的方差分析相似，都是将反应变量的变异分解成为两部分：一部分为两组间变异（组别因素的效应），一部分为组内变异（随机误差）。然后对这两部分变异进行进行比较，看是否组间变异大于组内变异。不同的是，后者都是对组间均方与组内均方进行比较，而前者是对组间方差协方差矩阵与组内方差协方差矩阵进行比较

2、多元方差分析条件

各因变量服从多元正态分布：只要一个反应变量不服从正态分布，则这几个反应变量的联合分布肯定不服从多元正态分布。
各观察对象之间相互独立。
各组观察对象反应变量的方差-协方差矩阵相等。
反应变量间的确存在一定的关系，这可以从专业或研究目的角度予以判断。

3、举例说明

当因变量（结果变量）不止一个时，可用多元方差分析（MANOVA）对它们同时进行分析。

以MASS包中的UScereal数据集为例（Venables， Ripley（1999）），我们将研究美国谷物中的卡路里、脂肪和糖含量是否会因为储存架位置的不同而发生变化；其中1代表底层货架， 2代表中层货架， 3代表顶层货架。卡路里、脂肪和糖含量是因变量，货架是三水平（1、 2、 3）的自变量。分析过程见代码清单如下：

>library(MASS)
>attach(UScereal)
>shelf <- factor(shelf)# cbind() 函数将三个因变量（卡路里、脂肪和糖）合并成一个矩阵
>y <- cbind(calories, fat, sugars)
>aggregate(y, by=list(shelf), FUN=mean)#cov()则输出各谷物间的方差和协方差
>cov(y)#manova()函数能对组间差异进行多元检验
>fit <- manova(y ~ shelf)
>summary(fit)
>summary.aov(fit)

注意shelf变量已经转成了因子变量，因此它可以代表一个分组变量。由于多元检验是显著的，可以使用summary.aov()函数对每一个变量做单因素方差分析

4、假设检验

单因素多元方差分析有两个前提假设，一个是多元正态性，一个是方差协方差矩阵同质性。第一个假设即指因变量组合成的向量服从一个多元正态分布。可以用Q-Q图来检验该假设条件

>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, main="QQ Plot Assessing Multivariate Normality",ylab="Mahalanobis D2")
>abline(a=0,b=1) #斜率为1、截距项为0的直线上
>identify(coord$x, coord$y, labels=row.names(UScereal))

若有一个p×1的多元正态随机向量x，均值为μ，协方差矩阵为Σ，那么x与μ的马氏距离的平方服从自由度为p的卡方分布。 Q-Q图展示卡方分布的分位数，横纵坐标分别是样本量与马氏距离平方值。

如果点全部落在斜率为1、截距项为0的直线上，则表明数据服从多元正态分布。

这张图可以认为成正态分布

检验多元离群点**

若数据服从多元正态分布，那么点将落在直线上。还可以使用mvoutlier包中的ap.plot()函数来检验多元离群点。

library(mvoutlier)
outliers <- aq.plot(y)
outliers

5、稳健多元方差分析

如果多元正态性或者方差协方差均值假设都不满足或者你担心多元离群点，那么可以考虑用稳健或非参数版本的MANOVA 检验。稳健单因素MANOVA 可通过rrcov 包中的Wilks.test()函数实现。 vegan包中的adonis()函数则提供了非参数MANOVA的等同形式。

library(rrcov)
Wilks.test(y,shelf, method="mcd") #w是大写

五、用回归来做ANOVA

ANOVA和回归都是广义线性模型的特例。因此，本章所有的设计都可以用 lm()函数来分析。但是，为了更好地理解输出结果，需要弄明白在拟合模型时， R是如何处理类别型变量的。以单因素ANOVA问题为例，即比较五种降低胆固醇药物疗法（trt）的影响。

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)
contrasts(cholesterol$trt)fit.lm1 <- lm(response ~ trt, data=cholesterol, contrasts="contr.helmert") # WRONG
summary(fit.lm1)fit.lm2 <- lm(response ~ trt, data=cholesterol, contrasts="contr.SAS") # WRONG
summary(fit.lm2)fit.lm3 <- lm(response ~ trt, data=cholesterol, contrasts=list(trt="contr.SAS"))
summary(fit.lm3)fit.lm4 <- lm(response ~ trt, data=cholesterol, contrasts=list(trt="contr.sum"))
summary(fit.lm4)

六、Kruskal-Wallis检验和Friedman检验

如果无法满足ANOVA设计的假设，那么可以使用非参数方法来评估组间的差异。如果各组独立，则Kruskal-Wallis检验将是一种实用的方法。如果各组不独立（如重复测量设计或随机区组设计），那么Friedman检验会更合适。

1、Kruskal-Wallis

Kruskal-Wallis检验的调用格式为：

kruskal.test(y ~ A, data)

其中的y是一个数值型结果变量， A是一个拥有两个或更多水平的分组变量（若有两个水平，则它与Mann-Whitney U检验等价）。

2、Friedman检验

Friedman检验的调用格式为：

friedman.test(y ~ A | B, data)

其中的y是数值型结果变量， A是一个分组变量，而B是一个用以认定匹配观测的区组变量（blocking variable）。在以上两例中， data皆为可选参数，它指定了包含这些变量的矩阵或数据框。

3、案例1

考虑state.x77数据集。它包含了美国各州的人口、收入、文盲率、预期寿命、谋杀率和高中毕业率数据。如果你想比较美国四个地区（东北部、南部、中北部和西部）的文盲率，应该怎么做呢？

states <- data.frame(state.region, state.x77)
kruskal.test(Illiteracy ~ state.region, data=states)

4、案例2

example14_11  <- read.table ("example14_11.csv", header=TRUE, sep=",")
attach(example14_11)
kruskal.test(rate ~ group)
library(nparcomp)
nparcomp(rate ~ group, data=example14_11, alternative = "two.sided")
detach(example14_11)

5、案例3

example14_18  <- read.table ("example14_18.csv", header=TRUE, sep=",")
attach(example14_18)
friedman.test (rate~ treat|block)
library(PMCMR)
posthoc.friedman.nemenyi.test(rate,treat,block)
detach(example14_18)