R语言学习系列27方差分析讲课教案Word文档下载推荐.docx
《R语言学习系列27方差分析讲课教案Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《R语言学习系列27方差分析讲课教案Word文档下载推荐.docx（20页珍藏版）》请在冰豆网上搜索。
试验的因素（factor）可以是品种、人员、方法、时间、地区等等，因素所处的状态叫水平（level）。
在每一个因素下面可以分成若干水平。
（3）因素间的交互影响
多因素的试验设计，有时需要分析因素间的交互影响（interaction），2个因素间的交互影响称为一级交互影响（A×
B）；
3个因素间的交互影响称为二级交互影响（A×
B×
C）。
当交互影响项呈现统计不显著时，表明各个因素独立，当呈现统计显著时，就需要列出这个交互影响项的效应，以助于作出正确的统计推断。
举例解释上述概念：
要考察焦虑症的治疗疗效，一个因素是治疗方案，有2种治疗方案，即该因素有2个水平；
（治疗方案称为组间因子，因为每个患者只能被分配到一个组别中，没有患者同时接受两种治疗）；
再考虑一个因素治疗时间，也有两个水平：
治疗5周和治疗6个月，同一患者在5周和6个月不止一次地被测量（两次），称为重复测量（治疗时间称为组内因子，因为每个患者在所有水平下都进行了测量）。
建立方差分析模型时，既要考虑两个因素治疗方案和治疗时间（主效应），又要考虑治疗方案和时间的交互影响（交互效应），此时即两因素混合模型方差分析。
当某个因素的各个水平下的因变量的均值呈现统计显著性差异时，必要时可作两两水平间的比较，称为均值间的两两比较。
二、R语言实现
方差分析对数据的要求：
满足正态性（来自同一正态总体）和方差齐性（各组方差相等），在这两个条件下，若各组有差异，则只可能是来自影响因素的不同水平。
用aov（）函数进行方差分析，基本格式为：
aov（formula,data=NULL,projections=FALSE,qr=TRUE,
contrasts=NULL,...）
其中，formula为方差分析公式；
data为数据框；
projection设置是否返回预测结果；
qr设置是否返回QR分解结果；
contrasts为公式中一些因子的列表。
formula公式的表示：
（y为因变量，ABC为分组因子）
符号
用法
~
分隔符号，左边为响应变量，右边为解释变量
eg：
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，代码y~.可展开为y~A+B+C
常见研究设计的表达式：
（小写字母表示定量变量，大写字母表
示组别因子，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）
注意：
非均衡设计时或存在协变量时，效应项的顺序对结果影响较大，越基础的效应越需要放在表达式前面，首先是协变量、然后是主效应、接着是双因素的交互项，再接着是三因素的交互项。
若研究不是正交的，一定要谨慎设置效应的顺序。
有三种类型的方法可以分解y~A+B+A:
B右边各效应对y所解释的方差：
类型I（序贯型）
效应根据表达式中先出现的效应做调整。
A不做调整，B根据A调整，A:
B交互项根据A和B调整。
类型II（分层型）
效应根据同水平或低水平的效应做调整。
A根据B调整，B依据A调整，A:
B交互项同时根据A和B调整。
类型III（边界型）
每个效应根据模型其他各效应做相应调整。
A根据B和A:
B做调整，A:
R默认调用类型I方法，其他软件（比如SAS和SPSS）默认调用类型III方法。
car包中的Anova（）函数（不要与标准anova（）函数混淆）提供了使用类型II或类型III方法的选项，而aov（）函数使用的是类型I方法。
若想使结果与其他软件（如SAS和SPSS）提供的结果保持一致，可以使用Anova（）函数。
三、单因素方差分析
1个因变量，1个影响因素：
总差异Yij=平均差异μ+因素差异αi+随机差异εij
例1比较4种品牌的胶合板的耐磨性，各抽取5个样品，相同转速磨损相同时间测得磨损深度（mm），比较4个品牌胶合板的耐磨性有无差异？
部分数据如下（ex27_ex1.Rdata）：
setwd（"
E:
/办公资料/R语言/R语言学习系列/codes"
）
load（"
ex27_ex1.Rdata"
head（datas）
wearbrand
12.30A
22.32A
32.40A
42.45A
52.58A
62.35B
attach（datas）
table（brand）#各组的样本数
brand
ABCD
5555
aggregate（wear,by=list（brand）,mean）#各组均值
Group.1x
1A2.410
2B2.404
3C2.046
4D2.572
aggregate（wear,by=list（brand）,sd）#各组标准差
1A0.11269428
2B0.11760102
3C0.11216060
4D0.03271085
library（car）
qqPlot（lm（wear~brand,data=datas）,simulate=TRUE）#用Q-Q图检验数据的正态性
leveneTest（wear~as.factor（brand）,data=datas）#方差齐性检验
Levene'
sTestforHomogeneityofVariance（center=median）
DfFvaluePr（>
F）
group30.69870.5664
16
fit<
-aov（wear~brand,data=datas）#单因素方差分析（检验组间差异）
summary（fit）
DfSumSqMeanSqFvaluePr（>
F）
brand30.73980.2466024.553.15e-06***
Residuals160.16070.01005
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
说明：
方差齐性检验，原假设H0：
方差齐，p值=0.5664>
0.05,故接受原假设，即方差齐。
单因素方差分析结果，brand是因素，Residuals是残差，各列依次为自由度、平方和、均方和、F统计量，p值=3.15e-06<
0.05,拒绝原假设，即不同品牌的磨损（均值）有显著差别。
library（gplots）
plotmeans（wear~brand,xlab="
品牌"
ylab="
磨损"
）#图形展示带95%置信区间的各组均值
通过前面的分析知道，不同品牌的磨损（均值）有显著差别，但并不知道哪个品牌与其它品牌有显著差别。
TukeyHSD（）函数提供了对各组均值差异的成对检验。
TukeyHSD（fit）
Tukeymultiplecomparisonsofmeans
95%family-wiseconfidencelevel
Fit:
aov（formula=wear~brand,data=datas）
$brand
difflwruprpadj
B-A-0.006-0.187353450.17535350.9996826
C-A-0.364-0.54535345-0.18264650.0001610
D-A0.162-0.019353450.34335350.0886142
C-B-0.358-0.53935345-0.17664650.0001929
D-B0.168-0.013353450.34935350.0744337
D-C0.5260.344646550.70735350.0000019
说明：
可以看出（H0：
无差异），B与A的差异非常不显著，C与A、C与B、D与C的差异非常显著。
multcomp包中的glht（）函数提供了更为全面的多重均值比较方法。
library（multcomp）
tuk<
-glht（fit,linfct=mcp（brand="
Tukey"
））
#注意datas$brand必须是因子型
summary（tuk）
SimultaneousTestsforGeneralLinearHypotheses
MultipleComparisonsofMeans:
TukeyContrasts
LinearHypotheses:
EstimateStd.ErrortvaluePr（>
|t|）
B-A==0-0.006000.06339-0.0950.9997
C-A==0-0.364000.06339-5.742<
0.001***
D-A==00.162000.063392.5560.0886.
C-B==0-0.358000.06339-5.648<
D-B==00.168000.063392.6500.0743.
D-C==00.526000.063398.298<
（Adjustedpvaluesreported--single-stepmethod）
plot（cld（tuk,level=0.05）,col="
lightgrey"
标记相同字母（标记b）的品牌ABD认为是无显著差异，在同一亚组，而品牌C（标记a）与另外三个品牌有显著差异。
另外，也可以进行多重t检验，使用函数：
pairwise.t.test（x,g,p.adjust.method=,...）
其中，x为因变量，g为因子型的分组变量；
p.adjust.method设置p值的修正方法，由于多次重复t检验会大大增加犯第一类错误的概率，为此要进行p值的修正，使用bonferroni法修正效果较好。
pairwise.t.test（wear,brand,p.adjust.method="
bonferroni"
PairwisecomparisonsusingttestswithpooledSD
data:
wearandbrand
ABC
B1.00000--
C0.000180.00022-
D0.126950.104742.1e-06
Pvalueadjustmentmethod:
bonferroni
原假设H0:
无差异，可见A与B无差异，C与ABD有显著差异。
最后，方差分析对离群点非常敏感，检验是否有离群点：
outlierTest（fit）
NoStudentizedresidualswithBonferonnip<
0.05
Largest|rstudent|:
rstudentunadjustedp-valueBonferonnip
92.5281030.0231820.46364
经检验无离群点。
三、两因素方差分析
1个因变量，2个影响因素：
总差异Yijk=平均差异μ+因素1差异αi+因素2差异βi
+因素1,2交互作用差异γij+随机差异εijk
例2研究60只豚鼠的牙齿生长数据，按2种喂食方法：
橙汁、维生素C，各喂食方法中抗坏血酸含量都有3个水平：
0.5mg、1mg、2mg，分配为6组，每组各10只，牙齿长度为因变量。
做两因素方差分析。
attach（ToothGrowth）
head（ToothGrowth）
lensuppdose
14.2VC0.5
211.5VC0.5
37.3VC0.5
45.8VC0.5
56.4VC0.5
610.0VC0.5
table（supp,dose）#各组样本数相同，即为均衡设计
dose
supp0.512
OJ101010
VC101010
aggregate（len,by=list（supp,dose）,mean）#计算各组均值
Group.1Group.2x
1OJ0.513.23
2VC0.57.98
3OJ1.022.70
4VC1.016.77
5OJ2.026.06
6VC2.026.14
aggregate（len,by=list（supp,dose）,sd）#计算各组标准差
1OJ0.54.459709
2VC0.52.746634
3OJ1.03.910953
4VC1.02.515309
5OJ2.02.655058
6VC2.04.797731
bartlett.test（len~supp,data=ToothGrowth）#关于因素supp的方差齐性检验
Bartletttestofhomogeneityofvariances
lenbysupp
Bartlett'
sK-squared=1.4217,df=1,p-value=0.2331
bartlett.test（len~dose,data=ToothGrowth）#关于因素dose的方差齐性检验
lenbydose
sK-squared=0.66547,df=2,p-value=0.717
fit<
-aov（len~supp*dose,data=ToothGrowth）
#做两因素方差分析，考虑全部效应
supp1205.4205.412.3170.000894***
dose12224.32224.3133.415<
2e-16***
supp:
dose188.988.95.3330.024631*
Residuals56933.616.7
Signif.codes:
可以看出，主效应supp和dose都非常显著（p值都远小于0.05），交互效应也显著（p值=0.0246<
0.05）。
若交互作用不显著，可以可以只做去掉交互效应的方差分析。
图形化展示两因素方差分析的交互效应：
par（mfrow=c（1,2））
interaction.plot（dose,supp,len,type="
b"
col=c（"
red"
"
blue"
）,pch=c（16,18）,main="
InteractionbetweenDoseandSupp"
interaction.plot（supp,dose,len,type="
col=c（"
InteractionbetweenSuppandDose"
有一个图的线有交叉，说明有交互作用。
可以看出随着橙汁和维生素C中的抗坏血酸剂量的增加，牙齿长度变长；
。
对于0.5mg
和1mg剂量，橙汁比维生素C更能促进牙齿生长；
对于2mg剂量的抗坏血酸，两种喂食方法下牙齿长度增长相同。
也可以用HH包中的interaction2wt（）函数（也适合三因素方差分析）来展示更全面的可视化结果：
library（HH）
interaction2wt（len~supp*dose）
三、重复测量方差分析
重复测量方差分析，即受试者被测量不止一次。
例3（1个组内1个组间因子的重复测量）在某浓度CO2的环境中，对寒带植物（来自魁北克）和非寒带植物的（来自密西西比）光合作用率进行比较。
因变量uptake为CO2吸收量，自变量Type（组间因子）为植物类型，自变量conc（组内因子）为七种水平的CO2浓度。
attach（CO2）
head（CO2）#注意CO2是长格式的数据
PlantTypeTreatmentconcuptake
1Qn1Quebecnonchilled9516.0
2Qn1Quebecnonchilled17530.4
3Qn1Quebecnonchilled25034.8
4Qn1Quebecnonchilled35037.2
5Qn1Quebecnonchilled50035.3
6Qn1Quebecnonchilled67539.2
w1b1<
-subset（CO2,Treatment=="
chilled"
）#先只考虑寒带植物
-aov（uptake~（conc*Type）+Error（Plant/conc）,data=w1b1）
Error:
Plant
Type12667.22667.260.410.00148**
Residuals4176.644.1
Plant:
conc
conc1888.6888.6215.460.000125***
conc:
Type1239.2239.258.010.001595**
Residuals416.54.1
Within
Residuals3086928.97
在0.01的显著水平下，主效应“类型”（p值=0.00148）和“浓度”（p值=0.000125）以及交叉效应“类型*浓度”（p值=0.001595）都非常显著。
attach（w1b1）
interaction.plot（conc,Type,uptake,type="
）,pch=c（16,18）,main="
InteractionPlotforPlantTypeandConcentration"
boxplot（uptake~Type*conc,data=w1b1,col=（c（"
gold"
green"
））,main="
升级会员 