方差分析（analysis of variation,简写为ANOVA）又称变异数分析或F检验,用于两个及两个以上样本均值差别的显著性检验,从函数的形式看，方差分析和回归都是广义线性模型的特例，回归分析lm()也能作方差分析。其目的是推断两组或多组数据的总体均值是否相同，检验两个或多个样本均值的差异是否有统计学意义。方差分析的基本思路为：将试验数据的总变异分解为来源于不同因素的相应变异，并作出数量估计，从而明确各个变异因素在总变异中所占的重要程度；也就是将试验数据的总变异方差分解成各变因方差，并以其中的误差方差作为和其他变因方差比较的标准，以推断其它变因所引起的变异量是否真实的一种统计分析方法。把对试验结果发生影响和起作用的自变量称为因素（factor），即我们所要检验的对象。如果方差分析研究的是一个因素对于试验结果的影响和作用，就称为单因素方差分析。因素的不同选择方案称之为因素的水平(level of factor)或处理(treatment)。因素的水平实际上就是因素的取值或者是因素的分组。样本数据之间差异如果是由于抽样的随机性造成的，称之为随机误差；如果是由于因素水平本身不同引起的差异，称之为系统误差。
方差分析的基本前提各组的观察数据，要能够被看作是从服从正态分布的总体中随机抽得的样本。各组的观察数据，是从具有相同方差的总体中抽取得到的。观察值是相互独立的。 方差分析的原假设：\(H_{0}\ \theta _{1}=\theta _{2}=...=\theta _{k}\) 即因素的不同水平对实验结果没有显著差异或影响。备择假设：不是所有的\(\theta _{i}\)都相等\((i=1,2,...,k)\)，即因素的不同水平对实验结果有显著差异或影响。
aov()函数的语法为aov(formula,data=dataframe),formula可使用的特殊符号如下，其中y为因变量，A、B、C为自变量。
| ～ |
分隔符，左边为因变量，右边为自变量。例y～A+B+C |
| + |
分隔自变量 |
| ： |
表示交互项，如y～A+B+A:B |
| * |
表示所有可能的交互项，如y～A * B *C等价于y～A+B+C+A:B+A:C+B:C+A:B:C |
| ^ |
表示交互项达到的某个次数，如y～(A+B+C)^2等价于y～A+B+C+A:B+A:C+B:C |
| . |
表示包含除因变量以外的所有变量。如y～. |
常用的方差设计表达式如下，其中小写字母表示定量变量，大写字母表示组别因子，Subject是被试着的标识变量。
| 单因素ANOVA |
Y~A |
| 含但个协变量的单因素ANCOVA |
Y~x+A |
| 双因素ANOVA |
Y~A*B |
| 含两个协变量的双因素ANCOVA |
Y~x1+x2+A*B |
| 随机化区组 |
y~B+A(B是区组因子) |
| 单因素组内ANOVA |
y~A+Error(Subject/A) |
| 含单个组内因子(w)和单个组间因子(b)的重复测量ANOVA |
Y~B*W+Error(Subject/W) |
组别间观测数相等的设计均衡设计(balanced design)，观测数不等的设计为非均衡设计(unbalanced design)。如果因子不止一个，且别是非平衡设计，或者存在协变量，表达式中的顺序会对结果造成影响。样本大小越不平衡，效应项的顺序对结果影响越大。通常，越基础的效应需要风在表达式的前面，如，先协变量，然后主效应，接着双因素的交互项，再接着是三因素的交互项。标准的anova()默认类型为序贯型，car包中的Anova()函数提供使用分层型和边界型(SAS和SPSS默认类型)的选项。
单因素方差分析（one-way ANOVA）
单因素方差分析是指对单因素试验结果进行分析，检验因素对试验结果有无显著性影响的方法。单因素方差分析是用来检验多个平均数之间的差异，从而确定因素对试验结果有无显著性影响的一种统计方法。对于完全随机设计试验且处理数大于2时可以用单因素方差分析（等于2 时用t检验）。离差平方和的分解公式为：SST(总和)=SSR(组间)+SSE(组内)，F统计量为MSR/MSE，MSR=SSR/k-1,MSE=SSE/n-k。其中SST为总离差、SSR为组间平方和、SSE为组内平方和或残差平方和、MSR为组间均方差、MSE为组内均方差。
例 某医院欲研究A、B、C三种降血脂药物对家兔血清肾素血管紧张素转化酶（ACE）的影响，将家兔随机分为三组，均喂以高脂饮食，分别给予不同的降血脂药物。一定时间后测定家兔血清ACE浓度（u/ml），A组（45 44 43 47 48 44 46 44 40 45 42 40 43 46 47 45 46 45 43 44），B组（45 48 47 43 46 47 48 46 43 49 46 43 47 46 47 46 45 46 44 45 46 44 43 42 45），c组（47 48 45 46 46 44 45 48 49 50 49 48 47 44 45 46 45 43 44 45 46 43 42），问三组家兔血清ACE浓度是否相同？
a <- c(45, 44, 43, 47, 48, 44, 46, 44, 40, 45, 42, 40, 43, 46, 47, 45,
46, 45, 43, 44)
b <- c(45, 48, 47, 43, 46, 47, 48, 46, 43, 49, 46, 43, 47, 46, 47, 46,
45, 46, 44, 45, 46, 44, 43, 42, 45)
c <- c(47, 48, 45, 46, 46, 44, 45, 48, 49, 50, 49, 48, 47, 44, 45, 46,
45, 43, 44, 45, 46, 43, 42)
dfCRp <- data.frame(length = c(a, b, c), site = factor(c(rep("1", 20),
rep("2", 25), rep("3", 23))))
boxplot(length ~ site, data = dfCRp, xlab = "Sites", ylab = "Length")
plot.design(length ~ site, fun = mean, data = dfCRp, main = "Group means")
箱形图中可观察到不同的因素对于因变量的影响。
假设检验
方差分析需要一定的假设，即数据集应该符合正态和各组的方差相等，可以分别用shapiro.test和bartlett.test检验从P值观察到这两个假设是符合的。对于不符合假设的情况，我们就要用到非参数方法，例如Kruskal-Wallis秩和检验
shapiro.test(dfCRp$length)
##
## Shapiro-Wilk normality test
##
## data: dfCRp$length
## W = 0.97397, p-value = 0.1654
bartlett.test(length ~ site,data = dfCRp) #Fligner-Killeen(fligner.test()函数)和Brown-Forsythe检验(HH包中的hov()函数)也可以用来检验方差齐性
##
## Bartlett test of homogeneity of variances
##
## data: length by site
## Bartlett's K-squared = 0.76406, df = 2, p-value = 0.6825
正态性检验和方差齐性检验P值均大于0.05，可以认为数据满足正态性和方差齐性的要求。
oneway.test()和aov()函数进行方差分析
aovCRp =aov(length ~ site, data = dfCRp)
summary(aovCRp)
## Df Sum Sq Mean Sq F value Pr(>F)
## site 2 26.29 13.146 3.244 0.0454 *
## Residuals 65 263.40 4.052
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#oneway.test(length ~ site, data=dfCRp, var.equal=TRUE)，与aov()结果基本相同。
plotmeans(length ~ site,data =dfCRp ) #绘制有置信区间的组均值图
par(mfrow=c(2,2))
plot(aovCRp)
par(mfrow=c(1,1))
用aov函数建立单因子方差模型，从结果的P值可看到各组均值有显著不同。Sum Sq = deviance (within groups, and residual)，总方差和（分别有groups和residual的），Mean Sq = variance (within groups, and residual)，平均方差和（分别有groups和residual的）。单因子方差分析结果显示F value = 3.24 ，Pr(>F) = 0.045，因此拒绝原假设，即认为三组组家兔血清ACE浓度在统计学上有显著差异。
模型比较
(anovaCRp <- anova(lm(length ~ site, data=dfCRp)))
## Analysis of Variance Table
##
## Response: length
## Df Sum Sq Mean Sq F value Pr(>F)
## site 2 26.292 13.1462 3.2442 0.0454 *
## Residuals 65 263.399 4.0523
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(length ~ 1, data=dfCRp), lm(length ~ site, data=dfCRp))
## Analysis of Variance Table
##
## Model 1: length ~ 1
## Model 2: length ~ site
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 67 289.69
## 2 65 263.40 2 26.293 3.2442 0.0454 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anovaCRp["Residuals", "Sum Sq"]
## [1] 263.3987
比较不含有自变量和含有一个自变量site模型，含有site变量的模型结果较好(残差的总方差和较小)。
效果大小(Effect size)
效果大小是指某个特定总体中的某种特殊的非零的数值,这个数值越大,就表明由研究者所处理的研究现象所造成的效果越大。效果大小本身可以被视为是一种参数:当原假设为真时,效果大小的值为零;当原假设为假时,效果大小为某种非零的值。因此,可以把效果大小视为某种与原假设分离程度的指标 。方差分析效果大小的含义也基本上与Z检验或t检验的效果大小的含义相同只不过它反映的是多组实验处理下不同组之间实验效果差异大小的指标。常用的指标如下\(\eta^{2}=\frac{SS}{SS_{total}}\),\(f=\sqrt{\frac{\eta^{2}}{1-\eta^{2}}}\),\(\omega ^{2}=\frac{SS-DF*MSE}{SS_{total}+MSE}\)。
dfSSb <- anovaCRp["site", "Df"]
SSb <- anovaCRp["site", "Sum Sq"]
MSb <- anovaCRp["site", "Mean Sq"]
SSw <- anovaCRp["Residuals", "Sum Sq"]
MSw <- anovaCRp["Residuals", "Mean Sq"]
(etaSq <- SSb / (SSb + SSw))# DescTools包中EtaSq(aovCRp, type=1)函数可以计算
## [1] 0.09076038
(omegaSq <- dfSSb * (MSb-MSw) / (SSb + SSw + MSw))
## [1] 0.06191765
(f <- sqrt(etaSq / (1-etaSq)))
## [1] 0.3159432
\(\eta^{2},\omega^{2},f^{2}\)值如上，如\(\eta^{2}\)实验处理之后各组间平方和在总体平方和中所占的比重，\(\eta^{2}\)越大反映实验效果大。一般\(\eta^{2}\)大于0.14，就认为效果有大的效果。
多重比较
方差分析只告诉我们这三组之间是不同的，但没有告诉哪两组之间有明显差别，此时需要使用TukeyHSD函数进行均值的多重比较分析，从结果中观察到有一个两两比较是不显著的。 ####计划好的多重比较(Planned comparisons - a-priori) 在收集数据之前就已确定。它与实验目的有关，反映了实验者的意图。可以直接进行计划好的多重比较，不用考虑基本的“均值相等的 F-test”。
cntrMat <- rbind("a-c" =c(1,0,-1),
"1/3*(a+b)-c"=c(1/3,1/3,-1),
"b-c" =c(0,1,-1))
summary(glht(aovCRp, linfct=mcp(site=cntrMat), alternative="less"),
test=adjusted("none"))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: User-defined Contrasts
##
##
## Fit: aov(formula = length ~ site, data = dfCRp)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(<t)
## a-c >= 0 -1.5196 0.6155 -2.469 0.00809 **
## 1/3*(a+b)-c >= 0 -1.1429 0.5331 -2.144 0.01790 *
## b-c >= 0 -0.3896 0.5816 -0.670 0.25268
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- none method)
#pairwise.t.test(dfCRp$length, dfCRp$site, p.adjust.method="bonferroni") #结果与glht()函数类似。
依据事先实验的目的，进行多重比较，a组和c组，a、b组和c组的差异有显著意义。
非计划的多重比较(Planned comparisons - post-hoc)
在查看数据之后，并且“均值相等的 F-test” 结果显著情况下才进行。它用于探究研究者感兴趣但头脑中没有特定假设。
#ScheffeTest检验
ScheffeTest(aovCRp, which="site", contrasts=t(cntrMat)) #DescTools包
##
## Posthoc multiple comparisons of means : Scheffe Test
## 95% family-wise confidence level
##
## $site
## diff lwr.ci upr.ci pval
## 1-3 -1.5195652 -3.061467 0.02233664 0.0543 .
## 1,2-3 -15.9262319 -17.092478 -14.75998618 <2e-16 ***
## 2-3 -0.3895652 -1.846661 1.06753062 0.7997
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Tukey HSD检验
(tHSD <- TukeyHSD(aovCRp))
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = length ~ site, data = dfCRp)
##
## $site
## diff lwr upr p adj
## 2-1 1.1300000 -0.31850529 2.578505 0.1552673
## 3-1 1.5195652 0.04333482 2.995796 0.0422495
## 3-2 0.3895652 -1.00547115 1.784602 0.7817904
plot(tHSD)
置信区间包含0说明差异不显著。
multcomp包中glht()函数提供了多重均值更全面的方法，适用于线性模型和广义线性模型。下面的代码重现Tukey HSD检验。
tukey <- glht(aovCRp, linfct=mcp(site="Tukey"))
summary(tukey)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = length ~ site, data = dfCRp)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 2 - 1 == 0 1.1300 0.6039 1.871 0.1552
## 3 - 1 == 0 1.5196 0.6155 2.469 0.0422 *
## 3 - 2 == 0 0.3896 0.5816 0.670 0.7817
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
plot(cld(tukey,level = .05),col="lightgrey") #cld()函数中level选项设置了使用显著水平0.05，即95%的置信区间
有相同字母的组（箱线图表示）说明均值差异不显著。
离群点检测
outlierTest(aovCRp)
##
## No Studentized residuals with Bonferonni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferonni p
## 9 -2.288156 0.025441 NA
离群点检测结果显示，数据中没有离群点（当p>1时产生NA）。
残差的相关检验
残差的正态性检验
Estud <- rstudent(aovCRp)
shapiro.test(Estud)
##
## Shapiro-Wilk normality test
##
## data: Estud
## W = 0.9876, p-value = 0.7381
qqnorm(Estud, pch=20, cex=2)
qqline(Estud, col="gray60", lwd=2)
残差满足正态性的要求。
残差的方差齐性检验，levene检验是对方差模型的残差进行组间齐性检验的，bartlett.test是对原始数据进行检验。
plot(Estud ~ dfCRp$site, main="Residuals per group")
leveneTest(aovCRp)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.3886 0.6795
## 65
对模型的残差进行组间方差齐性检验，P值大于0.05满足残差方差齐性的要求。
单因素协方差分析(Analysis of covariance ,ANCOVA)
单因素协方差分析在单因素方差分析的基础上包含一个或多个定量的协变量。
例 multcomp包中litter数据集是怀孕小白鼠被分为四个小组，每个小组接受不同剂量（0、5、50和500）的药物处理dose为自变量，产下幼崽的体重weigth均值为因变量，怀孕时间gesttime为协变量。
data(litter,package = "multcomp")
pander(head(litter))
| 0 |
28.05 |
22.5 |
15 |
| 0 |
33.33 |
22.5 |
14 |
| 0 |
36.37 |
22 |
14 |
| 0 |
35.52 |
22 |
13 |
| 0 |
36.77 |
21.5 |
15 |
| 0 |
29.6 |
23 |
5 |
ddply(.data = litter,.(dose),summarize,mean=mean(weight))
## dose mean
## 1 0 32.30850
## 2 5 29.30842
## 3 50 29.86611
## 4 500 29.64647
单因素协方差分析
shapiro.test(litter$weight)
##
## Shapiro-Wilk normality test
##
## data: litter$weight
## W = 0.9674, p-value = 0.05282
bartlett.test(weight~dose,data = litter)
##
## Bartlett test of homogeneity of variances
##
## data: weight by dose
## Bartlett's K-squared = 9.6497, df = 3, p-value = 0.02179
ancova <- aov(weight~gesttime+dose,data = litter)
summary(ancova)
## Df Sum Sq Mean Sq F value Pr(>F)
## gesttime 1 134.3 134.30 8.049 0.00597 **
## dose 3 137.1 45.71 2.739 0.04988 *
## Residuals 69 1151.3 16.69
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
数据满足正态性的要求,但不满足方差齐性的要求。ANCOVA检验结果表明怀孕时间gesttime与出生体重weight相关，在控制怀孕时间后，每种药物剂量dose下出生体重weight均值不同。
调整的组均值
去除协变量效用的组均值，可以使用effects包中的effect()函数计算。
effect("dose",ancova)
##
## dose effect
## dose
## 0 5 50 500
## 32.35367 28.87672 30.56614 29.33460
多重比较
contrast <- rbind("no drug vs drug"=c(3,-1,-1,-1)) #设定第一组和其他三组的均值进行比较
summary(glht(ancova,linfct=mcp(dose=contrast)))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: User-defined Contrasts
##
##
## Fit: aov(formula = weight ~ gesttime + dose, data = litter)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## no drug vs drug == 0 8.284 3.209 2.581 0.012 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
在未用药和用药条件下，出生体重有显著的不同。
检验回归斜率的同质性
ANCOVA模型假定回归斜率相同，如果ANCOVA模型包含交互项，则需要对回归斜率的同质性进行检验。本例中假定四个处理组通过怀孕时间来预测出生体重的回归斜率都相同。
summary(aov(weight~gesttime*dose,data = litter))
## Df Sum Sq Mean Sq F value Pr(>F)
## gesttime 1 134.3 134.30 8.289 0.00537 **
## dose 3 137.1 45.71 2.821 0.04556 *
## gesttime:dose 3 81.9 27.29 1.684 0.17889
## Residuals 66 1069.4 16.20
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
交互效应不显著，支持斜率相等的假设。如果交互效应显著，则意味怀孕时间和出生体重的关系依赖于药物剂量，需使用不需要假设回归斜率同质性的非参数ANCOVA方法，如sm包中的sm.ancova()函数。
结果可视化
library(HH)
## Loading required package: lattice
## Loading required package: grid
## Loading required package: latticeExtra
## Loading required package: RColorBrewer
## Loading required package: gridExtra
##
## Attaching package: 'HH'
## The following object is masked from 'package:gplots':
##
## residplot
## The following object is masked from 'package:DescTools':
##
## OddsRatio
## The following objects are masked from 'package:car':
##
## logit, vif
ancova(weight~gesttime+dose,data = litter)
## Analysis of Variance Table
##
## Response: weight
## Df Sum Sq Mean Sq F value Pr(>F)
## gesttime 1 134.30 134.304 8.0493 0.005971 **
## dose 3 137.12 45.708 2.7394 0.049883 *
## Residuals 69 1151.27 16.685
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
用怀孕时间预测出生体重的回归线相互平行，只是截距不同。随着怀孕时间的增加，出生体重也会增加。若用ancova(weight~gesttime*dose,data = litter)生成的图形将允许斜率和截距依据组别发生变化，对违背回归斜率同质性的实例比较有用。
I类型的平方和(Type I sum of squares)单因素协方差分析
I类型的平方和效应根据表达式中先出现的效应做调整。A不做调整,B根据A调整,A:B交互项根据A和B调整。
fitFull <- lm(weight~gesttime+dose,data = litter)
fitGrp <- lm(weight ~ dose, data=litter)
fitRegr <- lm(weight ~ gesttime, data=litter)
anova(fitFull)
## Analysis of Variance Table
##
## Response: weight
## Df Sum Sq Mean Sq F value Pr(>F)
## gesttime 1 134.30 134.304 8.0493 0.005971 **
## dose 3 137.12 45.708 2.7394 0.049883 *
## Residuals 69 1151.27 16.685
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
II/III类型的平方和(Type II/III sum of squares)单因素协方差分析
II类型的平方和效应根据同水平或低水平的效应做调整。A根据B调整,B依据A调整,A:B交互项同时根 据A和B调整。III类型的平方和每个效应根据模型其他各效应做相应调整。A根据B和A:B做调整,A:B交互项根据A和B 调整。
fitFiii <- lm(weight~gesttime+dose,contrasts=list(dose=contr.sum), data=litter)
Anova(fitFiii, type="III")
## Anova Table (Type III tests)
##
## Response: weight
## Sum Sq Df F value Pr(>F)
## (Intercept) 60.14 1 3.6045 0.06181 .
## gesttime 161.49 1 9.6789 0.00271 **
## dose 137.12 3 2.7394 0.04988 *
## Residuals 1151.27 69
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
基于II类型的平方和的模型比较
anova(fitRegr, fitFull)
## Analysis of Variance Table
##
## Model 1: weight ~ gesttime
## Model 2: weight ~ gesttime + dose
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 72 1288.4
## 2 69 1151.3 3 137.12 2.7394 0.04988 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(fitGrp, fitFull)
## Analysis of Variance Table
##
## Model 1: weight ~ dose
## Model 2: weight ~ gesttime + dose
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 70 1312.8
## 2 69 1151.3 1 161.49 9.6789 0.00271 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RSS值较小的模型较好，以怀孕时间gesttime为协变量的单因素模型比仅含有药物剂量dose和怀孕时间gesttime的模型校好。
回归系数(Test individual regression coefficients)
(sumRes <- summary(fitFull))
##
## Call:
## lm(formula = weight ~ gesttime + dose, data = litter)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.5649 -2.0068 0.1476 3.0755 7.3023
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -45.366 24.984 -1.816 0.07374 .
## gesttime 3.519 1.131 3.111 0.00271 **
## dose5 -3.477 1.318 -2.639 0.01027 *
## dose50 -1.788 1.344 -1.330 0.18780
## dose500 -3.019 1.352 -2.232 0.02883 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.085 on 69 degrees of freedom
## Multiple R-squared: 0.1908, Adjusted R-squared: 0.1439
## F-statistic: 4.067 on 4 and 69 DF, p-value: 0.005105
confint(fitFull) #95%置信区间
## 2.5 % 97.5 %
## (Intercept) -95.206425 4.4752052
## gesttime 1.262361 5.7749312
## dose5 -6.105369 -0.8485267
## dose50 -4.468125 0.8930650
## dose500 -5.716971 -0.3211661
效果大小(Effect size)
\(\omega^{2}\)基于II类型的平方和
anRes <- anova(fitRegr, fitFull)
dfGrp <- anRes[2, "Df"]
dfE <- anRes[2, "Res.Df"]
MSgrp <- anRes[2, "Sum of Sq"] / dfGrp
MSE <- anRes[2, "RSS"] / dfE
SST <- sum(anova(fitFull)[ , "Sum Sq"])
(omegaSqHat <- dfGrp*(MSgrp - MSE) / (SST + MSE))
## [1] 0.0604895
效应\(\omega^{2}\)值如上。
调整的组均值
aovAncova <- aov(weight~gesttime+dose,data = litter)
YMjAdj <- effect("dose", aovAncova)
summary(YMjAdj)
##
## dose effect
## dose
## 0 5 50 500
## 32.35367 28.87672 30.56614 29.33460
##
## Lower 95 Percent Confidence Limits
## dose
## 0 5 50 500
## 30.53131 26.98687 28.59369 27.34813
##
## Upper 95 Percent Confidence Limits
## dose
## 0 5 50 500
## 34.17604 30.76658 32.53860 31.32108
双因素方差分析（Two-way ANOVA）
研究两个因素的不同水平对试验结果的影响是否显著的问题就称作双因素方差分析，分别对两个因素进行检验，考察各自的作用，同时分析两个因素（因素A和因素 B）对试验结果的影响。如果因素A和因素B对试验结果的影响是相互独立的，则可以分别考察各自的影响，这种双因素方差分析称为无交互作用的双因素方差分析，也叫无重复双因素方差分析。无交互作用的双因素方差分析，相当于对每个因素分别进行单因素方差分析。如果因素A和因素B除了各自对试验结果的影响外，还产生额外的新影响，这种额外的影响称为交互作用，这时的双因素方差分析则称为有交互作用的双因素方差分析，也叫有重复双因素方差分析。可用于随机区组实验设计，用来分析两个因素的不同水平对结果是否有显著影响，以及两因素之间是否存在交互效应。
例 基础安装中的ToothGrowth数据集是随机分配60只豚鼠,分别采用两种喂食方法supp(橙汁或维生素C),各喂食方法中抗坏血酸含量有三种水平dose(0.5mg、1mg或2mg),每种处理方式组合都被分配10只豚鼠,牙齿长度len为因变量。
pander(head(ToothGrowth))
| 4.2 |
VC |
0.5 |
| 11.5 |
VC |
0.5 |
| 7.3 |
VC |
0.5 |
| 5.8 |
VC |
0.5 |
| 6.4 |
VC |
0.5 |
| 10 |
VC |
0.5 |
attach(ToothGrowth)
table(supp,dose)
## dose
## supp 0.5 1 2
## OJ 10 10 10
## VC 10 10 10
ddply(.data = ToothGrowth,.(supp,dose),summarise,mean=mean(len))
## supp dose mean
## 1 OJ 0.5 13.23
## 2 OJ 1.0 22.70
## 3 OJ 2.0 26.06
## 4 VC 0.5 7.98
## 5 VC 1.0 16.77
## 6 VC 2.0 26.14
ddply(.data = ToothGrowth,.(supp,dose),summarise,sd=sd(len))
## supp dose sd
## 1 OJ 0.5 4.459709
## 2 OJ 1.0 3.910953
## 3 OJ 2.0 2.655058
## 4 VC 0.5 2.746634
## 5 VC 1.0 2.515309
## 6 VC 2.0 4.797731
table语句的预处理表明该设计是均衡设计(各设计单元中样本大小都相同),ddply语句处理可获得各单元的均值和标准差。
I型双因素方差分析(SS type I)
aovCRFpq <- aov(len~ supp*dose, data=ToothGrowth)
summary(aovCRFpq)
## Df Sum Sq Mean Sq F value Pr(>F)
## supp 1 205.4 205.4 12.317 0.000894 ***
## dose 1 2224.3 2224.3 133.415 < 2e-16 ***
## supp:dose 1 88.9 88.9 5.333 0.024631 *
## Residuals 56 933.6 16.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
par(mfrow=c(2,2))
plot(aovCRFpq)
par(mfrow=c(1,1))
得到方差分析表,可以看到主效应(supp和dose)和交互效应都非常显著。
II/III型双因素方差分析(SS type II or III)
ToothGrowth$supp <- as.factor(ToothGrowth$supp) #转为因子
ToothGrowth$dose <- as.factor(ToothGrowth$dose)
fitIII <- lm(len ~ supp + dose + supp:dose, data=ToothGrowth,contrasts=list(supp=contr.sum, dose=contr.sum))
Anova(fitIII, type="III")
## Anova Table (Type III tests)
##
## Response: len
## Sum Sq Df F value Pr(>F)
## (Intercept) 21236.5 1 1610.393 < 2.2e-16 ***
## supp 205.4 1 15.572 0.0002312 ***
## dose 2426.4 2 92.000 < 2.2e-16 ***
## supp:dose 108.3 2 4.107 0.0218603 *
## Residuals 712.1 54
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
得到方差分析表,可以看到主效应(supp和dose)和交互效应都非常显著。
绘制边际均数及格均数图
plot.design(len ~ supp*dose, data=ToothGrowth, main="Marginal means")
interaction.plot(ToothGrowth$dose, ToothGrowth$supp,ToothGrowth$len,main="Cell means", col=c("red", "blue", "green"), lwd=2) #interaction.plot(f1, f2, y)展示双因素方差分析的交互效应，如果f1和f2是因子,作y的均值图,以f1的不同值作为x轴,而f2的不同值对应不同曲线;可以用选项fun指定y的其他的统计量(缺省计算均值,fun=mean)
gplots包中的plotmeans()函数来展示交互效应，HH包中interaction2wt()函数来展示交互效应
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) #interaction2wt() 函数,因为它能展示任意复杂度设计(双因素方差分析、三因素方差分析等)的主效应(箱线图)和交互效应
效果大小(Effect size estimate)
EtaSq(aovCRFpq, type=1)
## eta.sq eta.sq.part
## supp 0.05948365 0.18029211
## dose 0.64431327 0.70435310
## supp:dose 0.02575745 0.08695875
\(\eta_{p}^{2}\)值如上。
简单效应(Simple effects)
简单效应指一个因素的不同水平在另一个因素的某个水平上的变异。phia包中testInteractions()可计算简单效应。
testInteractions(aovCRFpq, fixed="dose", across="supp", adjustment="none")
## Warning in testInteractions(aovCRFpq, fixed = "dose", across = "supp",
## adjustment = "none"): Some factors with specified contrasts are not in the
## model and will be ignored.
## F Test:
## P-value adjustment method: none
## Value Df Sum of Sq F Pr(>F)
## Mean 3.7 1 205.35 12.317 0.0008936 ***
## Residuals 56 933.63
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#testInteractions(aovCRFpq, fixed="supp", across="dose", adjustment="none")
喂食方法supp在剂量dose的0.5mg和1mg水平上的变异显著。同样的交换dose变量和supp变量后，剂量dose在喂食方法OJ和VC的水平上变异显著。
多重比较
计划好的主效应(Main effects)多重比较,计划比较(0.5mg、1mg)和2mg剂量，0.5mg和2mg剂量之间是否有差别
cMat <- rbind("c1"=c(1/2, 1/2,-1),"c2"=c(-1,0,1))
aovCRFpq <- aov(len~ supp*dose, data=ToothGrowth)
summary(glht(aovCRFpq), linfct=mcp(dose=cMat), alternative="two.sided",test=adjusted("bonferroni"))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: aov(formula = len ~ supp * dose, data = ToothGrowth)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) == 0 13.230 1.148 11.521 2.66e-15 ***
## suppVC == 0 -5.250 1.624 -3.233 0.0126 *
## dose1 == 0 9.470 1.624 5.831 1.91e-06 ***
## dose2 == 0 12.830 1.624 7.900 8.58e-10 ***
## suppVC:dose1 == 0 -0.680 2.297 -0.296 1.0000
## suppVC:dose2 == 0 5.330 2.297 2.321 0.1446
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)
结果显示0.5mg、1mg和2mg剂量，0.5mg和2mg剂量之间的差异显著。
非计划的多重比较
非计划好的主效应多重比较,dose变量的两两比较。
aovCRF <- aov(len~ supp*dose, data=ToothGrowth)
TukeyHSD(aovCRF, which="dose")
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = len ~ supp * dose, data = ToothGrowth)
##
## $dose
## diff lwr upr p adj
## 1-0.5 9.130 6.362488 11.897512 0.0e+00
## 2-0.5 15.495 12.727488 18.262512 0.0e+00
## 2-1 6.365 3.597488 9.132512 2.7e-06
multcomp包中的glht()函数
tukey <- glht(aovCRF, linfct=mcp(dose="Tukey"))
## Warning in mcp2matrix(model, linfct = linfct): covariate interactions found
## -- default contrast might be inappropriate
summary(tukey)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = len ~ supp * dose, data = ToothGrowth)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1 - 0.5 == 0 9.470 1.624 5.831 <1e-04 ***
## 2 - 0.5 == 0 12.830 1.624 7.900 <1e-04 ***
## 2 - 1 == 0 3.360 1.624 2.069 0.106
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
confint(tukey) #95%置信区间
##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = len ~ supp * dose, data = ToothGrowth)
##
## Quantile = 2.4108
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 1 - 0.5 == 0 9.4700 5.5549 13.3851
## 2 - 0.5 == 0 12.8300 8.9149 16.7451
## 2 - 1 == 0 3.3600 -0.5551 7.2751
结果显示0.5mg和1mg剂量，0.5mg和2mg剂量之间的差异显著。
单元多重比较(Cell comparisons using the associated one-way ANOVA)
ToothGrowth$comb <- interaction(ToothGrowth$dose, ToothGrowth$supp)
aovCRFpqA <- aov(len ~ comb, data=ToothGrowth)
cntrMat <- rbind("c1"=c(-1/2, 1/4, -1/2, 1/4, 1/4, 1/4),
"c2"=c( 0, 0, -1, 0, 1, 0),
"c3"=c(-1/2, -1/2, 1/4, 1/4, 1/4, 1/4))
summary(glht(aovCRFpqA, linfct=mcp(comb=cntrMat), alternative="greater"),
test=adjusted("none"))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: User-defined Contrasts
##
##
## Fit: aov(formula = len ~ comb, data = ToothGrowth)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>t)
## c1 <= 0 -1.2475 0.9945 -1.254 0.892
## c2 <= 0 -9.2900 1.6240 -5.720 1.000
## c3 <= 0 1.2725 0.9945 1.280 0.103
## (Adjusted p values reported -- none method)
计划的单元的多重比较中，未发现显著的差异。
非计划Scheffe检验
ScheffeTest(aovCRFpqA, which="comb", contrasts=t(cntrMat)) #Post-hoc Scheffe tests using the associated one-way ANOVA
##
## Posthoc multiple comparisons of means : Scheffe Test
## 95% family-wise confidence level
##
## $comb
## diff lwr.ci upr.ci pval
## 1.OJ,0.5.VC,1.VC,2.VC-0.5.OJ,2.OJ -1.2475 -4.682547 2.187547 0.9020
## 1.VC-2.OJ -9.2900 -14.899408 -3.680592 8e-05 ***
## 2.OJ,0.5.VC,1.VC,2.VC-0.5.OJ,1.OJ 1.2725 -2.162547 4.707547 0.8943
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ScheffeTest(aovCRFpq, which="dose", contrasts=c(-1, 1/2, 1/2)) #Post-hoc Scheffe tests for marginal means
##
## Posthoc multiple comparisons of means : Scheffe Test
## 95% family-wise confidence level
##
## $dose
## diff lwr.ci upr.ci pval
## 1,2-0.5 12.3125 8.877453 15.74755 1.2e-14 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
残差的相关检验
正态性检验
Estud <- rstudent(aovCRFpq)
qqnorm(Estud, pch=20, cex=2)
qqline(Estud, col="gray60", lwd=2)
shapiro.test(Estud)
##
## Shapiro-Wilk normality test
##
## data: Estud
## W = 0.98457, p-value = 0.6478
P值大于0.05，可以认为残差满足正态性。
plot(Estud ~ ToothGrowth$comb, main="Residuals per group")
leveneTest(aovCRFpq)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 1.7086 0.1484
## 54
P值大于0.05，可以认为残差满足方差齐性。
重复测量方差分析
重复测量数据的方差分析是对同一因变量进行重复测量的一种试验设计技术。在给予一种或多种处理后，分别在不同的时间点上通过重复测量同一个受试对象获得的指标的观察值，或者是通过重复测量同一个个体的不同部位（或组织）获得的指标的观察值。重复测量数据在科学研究中十分常见，常用来分析该观察指标在不同时间点上的变化特点。分析前要对重复测量数据之间是否存在相关性进行球形检验。如果该检验结果为P﹥0.05，则说明重复测量数据之间不存在相关性，测量数据符合Huynh-Feldt条件，可以用单因素方差分析的方法来处理；如果检验结果P﹤0.05，则说明重复测量数据之间是存在相关性的，所以不能用单因素方差分析的方法处理数据。球形条件不满足时常有两种方法可供选择：（1）采用MANOVA（多变量方差分析方法）；（2）对重复测量ANOVA检验结果中与时间有关的F值的自由度进行调整。
在重复测量的方差分析中，实验对象被测量多次，所以会存在组内因子，组内因子要以下面的形式特别标明出来，其中B是组间因子，W是组内因子，subject是实验对象的ID， model=aov(Y ~ B * W + Error(Subject/W)) 上述方法的前提是对应组内因子不同水平的数据是等方差的，当传统方法的假设得不到满足时，则应用lme4包中lmer函数，利用混合效应模型来解决问题。
单因素重复测量方差分析(One-way repeated measures ANOVA)
单因素重复测量方差分析通常只有组内因素，没有组间因素。 例 将42名诊断为胎粪吸入综合症的新生儿患儿随机分为肺表面活性物质治疗组（PS组）和常规治疗组（对照组），每组各21例。PS组和对照组两组所有患儿均给予除用药外的其他相应的对症治疗。PS组患儿给予牛肺表面活性剂70mg/kg治疗。采集PS组及对照组患儿0小时，治疗后24小时和72小时静脉血2ml，离心并提取上清液后保存备用并记录血清中VEGF的含量变化情况。在治疗组，不同时间的记录的VEGF是否有差异？
MAS <- read.csv("MAS.csv",header = T)
pander(head(MAS))
| 1 |
time0 |
contrast |
1.03 |
| 2 |
time0 |
contrast |
1.772 |
| 3 |
time0 |
contrast |
0.094 |
| 4 |
time0 |
contrast |
0.596 |
| 5 |
time0 |
contrast |
1.314 |
| 6 |
time0 |
contrast |
1.516 |
传统的重复测量方差分析(Traditional univariate approach)
aov()函数在处理重复测量设计时,需要有长格式(long format)数据才能拟合模型，在长格式中,因变量的每次测量都要放到它独有的行中。reshape包可方便地将数据转换为相应的格式
dfRBpL <- subset(MAS,treatment=="ps")
dfRBpL$id <- as.factor(dfRBpL$id)
aovRBp <- aov(value ~ time + Error(id/time), data=dfRBpL) #id和time需为因子
summary(aovRBp)
##
## Error: id
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 20 0.6525 0.03263
##
## Error: id:time
## Df Sum Sq Mean Sq F value Pr(>F)
## time 2 0.8262 0.4131 17.59 3.31e-06 ***
## Residuals 40 0.9397 0.0235
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
在治疗组，不同时间记录的VEGF差异显著。
效果大小(Effect size estimate)
EtaSq(aovRBp, type=1)
## eta.sq eta.sq.part eta.sq.gen
## time 0.3416423 0.46788 0.3416423
\(\eta_{g}^{2}\)值如上。
重复测量的多变量方差分析(Multivariate approach)
\(Hotelling's T^{2}\)检验是单变量检验的推广，常用于两组均向量的比较。
DVw<- data.matrix(subset(dfRBpW,select=c("value.time0", "value.time24", "value.time72")))
diffMat <- combn(1:length(levels(dfRBpL$time)), 2, function(x) {DVw[ , x[1]] - DVw[ , x[2]]})
DVdiff<- diffMat[ , 1:(length(levels(dfRBpL$time))-1), drop=FALSE]
muH0 <- rep(0, ncol(DVdiff))
HotellingsT2Test(DVdiff, mu=muH0)
##
## Hotelling's one sample T2-test
##
## data: DVdiff
## T.2 = 14.925, df1 = 2, df2 = 19, p-value = 0.000127
## alternative hypothesis: true location is not equal to c(0,0)
P值小于0.05，可以认为不同时间记录的VEGF差异显著。
car包Anova()函数进行多变量方差分析
多元方差分析(multivariate analysis of variance, MANOVA)是单变量方差分析和\(Hotelling's T^{2}\)检验的推广，用于多组均向量间的比较。
summary(AnovaRBp, multivariate=TRUE, univariate=FALSE)
##
## Type III Repeated Measures MANOVA Tests:
##
## ------------------------------------------
##
## Term: (Intercept)
##
## Response transformation matrix:
## (Intercept)
## value.time0 1
## value.time24 1
## value.time72 1
##
## Sum of squares and products for the hypothesis:
## (Intercept)
## (Intercept) 241.8564
##
## Sum of squares and products for error:
## (Intercept)
## (Intercept) 1.957503
##
## Multivariate Tests: (Intercept)
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.99197 2471.071 1 20 < 2.22e-16 ***
## Wilks 1 0.00803 2471.071 1 20 < 2.22e-16 ***
## Hotelling-Lawley 1 123.55357 2471.071 1 20 < 2.22e-16 ***
## Roy 1 123.55357 2471.071 1 20 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------
##
## Term: time
##
## Response transformation matrix:
## time1 time2
## value.time0 1 0
## value.time24 0 1
## value.time72 -1 -1
##
## Sum of squares and products for the hypothesis:
## time1 time2
## time1 1.6046679 0.5625321
## time2 0.5625321 0.1972012
##
## Sum of squares and products for error:
## time1 time2
## time1 1.0566451 0.2309589
## time2 0.2309589 0.5838118
##
## Multivariate Tests: time
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.6110546 14.92502 2 19 0.00012704 ***
## Wilks 1 0.3889454 14.92502 2 19 0.00012704 ***
## Hotelling-Lawley 1 1.5710549 14.92502 2 19 0.00012704 ***
## Roy 1 1.5710549 14.92502 2 19 0.00012704 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
P值小于0.05，可以认为不同时间记录(time)的VEGF差异显著，截距(Intercept)差异显著，但无实际意义。
双因素重复测量方差分析(Two-way repeated-measures ANOVA)
双因素重复测定资料中的因素是指一个组间因素（处理因素）和一个组内因素（时间因素）。组间因素是指分组或分类变量，它把所有受试对象按分类变量的水平分为几个组。组内因素是指重复测定的时间变量。
例 将42名诊断为胎粪吸入综合症的新生儿患儿随机分为肺表面活性物质治疗组（PS组）和常规治疗组（对照组），每组各21例。PS组和对照组两组所有患儿均给予除用药外的其他相应的对症治疗。PS组患儿给予牛肺表面活性剂70mg/kg治疗。采集PS组及对照组患儿0小时，治疗后24小时和72小时静脉血2ml，离心并提取上清液后保存备用并记录血清中VEGF的含量变化情况。不同组间不同时间的记录的VEGF是否有差异？
传统的重复测量方差分析(Traditional univariate approach)
aov()同样需要长格式数据。
dfRBFpqL <- read.csv("MAS.csv",header = T)
id <- factor(rep(1:21, times=2*3)) #dfRBFpqL$id <- as.factor(dfRBFpqL$id)这种法方法每组因子水平不是21，而是42
dfRBFpqL$id <- id
aovRBFpq <- aov(value ~ treatment*time + Error(id/(treatment*time)), data=dfRBFpqL)
summary(aovRBFpq)
##
## Error: id
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 20 1.243 0.06216
##
## Error: id:treatment
## Df Sum Sq Mean Sq F value Pr(>F)
## treatment 1 0.1957 0.19572 2.486 0.131
## Residuals 20 1.5743 0.07871
##
## Error: id:time
## Df Sum Sq Mean Sq F value Pr(>F)
## time 2 3.154 1.5771 43.15 1.03e-10 ***
## Residuals 40 1.462 0.0365
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: id:treatment:time
## Df Sum Sq Mean Sq F value Pr(>F)
## treatment:time 2 0.2488 0.12440 3.27 0.0484 *
## Residuals 40 1.5216 0.03804
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
with(dfRBFpqL,interaction.plot(time,treatment, value, type = "b", col = c("red", "blue"), pch = c(16, 18), main = "Interaction Plot for treatment and time"))
boxplot(value ~ treatment*time, data = dfRBFpqL, col = (c("gold",
"green")), main = "treatment and time",
ylab = "value")
方差分析表明主效应time和交互效应treatment:time有显著性差异，主效应treatment无显著性差异。
宽格式数据
dfTemp <- reshape(dfRBFpqL, v.names="value", timevar="treatment",
idvar=c("id", "time"), direction="wide")
dfRBFpqW <- reshape(dfTemp, v.names=c("value.contrast", "value.ps"),
timevar="time", idvar="id", direction="wide")
fitRBFpq <- lm(cbind(value.contrast.time0,value.ps.time0,value.contrast.time24,value.ps.time24,value.contrast.time72,value.ps.time72) ~ 1,data=dfRBFpqW)
inRBFpq <- expand.grid(treatment=gl(2, 1), time=gl(3, 1))
AnovaRBFpq <- Anova(fitRBFpq, idata=inRBFpq, idesign=~treatment*time)
## Note: model has only an intercept; equivalent type-III tests substituted.
summary(AnovaRBFpq, multivariate=FALSE, univariate=TRUE)
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## SS num Df Error SS den Df F Pr(>F)
## (Intercept) 172.669 1 1.2433 20 2777.6126 < 2.2e-16 ***
## treatment 0.196 1 1.5743 20 2.4865 0.13051
## time 3.154 2 1.4619 40 43.1519 1.03e-10 ***
## treatment:time 0.249 2 1.5216 40 3.2704 0.04837 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## time 0.51779 0.001925
## treatment:time 0.70417 0.035722
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## time 0.67467 6.503e-08 ***
## treatment:time 0.77171 0.06299 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## time 0.7060381 3.485199e-08
## treatment:time 0.8238702 5.930375e-02
宽格式数据方差分析结果与长格式数据结果类似，但提供了mauchly.test检验，球形假设的条件不满足，同时给出了Greenhouse-Geisser和 Huynh-Feldt校正结果。
anova.mlm() 和 mauchly.test()
anova(fitRBFpq, M=~treatment, X=~1, idata=inRBFpq, test="Spherical")
## Analysis of Variance Table
##
##
## Contrasts orthogonal to
## ~1
##
##
## Contrasts spanned by
## ~treatment
##
## Greenhouse-Geisser epsilon: 1
## Huynh-Feldt epsilon: 1
##
## Df F num Df den Df Pr(>F) G-G Pr H-F Pr
## (Intercept) 1 2.4865 1 20 0.13051 0.13051 0.13051
## Residuals 20
anova(fitRBFpq, M=~treatment + time, X=~treatment, idata=inRBFpq, test="Spherical")
## Analysis of Variance Table
##
##
## Contrasts orthogonal to
## ~treatment
##
##
## Contrasts spanned by
## ~treatment + time
##
## Greenhouse-Geisser epsilon: 0.6747
## Huynh-Feldt epsilon: 0.7060
##
## Df F num Df den Df Pr(>F) G-G Pr H-F Pr
## (Intercept) 1 43.152 2 40 1.0301e-10 6.5031e-08 3.4852e-08
## Residuals 20
anova(fitRBFpq, M=~treatment + time + treatment:time, X=~treatment + time,idata=inRBFpq, test="Spherical")
## Analysis of Variance Table
##
##
## Contrasts orthogonal to
## ~treatment + time
##
##
## Contrasts spanned by
## ~treatment + time + treatment:time
##
## Greenhouse-Geisser epsilon: 0.7717
## Huynh-Feldt epsilon: 0.8239
##
## Df F num Df den Df Pr(>F) G-G Pr H-F Pr
## (Intercept) 1 3.2704 2 40 0.048365 0.062985 0.059304
## Residuals 20
mauchly.test(fitRBFpq, M=~treatment, X=~1, idata=inRBFpq)
##
## Mauchly's test of sphericity
## Contrasts orthogonal to
## ~1
##
## Contrasts spanned by
## ~treatment
##
##
## data: SSD matrix from lm(formula = cbind(value.contrast.time0, value.ps.time0, value.contrast.time24, SSD matrix from value.ps.time24, value.contrast.time72, value.ps.time72) ~ SSD matrix from 1, data = dfRBFpqW)
## W = 1, p-value = 1
mauchly.test(fitRBFpq, M=~treatment + time, X=~treatment, idata=inRBFpq)
##
## Mauchly's test of sphericity
## Contrasts orthogonal to
## ~treatment
##
## Contrasts spanned by
## ~treatment + time
##
##
## data: SSD matrix from lm(formula = cbind(value.contrast.time0, value.ps.time0, value.contrast.time24, SSD matrix from value.ps.time24, value.contrast.time72, value.ps.time72) ~ SSD matrix from 1, data = dfRBFpqW)
## W = 0.51779, p-value = 0.001925
mauchly.test(fitRBFpq, M=~treatment + time + treatment:time, X=~treatment + time, idata=inRBFpq)
##
## Mauchly's test of sphericity
## Contrasts orthogonal to
## ~treatment + time
##
## Contrasts spanned by
## ~treatment + time + treatment:time
##
##
## data: SSD matrix from lm(formula = cbind(value.contrast.time0, value.ps.time0, value.contrast.time24, SSD matrix from value.ps.time24, value.contrast.time72, value.ps.time72) ~ SSD matrix from 1, data = dfRBFpqW)
## W = 0.70417, p-value = 0.03572
效果大小(Effect size estimates)
EtaSq(aovRBFpq, type=1)
## eta.sq eta.sq.part eta.sq.gen
## treatment 0.02082200 0.1105771 0.03263793
## time 0.33556123 0.6833031 0.35221811
## treatment:time 0.02646937 0.1405380 0.04112598
简单效应(Simple effects)
summary(aov(value ~ treatment + Error(id/treatment), data=dfRBFpqL, subset=(time=="time0")))
##
## Error: id
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 20 1.837 0.09185
##
## Error: id:treatment
## Df Sum Sq Mean Sq F value Pr(>F)
## treatment 1 0.0219 0.0219 0.193 0.665
## Residuals 20 2.2646 0.1132
summary(aov(value ~ treatment + Error(id/treatment), data=dfRBFpqL, subset=(time=="time24")))
##
## Error: id
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 20 0.5558 0.02779
##
## Error: id:treatment
## Df Sum Sq Mean Sq F value Pr(>F)
## treatment 1 0.1670 0.16695 8.759 0.00775 **
## Residuals 20 0.3812 0.01906
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov(value ~ treatment + Error(id/treatment), data=dfRBFpqL, subset=(time=="time72")))
##
## Error: id
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 20 0.3124 0.01562
##
## Error: id:treatment
## Df Sum Sq Mean Sq F value Pr(>F)
## treatment 1 0.2557 0.2557 11.36 0.00304 **
## Residuals 20 0.4501 0.0225
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
多元方法(Multivariate approach)
summary(AnovaRBFpq, multivariate=TRUE, univariate=FALSE)
##
## Type III Repeated Measures MANOVA Tests:
##
## ------------------------------------------
##
## Term: (Intercept)
##
## Response transformation matrix:
## (Intercept)
## value.contrast.time0 1
## value.ps.time0 1
## value.contrast.time24 1
## value.ps.time24 1
## value.contrast.time72 1
## value.ps.time72 1
##
## Sum of squares and products for the hypothesis:
## (Intercept)
## (Intercept) 1036.012
##
## Sum of squares and products for error:
## (Intercept)
## (Intercept) 7.459729
##
## Multivariate Tests: (Intercept)
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.99285 2777.613 1 20 < 2.22e-16 ***
## Wilks 1 0.00715 2777.613 1 20 < 2.22e-16 ***
## Hotelling-Lawley 1 138.88063 2777.613 1 20 < 2.22e-16 ***
## Roy 1 138.88063 2777.613 1 20 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------
##
## Term: treatment
##
## Response transformation matrix:
## treatment1
## value.contrast.time0 1
## value.ps.time0 -1
## value.contrast.time24 1
## value.ps.time24 -1
## value.contrast.time72 1
## value.ps.time72 -1
##
## Sum of squares and products for the hypothesis:
## treatment1
## treatment1 1.174341
##
## Sum of squares and products for error:
## treatment1
## treatment1 9.445765
##
## Multivariate Tests: treatment
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.1105771 2.486492 1 20 0.13051
## Wilks 1 0.8894229 2.486492 1 20 0.13051
## Hotelling-Lawley 1 0.1243246 2.486492 1 20 0.13051
## Roy 1 0.1243246 2.486492 1 20 0.13051
##
## ------------------------------------------
##
## Term: time
##
## Response transformation matrix:
## time1 time2
## value.contrast.time0 1 0
## value.ps.time0 1 0
## value.contrast.time24 0 1
## value.ps.time24 0 1
## value.contrast.time72 -1 -1
## value.ps.time72 -1 -1
##
## Sum of squares and products for the hypothesis:
## time1 time2
## time1 11.956939 3.545731
## time2 3.545731 1.051457
##
## Sum of squares and products for error:
## time1 time2
## time1 3.1266691 -0.1861941
## time2 -0.1861941 1.0728858
##
## Multivariate Tests: time
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.840054 49.89489 2 19 2.74e-08 ***
## Wilks 1 0.159946 49.89489 2 19 2.74e-08 ***
## Hotelling-Lawley 1 5.252094 49.89489 2 19 2.74e-08 ***
## Roy 1 5.252094 49.89489 2 19 2.74e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------
##
## Term: treatment:time
##
## Response transformation matrix:
## treatment1:time1 treatment1:time2
## value.contrast.time0 1 0
## value.ps.time0 -1 0
## value.contrast.time24 0 1
## value.ps.time24 0 -1
## value.contrast.time72 -1 -1
## value.ps.time72 1 1
##
## Sum of squares and products for the hypothesis:
## treatment1:time1 treatment1:time2
## treatment1:time1 0.8544617 0.12687829
## treatment1:time2 0.1268783 0.01884005
##
## Sum of squares and products for error:
## treatment1:time1 treatment1:time2
## treatment1:time1 4.210046 1.088387
## treatment1:time2 1.088387 1.443103
##
## Multivariate Tests: treatment:time
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.1748240 2.012695 2 19 0.16114
## Wilks 1 0.8251760 2.012695 2 19 0.16114
## Hotelling-Lawley 1 0.2118626 2.012695 2 19 0.16114
## Roy 1 0.2118626 2.012695 2 19 0.16114
多元方法分析结果类似，主效应time具有显著性，交互效应交互效应treatment:time和treatment主效应无显著性。可以认为不同时间的记录的VEGF有差异，不同组见记录的VEGF无差异。
两级裂区设计（Two-way split-plot ANOVA）
在一个区组上，先按第一个因素（主因素或主处理）的水平数划分主因素的试验小区，主因素的小区称为主区或整区，用于安排主因素；在主区内再按第二个因素（副因素或副处理）的水平数划分小区，安排副因素，主区内的小区称副区或裂区。从整个试验所有处理组合来说，主区仅是一个不完全的区组，对第二个因素来讲，主区就是一个区组，这种设计将主区分裂成副区，称为裂区设计。
例 试验一种全身注射抗毒素对皮肤损伤的保护作用，将10只家兔随机等分两组，一组注射抗毒素，一组注射生理盐水作对照。之后，每只家兔取甲、乙两部位，随机分配分别注射低浓度毒素和高浓度毒素， 观察指标为皮肤受损直径。结果如下：
| 1 |
抗毒素A1 |
15.75 |
19.00 |
| 2 |
|
15.50 |
20.75 |
| 3 |
|
15.50 |
18.50 |
| 4 |
|
17.00 |
20.50 |
| 5 |
|
16.50 |
20.00 |
| 6 |
生理盐水A2 |
18.25 |
22.25 |
| 7 |
|
18.50 |
21.50 |
| 8 |
|
19.75 |
23.50 |
| 9 |
|
21.25 |
24.75 |
| 10 |
|
20.75 |
23.75 |
diameter<- c(15.75,15.50,15.50,17.00,16.50,19.00,20.75,18.50,
20.50,20.00,18.25,18.50,19.75,21.25,20.75,22.25,21.50,23.50,24.75,23.75)
dfSPFpqL <- data.frame(id=factor(rep(1:5, times=4)),
B=factor(rep(1:2,each=5,times=2)),
A=factor(rep(1:2,each=10)),
Diameter=diameter)
aovSPFpq <- aov(Diameter ~ A*B + Error(id/B), data=dfSPFpqL)
summary(aovSPFpq)
##
## Error: id
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 4 12.08 3.02
##
## Error: id:B
## Df Sum Sq Mean Sq F value Pr(>F)
## B 1 63.90 63.90 560.2 1.89e-05 ***
## Residuals 4 0.46 0.11
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 62.13 62.13 74.881 2.47e-05 ***
## A:B 1 0.08 0.08 0.094 0.767
## Residuals 8 6.64 0.83
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
主因素和副因素均有统计学意义。
宽数据格式
dfSPFpqW <- reshape(dfSPFpqL, v.names="Diameter", timevar="B",
idvar=c("id", "A"), direction="wide")
fitSPFpq <- lm(cbind(Diameter.1, Diameter.2) ~ A, data=dfSPFpqW)
inSPFpq <- data.frame(B=gl(2, 1))
AnovaSPFpq <- Anova(fitSPFpq, idata=inSPFpq, idesign=~B)
summary(AnovaSPFpq, multivariate=FALSE, univariate=TRUE)
##
## Univariate Type II Repeated-Measures ANOVA Assuming Sphericity
##
## SS num Df Error SS den Df F Pr(>F)
## (Intercept) 7732.3 1 17.1875 8 3599.0240 6.622e-12 ***
## A 62.1 1 17.1875 8 28.9178 0.0006636 ***
## B 63.9 1 1.9875 8 257.2201 2.291e-07 ***
## A:B 0.1 1 1.9875 8 0.3145 0.5903087
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
与长数据格式结果基本一致。
宽数据格式anova.mlm()和mauchly.test()
anova(fitSPFpq, M=~1, X=~0, idata=inSPFpq, test="Spherical")
## Analysis of Variance Table
##
##
## Contrasts orthogonal to
## ~0
##
##
## Contrasts spanned by
## ~1
##
## Greenhouse-Geisser epsilon: 1
## Huynh-Feldt epsilon: 1
##
## Df F num Df den Df Pr(>F) G-G Pr H-F Pr
## (Intercept) 1 3599.024 1 8 0.0000000 0.0000000 0.0000000
## A 1 28.918 1 8 0.0006636 0.0006636 0.0006636
## Residuals 8
anova(fitSPFpq, M=~B, X=~1, idata=inSPFpq, test="Spherical")
## Analysis of Variance Table
##
##
## Contrasts orthogonal to
## ~1
##
##
## Contrasts spanned by
## ~B
##
## Greenhouse-Geisser epsilon: 1
## Huynh-Feldt epsilon: 1
##
## Df F num Df den Df Pr(>F) G-G Pr H-F Pr
## (Intercept) 1 257.2201 1 8 0.00000 0.00000 0.00000
## A 1 0.3145 1 8 0.59031 0.59031 0.59031
## Residuals 8
mauchly.test(fitSPFpq, M=~B, X=~1, idata=inSPFpq)
##
## Mauchly's test of sphericity
## Contrasts orthogonal to
## ~1
##
## Contrasts spanned by
## ~B
##
##
## data: SSD matrix from lm(formula = cbind(Diameter.1, Diameter.2) ~ A, data = dfSPFpqW)
## W = 1, p-value = 1
效果大小(Effect size estimates)
EtaSq(aovSPFpq, type=1)
## eta.sq eta.sq.part eta.sq.gen
## B 0.4398485728 0.99291090 0.769193154
## A 0.4276311544 0.90347648 0.764154207
## A:B 0.0005377385 0.01163332 0.004057783
简单效应
#Between-subjects effect at a fixed level of the within-subjects factor
summary(aov(Diameter ~ A, data=dfSPFpqL, subset=(B==1)))
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 33.31 33.31 30.11 0.000583 ***
## Residuals 8 8.85 1.11
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov(Diameter ~ A, data=dfSPFpqL, subset=(B==2)))
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 28.90 28.900 22.39 0.00148 **
## Residuals 8 10.32 1.291
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Within-subjects effect at a fixed level of the between-subjects factor
summary(aov(Diameter ~ B + Error(id/B), data=dfSPFpqL,subset=(A==1)))
##
## Error: id
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 4 3.962 0.9906
##
## Error: id:B
## Df Sum Sq Mean Sq F value Pr(>F)
## B 1 34.23 34.23 86.24 0.000748 ***
## Residuals 4 1.59 0.40
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov(Diameter ~ B + Error(id/B), data=dfSPFpqL,subset=(A==2)))
##
## Error: id
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 4 13.22 3.306
##
## Error: id:B
## Df Sum Sq Mean Sq F value Pr(>F)
## B 1 29.76 29.76 297.6 6.63e-05 ***
## Residuals 4 0.40 0.10
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
计划的多重比较（Planned comparisons for the between-subjects factor）
mDf <- aggregate(Diameter ~ id + A, data=dfSPFpqL, FUN=mean)
aovRes <- aov(Diameter ~ A, data=mDf)
cMat <- rbind("1-2"=c(-1, 1))
summary(glht(aovRes, linfct=mcp(A=cMat), alternative="greater"),
test=adjusted("none"))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: User-defined Contrasts
##
##
## Fit: aov(formula = Diameter ~ A, data = mDf)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>t)
## 1-2 <= 0 3.5250 0.6555 5.378 0.000332 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- none method)
主因素之间如果有多个分区的情况，其差异是否具有统计学意义可采取上述方法。
再裂区设计（Three-way split-plot ANOVA）
在裂区设计中，若需再引进第三个因素时，可在副区内再分裂出第二副区，称为再裂区，然后将第三个因素的各个处理（称为副副处理），随机排列于再裂区内，这种设计称为再裂区设计（split-split plot design ）。3个以上的多因素试验采用裂区设计，试验起来很复杂，统计分析也麻烦，特别是因素之间有交互作用比较难以解释。
例 观察18例不同分化程度的贲门癌患者的癌组织、癌旁组织、远离组织中碱性磷酸酶（ALP）的变化，一级单位处理为分化度（低分化、中分化和高分化，记为A1、A2和A3），二级单位处理是组织部位（癌组织、癌旁组织、远癌组织，记为B1、B2和B3），三级单位处理是活性剂（加与不加，记为C1和C2），数据如下：
ALP <- read.csv("ALP.csv",header = T)
pander(ALP)
| A1 |
1 |
72.5 |
87.5 |
3.2 |
3.6 |
0.74 |
1.3 |
| A1 |
2 |
61.7 |
66.2 |
2.1 |
3.7 |
1.1 |
1.7 |
| A1 |
3 |
76.1 |
89.4 |
4.3 |
5 |
1.8 |
2.2 |
| A1 |
4 |
93 |
98 |
5.1 |
5.5 |
1 |
1.9 |
| A1 |
5 |
82.9 |
85.1 |
3.6 |
4.9 |
0.8 |
1.1 |
| A1 |
6 |
75.6 |
90.2 |
2.2 |
3.3 |
1 |
1.8 |
| A2 |
7 |
61.1 |
65.3 |
3.2 |
4 |
0.9 |
1.3 |
| A2 |
8 |
53.2 |
58.2 |
3.1 |
4.1 |
1 |
1.5 |
| A2 |
9 |
63.2 |
63.8 |
1.9 |
1.9 |
1.3 |
2 |
| A2 |
10 |
55.1 |
55.7 |
1.7 |
2.3 |
2.1 |
1.8 |
| A2 |
11 |
53.2 |
61.9 |
1.6 |
2.8 |
1.8 |
0.9 |
| A2 |
12 |
49.9 |
63.2 |
2.2 |
3.1 |
1.3 |
1.9 |
| A3 |
13 |
43.1 |
45.6 |
1.9 |
2.3 |
1.7 |
1.4 |
| A3 |
14 |
39.2 |
43.1 |
1.7 |
3.9 |
1.4 |
1.9 |
| A3 |
15 |
41.9 |
47.2 |
2 |
2.2 |
1.9 |
2.1 |
| A3 |
16 |
28.5 |
35.9 |
3.9 |
4.1 |
1.8 |
2.5 |
| A3 |
17 |
36.3 |
41.2 |
1.3 |
2 |
1 |
0.8 |
| A3 |
18 |
34.9 |
40 |
3.9 |
3.9 |
2.2 |
1.5 |
SPF-pq⋅r
dfSPFpq.rL <-read.csv("ALP2.csv",header = T)
dfSPFpq.rL$id <- as.factor(dfSPFpq.rL$id)
dfSPFpq.rL$B <- as.factor(dfSPFpq.rL$B)
dfSPFpq.rL$C <- as.factor(dfSPFpq.rL$C)
dfSPFpq.rL$A <- as.factor(dfSPFpq.rL$A)
aovSPFpq.r <- aov(DV ~ C*B*A + Error(id/C), data=dfSPFpq.rL)
summary(aovSPFpq.r)
##
## Error: id
## Df Sum Sq Mean Sq F value Pr(>F)
## A 2 3645 1822.5 54.76 1.28e-07 ***
## Residuals 15 499 33.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: id:C
## Df Sum Sq Mean Sq F value Pr(>F)
## C 1 167.65 167.65 48.487 4.55e-06 ***
## C:A 2 15.07 7.54 2.179 0.148
## Residuals 15 51.86 3.46
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## B 2 79851 39925 2338.301 < 2e-16 ***
## C:B 2 213 106 6.226 0.00349 **
## B:A 4 6869 1717 100.571 < 2e-16 ***
## C:B:A 4 18 4 0.258 0.90354
## Residuals 60 1024 17
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
C、B级处理、C和B交互、A和B的交互均有统计学意义。
效果大小(Effect size estimates)
EtaSq(aovSPFpq.r, type=1)
## eta.sq eta.sq.part eta.sq.gen
## A 0.0394681461 0.87954182 0.698203907
## C 0.0018153406 0.76373261 0.096175442
## C:A 0.0001631849 0.22515207 0.009474745
## B 0.8646245925 0.98733269 0.980650719
## C:B 0.0023021900 0.17186708 0.118901553
## B:A 0.0743756486 0.87020985 0.813421028
## C:B:A 0.0001909351 0.01692099 0.011068126
宽格式数据
dfSPFpq.rW <- reshape(dfSPFpq.rL, v.names="DV", timevar="C",
idvar=c("id", "B", "A"), direction="wide")
fitSPFpq.r <- lm(cbind(DV.1, DV.2) ~ A*B, data=dfSPFpq.rW)
inSPFpq.r <- data.frame(C=gl(2, 1))
AnovaSPFpq.r <- Anova(fitSPFpq.r, idata=inSPFpq.r, idesign=~C)
summary(AnovaSPFpq.r, multivariate=FALSE, univariate=TRUE)
##
## Univariate Type II Repeated-Measures ANOVA Assuming Sphericity
##
## SS num Df Error SS den Df F Pr(>F)
## (Intercept) 50045 1 1419.79 45 1586.1761 < 2.2e-16 ***
## A 3645 2 1419.79 45 57.7638 3.736e-13 ***
## B 79851 2 1419.79 45 1265.4258 < 2.2e-16 ***
## A:B 6869 4 1419.79 45 54.4264 < 2.2e-16 ***
## C 168 1 155.75 45 48.4395 1.169e-08 ***
## A:C 15 2 155.75 45 2.1772 0.1252
## B:C 213 2 155.75 45 30.7152 3.875e-09 ***
## A:B:C 18 4 155.75 45 1.2737 0.2944
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A、C、B级处理、C和B交互、A和B的的交互均有统计学意义。
宽数据格式anova.mlm()和mauchly.test()
anova(fitSPFpq.r, M=~1, X=~0, idata=inSPFpq.r, test="Spherical")
## Analysis of Variance Table
##
##
## Contrasts orthogonal to
## ~0
##
##
## Contrasts spanned by
## ~1
##
## Greenhouse-Geisser epsilon: 1
## Huynh-Feldt epsilon: 1
##
## Df F num Df den Df Pr(>F) G-G Pr H-F Pr
## (Intercept) 1 1586.176 1 45 0.000e+00 0.000e+00 0.000e+00
## A 2 57.764 2 45 3.736e-13 3.736e-13 3.736e-13
## B 2 1265.426 2 45 0.000e+00 0.000e+00 0.000e+00
## A:B 4 54.426 4 45 1.100e-16 1.100e-16 1.100e-16
## Residuals 45
anova(fitSPFpq.r, M=~C, X=~1, idata=inSPFpq.r, test="Spherical")
## Analysis of Variance Table
##
##
## Contrasts orthogonal to
## ~1
##
##
## Contrasts spanned by
## ~C
##
## Greenhouse-Geisser epsilon: 1
## Huynh-Feldt epsilon: 1
##
## Df F num Df den Df Pr(>F) G-G Pr H-F Pr
## (Intercept) 1 48.4395 1 45 0.00000 0.00000 0.00000
## A 2 2.1772 2 45 0.12516 0.12516 0.12516
## B 2 30.7152 2 45 0.00000 0.00000 0.00000
## A:B 4 1.2737 4 45 0.29439 0.29439 0.29439
## Residuals 45
mauchly.test(fitSPFpq.r, M=~C, X=~1, idata=inSPFpq.r)
##
## Mauchly's test of sphericity
## Contrasts orthogonal to
## ~1
##
## Contrasts spanned by
## ~C
##
##
## data: SSD matrix from lm(formula = cbind(DV.1, DV.2) ~ A * B, data = dfSPFpq.rW)
## W = 1, p-value = 1
SPF-p⋅qr
aovSPFp.qr <- aov(DV ~ C*B*A + Error(id/(A*B)), data=dfSPFpq.rL)
## Warning in aov(DV ~ C * B * A + Error(id/(A * B)), data = dfSPFpq.rL):
## Error() model is singular
summary(aovSPFp.qr)
##
## Error: id
## Df Sum Sq Mean Sq F value Pr(>F)
## A 2 3645 1822.5 54.76 1.28e-07 ***
## Residuals 15 499 33.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: id:B
## Df Sum Sq Mean Sq F value Pr(>F)
## B 2 79851 39925 1301.08 < 2e-16 ***
## B:A 4 6869 1717 55.96 1.74e-13 ***
## Residuals 30 921 31
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## C 1 167.65 167.65 48.440 1.17e-08 ***
## C:B 2 212.61 106.31 30.715 3.88e-09 ***
## C:A 2 15.07 7.54 2.177 0.125
## C:B:A 4 17.63 4.41 1.274 0.294
## Residuals 45 155.75 3.46
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
C、C和B交互均有统计学意义。
效果大小(Effect size estimates)
EtaSq(aovSPFp.qr, type=1)
## eta.sq eta.sq.part eta.sq.gen
## A 0.0394681461 0.87954182 0.698203907
## B 0.8646245925 0.98860254 0.980650719
## B:A 0.0743756486 0.88181536 0.813421028
## C 0.0018153406 0.51840521 0.096175442
## C:B 0.0023021900 0.57718835 0.118901553
## C:A 0.0001631849 0.08822602 0.009474745
## C:B:A 0.0001909351 0.10170330 0.011068126
宽格式数据
dfW1 <- reshape(dfSPFpq.rL, v.names="DV", timevar="C",
idvar=c("id", "B", "A"), direction="wide")
dfSPFp.qrW <- reshape(dfW1, v.names=c("DV.1", "DV.2"),
timevar="B", idvar=c("id", "A"), direction="wide")
fitSPFp.qr <- lm(cbind(DV.1.1, DV.2.1, DV.1.2, DV.2.2, DV.1.3, DV.2.3) ~ A,
data=dfSPFp.qrW)
inSPFp.qr <- expand.grid(B=gl(3, 1), C=gl(2, 1))
AnovaSPFp.qr <- Anova(fitSPFp.qr, idata=inSPFp.qr, idesign=~B*C)
summary(AnovaSPFp.qr, multivariate=FALSE, univariate=TRUE)
##
## Univariate Type II Repeated-Measures ANOVA Assuming Sphericity
##
## SS num Df Error SS den Df F Pr(>F)
## (Intercept) 50045 1 499.20 15 1503.755 < 2.2e-16 ***
## A 3645 2 499.20 15 54.762 1.277e-07 ***
## B 20295 2 299.74 30 1015.640 < 2.2e-16 ***
## A:B 1704 4 299.74 30 42.633 5.788e-12 ***
## C 40475 1 481.24 15 1261.578 6.827e-16 ***
## A:C 3528 2 481.24 15 54.976 1.245e-07 ***
## B:C 19461 2 295.35 30 988.349 < 2.2e-16 ***
## A:B:C 1670 4 295.35 30 42.410 6.185e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## B 0.77596 0.16939
## A:B 0.77596 0.16939
## B:C 0.76557 0.15413
## A:B:C 0.76557 0.15413
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## B 0.81697 < 2.2e-16 ***
## A:B 0.81697 4.007e-10 ***
## B:C 0.81009 < 2.2e-16 ***
## A:B:C 0.81009 4.960e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## B 0.9031466 9.448097e-26
## A:B 0.9031466 5.438711e-11
## B:C 0.8939930 2.349101e-25
## A:B:C 0.8939930 7.133228e-11
宽数据格式anova.mlm()和mauchly.test()
anova(fitSPFp.qr, M=~1, X=~0,idata=inSPFp.qr, test="Spherical")
## Analysis of Variance Table
##
##
## Contrasts orthogonal to
## ~0
##
##
## Contrasts spanned by
## ~1
##
## Greenhouse-Geisser epsilon: 1
## Huynh-Feldt epsilon: 1
##
## Df F num Df den Df Pr(>F) G-G Pr H-F Pr
## (Intercept) 1 1503.755 1 15 0.0000e+00 0.0000e+00 0.0000e+00
## A 2 54.762 2 15 1.2772e-07 1.2772e-07 1.2772e-07
## Residuals 15
anova(fitSPFp.qr, M=~B, X=~1,idata=inSPFp.qr, test="Spherical")
## Analysis of Variance Table
##
##
## Contrasts orthogonal to
## ~1
##
##
## Contrasts spanned by
## ~B
##
## Greenhouse-Geisser epsilon: 0.8170
## Huynh-Feldt epsilon: 0.9031
##
## Df F num Df den Df Pr(>F) G-G Pr H-F Pr
## (Intercept) 1 1015.640 2 30 0.0000e+00 0.0000e+00 0.0000e+00
## A 2 42.633 4 30 5.7885e-12 4.0068e-10 5.4387e-11
## Residuals 15
anova(fitSPFp.qr, M=~B + C, X=~B,idata=inSPFp.qr, test="Spherical")
## Analysis of Variance Table
##
##
## Contrasts orthogonal to
## ~B
##
##
## Contrasts spanned by
## ~B + C
##
## Greenhouse-Geisser epsilon: 1
## Huynh-Feldt epsilon: 1
##
## Df F num Df den Df Pr(>F) G-G Pr H-F Pr
## (Intercept) 1 1261.578 1 15 0.0000e+00 0.0000e+00 0.0000e+00
## A 2 54.976 2 15 1.2448e-07 1.2448e-07 1.2448e-07
## Residuals 15
anova(fitSPFp.qr, M=~B + C + B:C, X=~B + C,idata=inSPFp.qr, test="Spherical")
## Analysis of Variance Table
##
##
## Contrasts orthogonal to
## ~B + C
##
##
## Contrasts spanned by
## ~B + C + B:C
##
## Greenhouse-Geisser epsilon: 0.8101
## Huynh-Feldt epsilon: 0.8940
##
## Df F num Df den Df Pr(>F) G-G Pr H-F Pr
## (Intercept) 1 988.35 2 30 0.0000e+00 0.0000e+00 0.0000e+00
## A 2 42.41 4 30 6.1847e-12 4.9597e-10 7.1332e-11
## Residuals 15
mauchly.test(fitSPFp.qr, M=~B, X=~1,idata=inSPFp.qr)
##
## Mauchly's test of sphericity
## Contrasts orthogonal to
## ~1
##
## Contrasts spanned by
## ~B
##
##
## data: SSD matrix from lm(formula = cbind(DV.1.1, DV.2.1, DV.1.2, DV.2.2, DV.1.3, DV.2.3) ~ SSD matrix from A, data = dfSPFp.qrW)
## W = 0.77596, p-value = 0.1694
mauchly.test(fitSPFp.qr, M=~B + C, X=~B,
idata=inSPFp.qr)
##
## Mauchly's test of sphericity
## Contrasts orthogonal to
## ~B
##
## Contrasts spanned by
## ~B + C
##
##
## data: SSD matrix from lm(formula = cbind(DV.1.1, DV.2.1, DV.1.2, DV.2.2, DV.1.3, DV.2.3) ~ SSD matrix from A, data = dfSPFp.qrW)
## W = 1, p-value = 1
mauchly.test(fitSPFp.qr, M=~B + C + B:C, X=~B + C,
idata=inSPFp.qr)
##
## Mauchly's test of sphericity
## Contrasts orthogonal to
## ~B + C
##
## Contrasts spanned by
## ~B + C + B:C
##
##
## data: SSD matrix from lm(formula = cbind(DV.1.1, DV.2.1, DV.1.2, DV.2.2, DV.1.3, DV.2.3) ~ SSD matrix from A, data = dfSPFp.qrW)
## W = 0.76557, p-value = 0.1541
混合模型重复测量方差分析（Mixed-effects models for repeated-measures ANOVA）
在分析数据时，考虑一个因素和它的不同水平对结果变量的影响，称之为这个因素不同水平对因变量的效应。这种效应不是固定效应就是随机效应，当参数能被认为是固定的常数时，这种因素所产生的效应为固定效应，当参数有随机变量的特征时，称之为随机效应。当模型中有多个因素，一部分产生固定效应，一部分产生随机效应，这样的模型就称为混合效应模型。重复测量中的单次测量为低水平，个体为高水平，建立的模型如下：\[Y=X\beta +Z\gamma +\epsilon\],\(X\)为已知设计矩阵，\(\beta\)为固定效应参数构成的未知向量，\(\epsilon\)为未知的随机误差向量，其元素不必为独立同分布。\(Y\)和\(\gamma\)均为正态随机变量。
例 d1为长格式的重复观测数据，因变量为Y，自变量为Xw1、Xb1和Xb2，w表示组内因子，b表示组间因子，id为标示变量。d2为长格式的重复观测数据，因变量为Y，自变量为Xw1、Xw2、Xb1和Xb2，w表示组内因子，b表示组间因子，id为标示变量。
d1 <- read.csv("d1.csv",header = T)
d2 <- read.csv("d2.csv",header = T)
单因素重复测量方差分析(One-way repeated measures ANOVA, RB-p design)
常规分析(Conventional analysis using aov())
summary(aov(Y ~ Xw1 + Error(id/Xw1), data=d1))
##
## Error: id
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 79 75040 949.9
##
## Error: id:Xw1
## Df Sum Sq Mean Sq F value Pr(>F)
## Xw1 2 5756 2878.2 3.363 0.0371 *
## Residuals 158 135211 855.8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
结果显示自变量Xw1有统计学意义。
混合效应分析(Mixed-effects analysis)
对重复测量的数据有个假设就是重复测量的数据间的关系是相同的，这就是我们所说的compound symmetry。但在实际中，往往会违背这个假设，特别是当临床试验的时间特别长或各个测量的时间点的间隔不相同时，这是因为间隔时间长的两个点的测量值之间的关系往往不如间隔时间短的两个点的测量值之间的关系紧密。
#没有明确是否符合compound symmetry假设
anova(lme(Y ~ Xw1, random=~1 | id, method="ML", data=d1))
## numDF denDF F-value p-value
## (Intercept) 1 158 2554.7564 <.0001
## Xw1 2 158 3.3633 0.0371
#符合compound symmetry假设
lmeFit <- lme(Y ~ Xw1, random=~1 | id, correlation=corCompSymm(form=~1|id),
method="ML", data=d1)
anova(lmeFit)
## numDF denDF F-value p-value
## (Intercept) 1 158 2554.7627 <.0001
## Xw1 2 158 3.3633 0.0371
anova(lme(Y ~ Xw1, random=list(id=pdCompSymm(~Xw1-1)), method="REML", data=d1))
## numDF denDF F-value p-value
## (Intercept) 1 158 2554.7564 <.0001
## Xw1 2 158 3.3633 0.0371
结果显示自变量Xw1有统计学意义。
lme4包lmer()方法
fitF <- lmer(Y ~ Xw1 + (1|id), data=d1)
anova(fitF)
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Xw1 2 5756.4 2878.2 3.3633
fitR <- lmer(Y ~ 1 + (1|id), data=d1)
library(pbkrtest)
KRmodcomp(fitF, fitR)
## F-test with Kenward-Roger approximation; computing time: 0.11 sec.
## large : Y ~ Xw1 + (1 | id)
## small : Y ~ 1 + (1 | id)
## stat ndf ddf F.scaling p.value
## Ftest 3.3633 2.0000 158.0000 1 0.03712 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
结果显示自变量Xw1有统计学意义。
AIC值比较
library(AICcmodavg)
AICc(fitF)
## [1] 2304.445
aictab(cand.set=list(fitR, fitF),
modnames=c("restricted", "full"),
sort=FALSE, second.ord=FALSE)
## Warning in aictab.AIClmerMod(cand.set = list(fitR, fitF), modnames = c("restricted", :
## Model selection for fixed effects is only appropriate with ML estimation:
## REML (default) should only be used to select random effects for a constant set of fixed effects
##
## Model selection based on AIC :
##
## K AIC Delta_AIC AICWt Res.LL
## restricted 3 2316.36 12.17 0 -1155.18
## full 5 2304.19 0.00 1 -1147.09
多重比较(基于multcomp包的glht()方法)
contr <- glht(lmeFit, linfct=mcp(Xw1="Tukey"))
summary(contr)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lme.formula(fixed = Y ~ Xw1, data = d1, random = ~1 | id, correlation = corCompSymm(form = ~1 |
## id), method = "ML")
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## B - A == 0 10.36375 4.59638 2.255 0.0624 .
## C - A == 0 10.41417 4.59638 2.266 0.0607 .
## C - B == 0 0.05042 4.59638 0.011 0.9999
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
confint(contr) #置信区间
##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lme.formula(fixed = Y ~ Xw1, data = d1, random = ~1 | id, correlation = corCompSymm(form = ~1 |
## id), method = "ML")
##
## Quantile = 2.3437
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## B - A == 0 10.36375 -0.40868 21.13618
## C - A == 0 10.41417 -0.35827 21.18660
## C - B == 0 0.05042 -10.72202 10.82285
Xw1变量B-A和C-A比较有差异。
双因素重复测量方差分析(Two-way repeated measures ANOVA ,RBF-pq design)
常规分析（Conventional analysis using aov()）
summary(aov(Y ~ Xw1*Xw2 + Error(id/(Xw1*Xw2)), data=d2))
##
## Error: id
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 79 225120 2850
##
## Error: id:Xw1
## Df Sum Sq Mean Sq F value Pr(>F)
## Xw1 2 17269 8635 3.363 0.0371 *
## Residuals 158 405633 2567
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: id:Xw2
## Df Sum Sq Mean Sq F value Pr(>F)
## Xw2 2 9859 4929 1.616 0.202
## Residuals 158 481938 3050
##
## Error: id:Xw1:Xw2
## Df Sum Sq Mean Sq F value Pr(>F)
## Xw1:Xw2 4 6118 1530 0.603 0.661
## Residuals 316 802069 2538
结果显示自变量Xw1有统计学意义。
Mixed-effects analysis
nlme包lme方法（Using lme() from package nlme）
anova(lme(Y ~ Xw1*Xw2, random=list(id=pdBlocked(list(~1, pdIdent(~Xw1-1), pdIdent(~Xw2-1)))),
method="ML", data=d2))
## numDF denDF F-value p-value
## (Intercept) 1 632 2440.2286 <.0001
## Xw1 2 632 3.3889 0.0344
## Xw2 2 632 1.6522 0.1924
## Xw1:Xw2 4 632 0.6003 0.6626
#符合compound symmetry假设
anova(lme(Y ~ Xw1*Xw2,
random=list(id=pdBlocked(list(~1, pdCompSymm(~Xw1-1), pdCompSymm(~Xw2-1)))),
method="ML", data=d2))
## numDF denDF F-value p-value
## (Intercept) 1 632 2554.7545 <.0001
## Xw1 2 632 3.3633 0.0352
## Xw2 2 632 1.6160 0.1995
## Xw1:Xw2 4 632 0.6026 0.6609
结果显示自变量Xw1有统计学意义。
lme4包lmer()方法(Using lmer() from package lme4)
anova(lmer(Y ~ Xw1*Xw2 + (1|id) + (1|Xw1:id) + (1|Xw2:id), data=d2))
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Xw1 2 17269.2 8634.6 3.3889
## Xw2 2 8419.5 4209.7 1.6522
## Xw1:Xw2 4 6118.0 1529.5 0.6003
根据F值判断自变量Xw1有统计学意义。
两级裂区设计的方差分析(Two-way split-plot-factorial ANOVA ,SPF-p⋅q design)
常规分析(Conventional analysis using aov())
summary(aov(Y ~ Xb1*Xw1 + Error(id/Xw1), data=d1))
##
## Error: id
## Df Sum Sq Mean Sq F value Pr(>F)
## Xb1 1 5335 5335 5.97 0.0168 *
## Residuals 78 69705 894
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: id:Xw1
## Df Sum Sq Mean Sq F value Pr(>F)
## Xw1 2 5756 2878 3.541 0.03133 *
## Xb1:Xw1 2 8414 4207 5.176 0.00666 **
## Residuals 156 126797 813
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
混合效应分析(Mixed-effects)
nlme包lme方法(Using lme() from package nlme)
#没有明确是否符合compound symmetry假设
anova(lme(Y ~ Xb1*Xw1, random=~1 | id, method="ML", data=d1))
## numDF denDF F-value p-value
## (Intercept) 1 156 2715.4696 <.0001
## Xb1 1 78 5.9697 0.0168
## Xw1 2 156 3.5411 0.0313
## Xb1:Xw1 2 156 5.1762 0.0067
#符合compound symmetry假设
anova(lme(Y ~ Xb1*Xw1, random=~1 | id, correlation=corCompSymm(form=~1|id),
method="ML", data=d1))
## numDF denDF F-value p-value
## (Intercept) 1 156 2715.4765 <.0001
## Xb1 1 78 5.9697 0.0168
## Xw1 2 156 3.5411 0.0313
## Xb1:Xw1 2 156 5.1762 0.0067
anova(lme(Y ~ Xb1*Xw1, random=list(id=pdCompSymm(~Xw1-1)), method="REML", data=d1))
## numDF denDF F-value p-value
## (Intercept) 1 156 2715.4698 <.0001
## Xb1 1 78 5.9697 0.0168
## Xw1 2 156 3.5411 0.0313
## Xb1:Xw1 2 156 5.1762 0.0067
自变量Xb1、Xw1和其交互效应均有统计学意义。
lme4包lmer()方法（Using lmer() from package lme4）
anova(lmer(Y ~ Xb1*Xw1 + (1|id), data=d1))
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Xb1 1 4852.2 4852.2 5.9697
## Xw1 2 5756.4 2878.2 3.5411
## Xb1:Xw1 2 8414.4 4207.2 5.1762
三级裂区设计的方差分析(Three-way split-plot-factorial ANOVA ,SPF-pq⋅r design)
常规分析(Conventional analysis using aov())
summary(aov(Y ~ Xb1*Xb2*Xw1 + Error(id/Xw1), data=d1))
##
## Error: id
## Df Sum Sq Mean Sq F value Pr(>F)
## Xb1 1 5335 5335 7.468 0.00781 **
## Xb2 1 7246 7246 10.144 0.00210 **
## Xb1:Xb2 1 8169 8169 11.436 0.00114 **
## Residuals 76 54290 714
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: id:Xw1
## Df Sum Sq Mean Sq F value Pr(>F)
## Xw1 2 5756 2878 4.116 0.018166 *
## Xb1:Xw1 2 8414 4207 6.016 0.003058 **
## Xb2:Xw1 2 11336 5668 8.105 0.000452 ***
## Xb1:Xb2:Xw1 2 9167 4583 6.554 0.001861 **
## Residuals 152 106294 699
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
自变量Xb1、Xb2、Xw1和其交互效应均有统计学意义。
混合效应模型(Mixed-effects analysis)
nlme包lme方法(Using lme() from package nlme)
#没有明确是否符合compound symmetry假设
anova(lme(Y ~ Xb1*Xb2*Xw1, random=~1 | id, method="ML", data=d1))
## numDF denDF F-value p-value
## (Intercept) 1 152 3397.096 <.0001
## Xb1 1 76 7.468 0.0078
## Xb2 1 76 10.144 0.0021
## Xw1 2 152 4.116 0.0182
## Xb1:Xb2 1 76 11.436 0.0011
## Xb1:Xw1 2 152 6.016 0.0031
## Xb2:Xw1 2 152 8.105 0.0005
## Xb1:Xb2:Xw1 2 152 6.554 0.0019
#符合compound symmetry假设
anova(lme(Y ~ Xb1*Xb2*Xw1, random=~1 | id,
correlation=corCompSymm(form=~1 | id), method="ML", data=d1))
## numDF denDF F-value p-value
## (Intercept) 1 152 3397.098 <.0001
## Xb1 1 76 7.468 0.0078
## Xb2 1 76 10.144 0.0021
## Xw1 2 152 4.116 0.0182
## Xb1:Xb2 1 76 11.436 0.0011
## Xb1:Xw1 2 152 6.016 0.0031
## Xb2:Xw1 2 152 8.105 0.0005
## Xb1:Xb2:Xw1 2 152 6.554 0.0019
anova(lme(Y ~ Xb1*Xb2*Xw1,
random=list(id=pdBlocked(list(~1, pdCompSymm(~Xw1-1)))),
method="ML", data=d1))
## Warning in lme.formula(Y ~ Xb1 * Xb2 * Xw1, random = list(id =
## pdBlocked(list(~1, : fewer observations than random effects in all level 1
## groups
## numDF denDF F-value p-value
## (Intercept) 1 152 3397.098 <.0001
## Xb1 1 76 7.468 0.0078
## Xb2 1 76 10.144 0.0021
## Xw1 2 152 4.116 0.0182
## Xb1:Xb2 1 76 11.436 0.0011
## Xb1:Xw1 2 152 6.016 0.0031
## Xb2:Xw1 2 152 8.105 0.0005
## Xb1:Xb2:Xw1 2 152 6.554 0.0019
自变量Xb1、Xb2、Xw1和其交互效应均有统计学意义。
lme4包lmer()方法（Using lmer() from package lme4）
anova(lmer(Y ~ Xb1*Xb2*Xw1 + (1|id), data=d1))
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Xb1 1 5222.5 5222.5 7.4682
## Xb2 1 7093.5 7093.5 10.1437
## Xw1 2 5756.4 2878.2 4.1158
## Xb1:Xb2 1 7997.0 7997.0 11.4356
## Xb1:Xw1 2 8414.4 4207.2 6.0163
## Xb2:Xw1 2 11335.9 5668.0 8.1052
## Xb1:Xb2:Xw1 2 9166.6 4583.3 6.5541
根据F值判断自变量Xb1、Xb2、Xw1和其交互效应均有统计学意义。
三级裂区设计的方差分析Three-way split-plot-factorial ANOVA (SPF-p⋅qr design)
常规分析(Conventional analysis using aov())
summary(aov(Y ~ Xb1*Xw1*Xw2 + Error(id/(Xw1*Xw2)), data=d2))
##
## Error: id
## Df Sum Sq Mean Sq F value Pr(>F)
## Xb1 1 16005 16005 5.97 0.0168 *
## Residuals 78 209116 2681
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: id:Xw1
## Df Sum Sq Mean Sq F value Pr(>F)
## Xw1 2 17269 8635 3.541 0.03133 *
## Xb1:Xw1 2 25243 12622 5.176 0.00666 **
## Residuals 156 380390 2438
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: id:Xw2
## Df Sum Sq Mean Sq F value Pr(>F)
## Xw2 2 9859 4929 1.604 0.204
## Xb1:Xw2 2 2462 1231 0.400 0.671
## Residuals 156 479476 3074
##
## Error: id:Xw1:Xw2
## Df Sum Sq Mean Sq F value Pr(>F)
## Xw1:Xw2 4 6118 1530 0.601 0.662
## Xb1:Xw1:Xw2 4 7609 1902 0.747 0.561
## Residuals 312 794460 2546
Xb1、Xw1和Xw2无统计学意义。
混合效应模型(Mixed-effects analysis)
nlme包lme方法(Using lme() from package nlme)
#没有明确是否符合compound symmetry假设
anova(lme(Y ~ Xb1*Xw1*Xw2,
random=list(id=pdBlocked(list(~1, pdIdent(~Xw1-1), pdIdent(~Xw2-1)))),
method="ML", data=d2))
## numDF denDF F-value p-value
## (Intercept) 1 624 2473.9526 <.0001
## Xb1 1 78 5.4387 0.0223
## Xw1 2 624 3.4396 0.0327
## Xw2 2 624 1.6751 0.1881
## Xb1:Xw1 2 624 5.0278 0.0068
## Xb1:Xw2 2 624 0.4183 0.6584
## Xw1:Xw2 4 624 0.6093 0.6561
## Xb1:Xw1:Xw2 4 624 0.7577 0.5531
#符合compound symmetry假设
anova(lme(Y ~ Xb1*Xw1*Xw2,
random=list(id=pdBlocked(list(~1, pdCompSymm(~Xw1-1), pdCompSymm(~Xw2-1)))),
method="ML", data=d2))
## numDF denDF F-value p-value
## (Intercept) 1 624 2715.3731 <.0001
## Xb1 1 78 5.9695 0.0168
## Xw1 2 624 3.4396 0.0327
## Xw2 2 624 1.6038 0.2020
## Xb1:Xw1 2 624 5.0278 0.0068
## Xb1:Xw2 2 624 0.4005 0.6702
## Xw1:Xw2 4 624 0.6093 0.6561
## Xb1:Xw1:Xw2 4 624 0.7577 0.5531
Xb1、Xw1和其交互作用有统计学意义。
lme4包lmer()方法（Using lmer() from package lme4）
anova(lmer(Y ~ Xb1*Xw1*Xw2 + (1|id) + (1|Xw1:id) + (1|Xw2:id), data=d2))
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Xb1 1 13653.2 13653.2 5.4387
## Xw1 2 17269.2 8634.6 3.4396
## Xw2 2 8410.1 4205.1 1.6751
## Xb1:Xw1 2 25243.2 12621.6 5.0278
## Xb1:Xw2 2 2100.1 1050.1 0.4183
## Xw1:Xw2 4 6118.0 1529.5 0.6093
## Xb1:Xw1:Xw2 4 7608.8 1902.2 0.7577
四级裂区设计的方差分析(Four-way split-plot-factorial ANOVA ,SPF-pq⋅rs design)
常规分析(Conventional analysis using aov())
summary(aov(Y ~ Xb1*Xb2*Xw1*Xw2 + Error(id/(Xw1*Xw2)), data=d2))
##
## Error: id
## Df Sum Sq Mean Sq F value Pr(>F)
## Xb1 1 16005 16005 7.468 0.00781 **
## Xb2 1 21738 21738 10.144 0.00210 **
## Xb1:Xb2 1 24507 24507 11.436 0.00114 **
## Residuals 76 162871 2143
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: id:Xw1
## Df Sum Sq Mean Sq F value Pr(>F)
## Xw1 2 17269 8635 4.116 0.018166 *
## Xb1:Xw1 2 25243 12622 6.016 0.003058 **
## Xb2:Xw1 2 34008 17004 8.105 0.000452 ***
## Xb1:Xb2:Xw1 2 27500 13750 6.554 0.001861 **
## Residuals 152 318882 2098
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: id:Xw2
## Df Sum Sq Mean Sq F value Pr(>F)
## Xw2 2 9859 4929 1.728 0.18110
## Xb1:Xw2 2 2462 1231 0.432 0.65031
## Xb2:Xw2 2 11822 5911 2.072 0.12945
## Xb1:Xb2:Xw2 2 34080 17040 5.974 0.00318 **
## Residuals 152 433574 2852
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: id:Xw1:Xw2
## Df Sum Sq Mean Sq F value Pr(>F)
## Xw1:Xw2 4 6118 1529 0.610 0.6558
## Xb1:Xw1:Xw2 4 7609 1902 0.759 0.5530
## Xb2:Xw1:Xw2 4 24545 6136 2.447 0.0465 *
## Xb1:Xb2:Xw1:Xw2 4 7595 1899 0.757 0.5539
## Residuals 304 762320 2508
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Xb1、Xb2、Xb1:Xb2、Xw1、Xb1:Xw1、Xb2:Xw1和Xb1:Xb2:Xw1均有统计学意义。
混合效应模型(Mixed-effects analysis)
nlme包lme方法(Using lme() from package nlme)
#没有明确是否符合compound symmetry假设
anova(lme(Y ~ Xb1*Xb2*Xw1*Xw2,
random=list(id=pdBlocked(list(~1, pdIdent(~Xw1-1), pdIdent(~Xw2-1)))),
method="ML", data=d2))
## numDF denDF F-value p-value
## (Intercept) 1 608 2782.9276 <.0001
## Xb1 1 76 6.1180 0.0156
## Xb2 1 76 8.3098 0.0051
## Xw1 2 608 3.6417 0.0268
## Xw2 2 608 1.8843 0.1528
## Xb1:Xb2 1 76 9.3682 0.0031
## Xb1:Xw1 2 608 5.3232 0.0051
## Xb2:Xw1 2 608 7.1714 0.0008
## Xb1:Xw2 2 608 0.4705 0.6249
## Xb2:Xw2 2 608 2.2595 0.1053
## Xw1:Xw2 4 608 0.6451 0.6305
## Xb1:Xb2:Xw1 2 608 5.7991 0.0032
## Xb1:Xb2:Xw2 2 608 6.5138 0.0016
## Xb1:Xw1:Xw2 4 608 0.8023 0.5240
## Xb2:Xw1:Xw2 4 608 2.5880 0.0359
## Xb1:Xb2:Xw1:Xw2 4 608 0.8008 0.5249
#符合compound symmetry假设
anova(lme(Y ~ Xb1*Xb2*Xw1*Xw2,
random=list(id=pdBlocked(list(~1, pdCompSymm(~Xw1-1), pdCompSymm(~Xw2-1)))),
method="ML", data=d2))
## numDF denDF F-value p-value
## (Intercept) 1 608 3113.1727 <.0001
## Xb1 1 76 6.8440 0.0107
## Xb2 1 76 9.2959 0.0032
## Xw1 2 608 3.6924 0.0255
## Xw2 2 608 1.7281 0.1785
## Xb1:Xb2 1 76 10.4799 0.0018
## Xb1:Xw1 2 608 5.3973 0.0047
## Xb2:Xw1 2 608 7.2713 0.0008
## Xb1:Xw2 2 608 0.4315 0.6497
## Xb2:Xw2 2 608 2.0722 0.1268
## Xw1:Xw2 4 608 0.6541 0.6242
## Xb1:Xb2:Xw1 2 608 5.8798 0.0030
## Xb1:Xb2:Xw2 2 608 5.9737 0.0027
## Xb1:Xw1:Xw2 4 608 0.8134 0.5168
## Xb2:Xw1:Xw2 4 608 2.6240 0.0339
## Xb1:Xb2:Xw1:Xw2 4 608 0.8119 0.5178
Xb1、Xb2、Xb1:Xb2、Xw1、Xb1:Xw1、Xb2:Xw1和Xb1:Xb2:Xw1均有统计学意义。
lme4包lmer()方法（Using lmer() from package lme4）
anova(lmer(Y ~ Xb1*Xb2*Xw1*Xw2 + (1|id) + (1|Xw1:id) + (1|Xw2:id), data=d2))
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Xb1 1 14506 14506.1 6.1180
## Xb2 1 19703 19702.9 8.3098
## Xw1 2 17269 8634.6 3.6417
## Xw2 2 8935 4467.7 1.8843
## Xb1:Xb2 1 22212 22212.5 9.3682
## Xb1:Xw1 2 25243 12621.6 5.3232
## Xb2:Xw1 2 34008 17003.9 7.1714
## Xb1:Xw2 2 2231 1115.7 0.4705
## Xb2:Xw2 2 10715 5357.5 2.2595
## Xw1:Xw2 4 6118 1529.5 0.6451
## Xb1:Xb2:Xw1 2 27500 13749.9 5.7991
## Xb1:Xb2:Xw2 2 30889 15444.5 6.5138
## Xb1:Xw1:Xw2 4 7609 1902.2 0.8023
## Xb2:Xw1:Xw2 4 24545 6136.2 2.5880
## Xb1:Xb2:Xw1:Xw2 4 7595 1898.7 0.8008
根据F值判断Xb1、Xb2、Xb1:Xb2、Xw1、Xb1:Xw1、Xb2:Xw1和Xb1:Xb2:Xw1均有统计学意义。