另一是研商中施加的对结果产生影响的可控因素,效应依据悉明式中先出现的成效做调解ca88亚洲城网站

方差深入分析(Analysis
of Variance,简称ANOVA),又称“
变异数剖判”,是Escort.A.Fisher发明的,用于四个及七个以上
样本均数差其他 分明性查验
由于各样因素的震慑,钻探所得的多少表现波动状。变成波动的缘故可分为两类,一是不可控的随机因素,另一是钻探中施加的对结果产生震慑的可控因素。方差深入分析是从观测变量的方差出手,商量广大 调节变量中哪些变量是对考查变量有分明影响的变量。

第楚辞方差深入分析

方差分析有七个基本点的前提:

9.2 ANOVA 模型拟合

  • 调控变量分歧档案的次序下考查变量的总体布满为正态遍及
  • 调整变量分化档期的顺序下考查变量的完整具有同样的方差

9.2.1 aov()函数

1.单要素方差分析

aov(formula, data = NULL, projections =FALSE, qr = TRUE,

1.1概述

二个调节变量的两样程度是不是对考查变量产生了断定影响。

ca88亚洲城网站 1

ca88亚洲城网站 2

F>>1–调节变量对观测值爆发刚烈影响。

F分布密度曲线

x<-rf(1000,df1[i],df2[i])

set.seed(12345)
x<-rnorm(1000,0,1)
Ord<-order(x,decreasing=FALSE)
#order()的返回值是对应“排名”的元素所在向量中的位置
x<-x[Ord]
y<-dnorm(x,0,1)
plot(x,y,xlim=c(-1,5),ylim=c(0,2),type="l",ylab="密度",main="标准正态分布与不同自由度下的F分布密度函数",lwd=1.5)
#######不同自由度的F分布
df1<-c(10,15,30,100)
df2<-c(10,20,25,110)
for(i in 1:4){
 x<-rf(1000,df1[i],df2[i])
 Ord<-order(x,decreasing=FALSE)
 x<-x[Ord]
 y<-df(x,df1[i],df2[i])
 lines(x,y,lty=i+1)
}
legend("topright",title="自由度",c("标准正态分布",paste(df1,df2,sep="-")),lty=1:5)

  

ca88亚洲城网站 3

 

contrasts = NULL, …)

1.2数学模型

ca88亚洲城网站 4

 

1.3 R程序

ca88亚洲城网站 5

ca88亚洲城网站 6

CarData<-read.table(file="CarData.txt",header=TRUE)
CarData$ModelYear<-as.factor(CarData$ModelYear)
aov(MPG~ModelYear,data=CarData)
OneWay<-aov(MPG~ModelYear,data=CarData)
anova(OneWay)
summary(OneWay)

  

> aov(MPG~ModelYear,data=CarData)
Call:
aov(formula = MPG ~ ModelYear, data = CarData)

Terms:
ModelYear Residuals
Sum of Squares 10401.78 13850.79
Deg. of Freedom 12 385

Residual standard error: 5.998007
Estimated effects may be unbalanced

> OneWay<-aov(MPG~ModelYear,data=CarData)
> anova(OneWay)
Analysis of Variance Table

Response: MPG
Df Sum Sq Mean Sq F value Pr(>F)
ModelYear 12 10402 866.82 24.094 < 2.2e-16 ***

ca88亚洲城网站 7

Residuals 385 13851 35.98

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> summary(OneWay)
Df Sum Sq Mean Sq F value Pr(>F)
ModelYear 12 10402 866.8 24.09 <2e-16 ***

ca88亚洲城网站 8

Residuals 385 13851 36.0

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

F=24.09,拒绝原若是,分裂年份车的型号的MPG总体均值存在明显差距。

 

1.4各总体均值的可视化 

plotmeans()

install.packages("gplots")
library("gplots")
plotmeans(MPG~ModelYear,data=CarData,p=0.95,use.t=TRUE,xlab="年代车型",ylab="平均MPG",main="不同年代车型MPG总体均值变化折线图(95%置信区间)")

  ca88亚洲城网站 9

9.2.2 表达式中各个的逐一

1.5单成分方差深入分析的前提如若

  • 调整变量区别水平下考查变量的总体布满为正态布满
  • 调控变量不一致水平下考查变量的完整具有一样的方差

y ~ A + B + A:B

完全正态核算

Q-Qplot / K-S test

ca88亚洲城网站 10

 unique(CarData$ModelYear)

###########检验方差分析的前提假设(正态性检验一)
par(mfrow=c(3,5),mar=c(4,4,4,4))
for(i in unique(CarData$ModelYear)){
 T<-subset(CarData,CarData$ModelYear==i)
 qqnorm(T$MPG,main=paste(i,"年车型mpg Q-Q图"),cex=0.7)
 qqline(T$MPG,distribution = qnorm)
}

 

############或者
library("lattice")
qqmath(~MPG|ModelYear,data=CarData)

  

 

ca88亚洲城网站 11

 

K-S test

有三种档案的次序的办法能够表达等式右侧各职能对y所讲授的方差。陆风X8私下认可类型I

全部方差齐性核实

ca88亚洲城网站 12

ks.test(数值型向量名,”pnorm”)

###########检验方差分析的前提假设(正态性检验二)
for(i in unique(CarData$ModelYear)){
 T<-subset(CarData,CarData$ModelYear==i)
 R<-ks.test(T$MPG,"pnorm")
 print(R)
}

One-sample Kolmogorov-Smirnov test

data: T$MPG
D = 1, p-value < 2.2e-16
alternative hypothesis: two-sided

李樯态遍布有拨云见日差别。

类型I(序贯型)

到处差齐性核算

#########检验方差分析的前提假设(方差齐性性检验)
library("car")
leveneTest(CarData$MPG,CarData$ModelYear, center=mean)

  

Levene’s Test for Homogeneity of Variance (center = mean)
Df F value Pr(>F)
group 12 1.7173 0.06103 .

效用依据表明式中先出现的功用做调度。A不做调节,B依照A调治,A:B交互项依据A和

385

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

不能够拒绝原借使,具备同样总体方差。

B调整。

1.6单成分方差查证中的多重相比较印证

类型II(分层型)

LSD检验

ca88亚洲城网站 13

ca88亚洲城网站 14

 

功用根据同程度或低端次的效能做调节。A依据B调度,B依附A调度,A:B交互项同期根

Tukey HSD检验

ca88亚洲城网站 15

ca88亚洲城网站 16

对差异车型MPG的多级比较印证

##############多重比较检验
CarData<-read.table(file="CarData.txt",header=TRUE)
CarData$ModelYear<-as.factor(CarData$ModelYear)
OneWay<-aov(MPG~ModelYear,data=CarData)
OneWay$coefficients
TukeyHSD(OneWay,ordered=FALSE,conf.level=0.95)
Result<-TukeyHSD(OneWay,ordered=TRUE,conf.level=0.95)
LineCol<-vector()
LineCol[Result[[1]][,4]<0.05]<-2
LineCol[Result[[1]][,4]>=0.05]<-1
#检验显著时,设置为红色;否则为黑色
par(las=2)
par(mar=c(5,8,4,2))
plot(Result,cex.axis=0.5,col=LineCol)

 

> TukeyHSD(OneWay,ordered=FALSE,conf.level=0.95)
Tukey multiple comparisons of means
95% family-wise confidence level

Fit: aov(formula = MPG ~ ModelYear, data = CarData)

$ModelYear
diff lwr upr p adj
71-70 3.5603448 -1.737584676 8.8582743 0.5621668
72-70 1.0246305 -4.273298962 6.3225600 0.9999872
73-70 -0.5896552 -5.466537903 4.2872276 0.9999999
74-70 5.0140485 -0.333563533 10.3616606 0.0912557
75-70 2.5770115 -2.630295013 7.7843180 0.9127501
76-70 3.8838742 -1.170630163 8.9383786 0.3375465
77-70 5.6853448 0.387415324 10.9832743 0.0229529
78-70 6.3714559 1.382000010 11.3609119 0.0018064
79-70 7.4034483 2.152197475 12.6546991 0.0002657
80-70 16.0068966 10.755645751 21.2581474 0.0000000
81-70 12.6448276 7.393576785 17.8960784 0.0000000
82-70 14.0200222 8.854163329 19.1858812 0.0000000

71,72,73,75…..-70无差异

 

 ca88亚洲城网站 17

 

据A和B调整。

1.7效益深入分析

ca88亚洲城网站 18

###################单因素方差分析的功效分析
library("pwr")
pwr.anova.test(k=13,f=0.25,sig.level=0.05,power=0.8)

  

Balanced one-way analysis of variance power calculation

k = 13
n = 22.15691
f = 0.25
sig.level = 0.05
power = 0.8

NOTE: n is number in each group

 

效应量是指由于因素引起的差距,是衡量管理功效大小的指标。与明显性核准昨今分化,这几个目的不受样本容积潜移默化。它代表差别管理下的总体均值之间差异的轻重,能够在分歧研商时期展开相比较。

#############效应量和样本量的关系曲线
library("pwr")
ES<-seq(from=0.1,to=0.8,by=0.01)
SampleSize<-matrix(nrow=length(ES),ncol=8)
for(i in 3:10){
 for(j in 1:length(ES)){
  result<-pwr.anova.test(k=i,f=ES[j],sig.level=0.05,power=0.8)
  SampleSize[j,i-2]<-ceiling(result$n)
  }
 }
plot(SampleSize[,1],ES,type="l",ylab="效应量",xlab="样本量(每个水平)",main="单因素方差分析(Alpha=0.05,Power=0.8)")
for(i in 2:8){
 lines(SampleSize[,i],ES,type="l",col=i)
}
legend("topright",title="水平数",paste("k",3:10,sep="="),lty=1,col=1:8)

  ca88亚洲城网站 19

ca88亚洲城网站 20

类型III(边界型)

1.8置换查验

###################单因素方差分析的置换检验
install.packages("lmPerm")
library("lmPerm")
CarData<-read.table(file="CarData.txt",header=TRUE)
CarData$ModelYear<-as.factor(CarData$ModelYear)
OneWay<-aov(MPG~ModelYear,data=CarData)
anova(OneWay)#
Fit<-aovp(MPG~ModelYear,data=CarData,perm="Prob")
anova(Fit)#两者结果一致

  

种种效应依据模型其余各职能做相应调解。A根据B和A:B做调节,A:B交互项根据A和B

2.单因素协方差剖析

调整。

2.1概述

除去调整变量外,其余变量也会对观测值产生潜移默化。为了更确切的研商调整变量对观衡量的震慑,应解除别的变量的影响。协方差分析将其余变量作为协变量,并在解除协变量对观测值影响的规格下研商调整变量的震慑。

9.3 单因素方差深入分析

2.2数学模型

ca88亚洲城网站 21

> library(multcomp)

2.3Evoque函数和示范

 

 

ca88亚洲城网站 22

比方说,在解除weight那么些体协会变量的熏陶下,核准车型对MPG的震慑

################单因素协方差分析
CarData<-read.table(file="CarData.txt",header=TRUE)
CarData$ModelYear<-as.factor(CarData$ModelYear)
Result<-aov(MPG~weight+ModelYear,data=CarData)
anova(Result)

e<-CarData$MPG-Result$fitted.values   #剔除协变量影响后的残差
anova(aov(e~CarData$ModelYear))

  

Analysis of Variance Table

Response: MPG
Df Sum Sq Mean Sq F value Pr(>F)
weight 1 16777.8 16777.8 1665.526 < 2.2e-16 ***
ModelYear 12 3606.6 300.5 29.835 < 2.2e-16 ***

> attach(cholesterol)

Residuals 384 3868.2 10.1

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Analysis of Variance Table

Response: e
Df Sum Sq Mean Sq F value Pr(>F)
CarData$ModelYear 12 0.0 0.000 0 1
Residuals 385 3868.2 10.047

ca88亚洲城网站 23

多种相比印证

install.packages("effects")
library("effects")
effect("ModelYear",Result)#调整后的MPG值
plot(effect("ModelYear",Result))
tapply(CarData$MPG,INDEX=CarData$ModelYear,FUN=mean)#调整前的MPG值

  ca88亚洲城网站 24

前提核查

对调整变量在分歧等级次序下与协变量的涉嫌是不是同样且无鲜明差别

coplot(MPG~weight|ModelYear,data=CarData)

  ca88亚洲城网站 25

 

> table(trt) #各组样本大小

3.多因素方差深入分析

trt

3.1概述

ca88亚洲城网站 26

ca88亚洲城网站 27

ca88亚洲城网站 28

1time 2times4times drugD drugE

3.3奥迪Q5函数和演示

ca88亚洲城网站 29

  • x~A+B+A:B 双要素方差,个中X~A+B中A和B是见仁见智因素的等级次序因子(不思索交互功能),A:B代表交互成效生成的因数
  • x~A*B*C==A+B+C+A:B+A:C+B:C+A:B:C
  • y~(A+B+C)^2==A+B+C+A:B+A:C+B:C
  • y~. 解析除了y以外别的任何因素对调查变量的影响

注意:~右侧调整变量的排列顺序很要紧。 ep, A+B与B+A是分化的

ca88亚洲城网站 30

 

################多因素方差分析
CarData<-read.table(file="CarData.txt",header=TRUE)
CarData$ModelYear<-as.factor(CarData$ModelYear)
CarData$cylinders<-as.factor(CarData$cylinders)
table(CarData$cylinders)

Result<-aov(MPG~cylinders+ModelYear+cylinders:ModelYear,data=CarData)
anova(Result)

  

Response: MPG
Df Sum Sq Mean Sq F value Pr(>F)
cylinders 4 15454.8 3863.7 280.3702 <2e-16 ***
ModelYear 12 3423.7 285.3 20.7036 <2e-16 ***
cylinders:ModelYear 26 482.0 18.5 1.3451 0.1236

10 10 10 10 10

Residuals 355 4892.1 13.8

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

互相项对其无影响–重新组织模型

############多因素方差分析的非饱和模型
Result<-aov(MPG~cylinders+ModelYear,data=CarData)
anova(Result)

  交互可视化

ca88亚洲城网站 31

######可视化交互效应
interaction.plot(CarData$ModelYear,CarData$cylinders,CarData$MPG,type="b",main="气缸数和车型对MPG的交互效应",xlab="车型",ylab="MPG均值")

  ca88亚洲城网站 32

ca88亚洲城网站 33

换到核准

CarData<-read.table(file="CarData.txt",header=TRUE)
CarData$ModelYear<-as.factor(CarData$ModelYear)
CarData$cylinders<-as.factor(CarData$cylinders)
Result<-aov(MPG~cylinders+ModelYear,data=CarData)
anova(Result)
library("lmPerm")

Fit<-aovp(MPG~cylinders+ModelYear,data=CarData)
anova(Fit)

  

Analysis of Variance Table

Response: MPG
Df Sum Sq Mean Sq F value Pr(>F)
cylinders 4 15454.8 3863.7 273.919 < 2.2e-16 ***
ModelYear 12 3423.7 285.3 20.227 < 2.2e-16 ***

> aggregate(response,by=list(trt),FUN=mean)#各组均值

Residuals 381 5374.1 14.1

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

Analysis of Variance Table

Response: MPG
Df R Sum Sq R Mean Sq Iter Pr(Prob)
cylinders 4 8476.7 2119.17 5000 < 2.2e-16 ***
ModelYear 12 3423.7 285.31 5000 < 2.2e-16 ***

Group.1 x

Residuals 381 5374.1 14.11

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

4.小结

1 1time 5.78197

 ca88亚洲城网站 34

2 2times 9.22497

3 4times 12.37478

4 drugD 15.36117

5 drugE 20.94752

> aggregate(response,by=list(trt),FUN=sd) #各组规范差

Group.1 x

1 1time 2.878113

2 2times 3.483054

3 4times 2.923119

4 drugD 3.454636

5 drugE 3.345003

 

> fit<-aov(response~trt)

> summary(fit) #查实组间差距(ANOVA)

Df SumSq Mean Sq F value Pr(>F)

trt 41351.4 337.8 32.43 9.82e-13 ***

Residuals 45 468.8 10.4


Signif. codes:

0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> library(gplots)

>
plotmeans(response~trt,xlab=”treatment”,ylab=”response”,main=”meanplot\nwith
95% CI”)

#绘图各组均值及其置信区间的图纸

> detach(cholesterol)

ca88亚洲城网站 35

9.3.1 多种比较

TukeyHSD()函数提供了对各组均值差异的成对核准。但要注意TukeyHSD()函数与HH包存在包容性难点:若载入HH包,TukeyHSD()函数将会失灵。使用detach(“package::HH”)将它从寻觅路线中删去,然后再调用TukeyHSD()

> detach(“package:HH”, unload=TRUE)

> TukeyHSD(fit)

Tukey multiplecomparisons of means

95% family-wise confidence level

Fit: aov(formula = response ~ trt)

$trt

diff lwr upr p adj

2times-1time 3.44300 -0.6582817 7.5442820.1380949

4times-1time 6.59281 2.4915283 10.6940920.0003542

drugD-1time 9.57920 5.4779183 13.680482 0.0000003

drugE-1time 15.16555 11.0642683 19.266832 0.0000000

4times-2times 3.14981 -0.9514717 7.2510920.2050382

drugD-2times 6.13620 2.0349183 10.2374820.0009611

drugE-2times 11.72255 7.6212683 15.8238320.0000000

drugD-4times 2.98639 -1.1148917 7.0876720.2512446

drugE-4times 8.57274 4.4714583 12.6740220.0000037

drugE-drugD 5.58635 1.4850683 9.687632 0.0030633

> par(las=2)

> par(mar=c(5,8,4,2))

> plot(TukeyHSD(fit))

ca88亚洲城网站 36

multcomp包中的glht()函数提供了多种均值相比较更为周到的章程,既适用于线性模型,也适用于广义线性模型.

> library(multcomp)

> par(mar=c(5,4,6,2))

> tuk<-glht(fit,linfct=mcp(trt=”Tukey”))

> plot(cld(tuk,level=.5),col=”lightgrey”)

ca88亚洲城网站 37

9.3.2 评估查证的借使规范

> library(car)

>
qqPlot(lm(response~trt,data=cholesterol),simulate=TRUE,main=”Q-Qplot”,labels=FALSE)

ca88亚洲城网站 38

Bartlett检验:

> bartlett.test(response~trt,data=cholesterol)

 

Bartlett test ofhomogeneity of variances

 

data: response by trt

Bartlett’s K-squared = 0.5797, df = 4,

p-value = 0.9653

方差齐性分析对离群点特别敏感。可应用car包中的outlierTest()函数来检查实验离群点:

> library(car)

> outlierTest(fit)

 

No Studentized residuals withBonferonni p < 0.05

Largest |rstudent|:

rstudent unadjusted p-value Bonferonni p

19 2.251149 0.029422 NA

9.4 单因素协方差深入分析

单元素协方差深入分析(ANCOVA)扩张了单因素方差剖析(ANOVA),包罗贰个或五个定量的

协变量

> data(litter,package=”multcomp”)

> attach(litter)

> table(dose)

dose

0 5 50 500

20 19 18 17

> aggregate(weight,by=list(dose),FUN=mean)

Group.1 x

1 0 32.30850

2 5 29.30842

3 50 29.86611

4 500 29.64647

> fit<-aov(weight~gesttime+dose)

> summary(fit)

Df Sum Sq Mean Sq Fvalue 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

出于采纳了协变量,要博得调解的组均值——即去除协变量效应后的组均值。可使

用effects包中的effects()函数来总括调治的均值:

> library(effects)

> effect(“dose”,fit)

dose effect

dose

0 5 50 500

32.35367 28.87672 30.56614 29.33460

对用户定义的比较的层层相比较

> library(multcomp)

> contrast<-rbind(“no drugvs. drug”=c(3,-1,-1,-1))

> summary(glht(fit,linfct=mcp(dose=contrast)))

 

Simultaneous Tests for General LinearHypotheses

 

Multiple Comparisons of Means: User-definedContrasts

 

 

Fit: aov(formula = weight ~ gesttime +dose)

 

Linear Hypotheses:

Estimate Std. Error tvalue

no drug vs. drug == 0 8.284 3.209 2.581

Pr(>|t|)

no drug vs. drug == 0 0.012 *


Signif. codes:

0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’0.1 ‘ ’ 1

(Adjusted p values reported –single-step method)

9.4.1 评估核查的要是标准

稽查回归斜率的同质性

> library(multcomp)

> fit2<-aov(weight~gesttime*dose,data=litter)

> summary(fit2)

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

9.4.2 结果可视化

HH包中的ancova()函数能够绘制因变量、协变量和因子之间的关系图。

> library(HH)

> 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

ca88亚洲城网站 39

9.5 双成分方差解析

双因素ANOVA

> attach(ToothGrowth)

> table(supp,dose)

dose

supp 0.5 1 2

OJ 10 10 10

VC 10 10 10

> aggregate(len,by=list(supp,dose),FUN=mean)

Group.1 Group.2 x

1 OJ 0.5 13.23

2 VC 0.5 7.98

3 OJ 1.0 22.70

4 VC 1.0 16.77

5 OJ 2.0 26.06

6 VC 2.0 26.14

> aggregate(len,by=list(supp,dose),FUN=sd)

Group.1 Group.2 x

1 OJ 0.5 4.459709

2 VC 0.5 2.746634

3 OJ 1.0 3.910953

4 VC 1.0 2.515309

5 OJ 2.0 2.655058

6 VC 2.0 4.797731

> fit<-aov(len~supp*dose)

> summary(fit)

Df Sum Sq Mean Sq F value Pr(>F)

supp 1 205.3 205.3 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

table语句的预管理申明该安排是平衡设计(各安插单元中样本大小都同样),aggregate

言辞管理可获得各单元的均值和标准差。用summary()函数得到方差深入分析表,能够见见主效应

(supp和dose)和互相功用都充足鲜明。有多种艺术对结果开始展览可视化处理。interaction.plot()函数来展现双成分方差剖析的竞相成效。

>
interaction.plot(dose,supp,len,type=”b”,col=c(“red”,”blue”),pch=c(16,18),main=”interactionbetween
dose and supplement type”)

ca88亚洲城网站 40

还足以用gplots包中的plotmeans()函数来显示交互功效。

> library(gplots)

>
plotmeans(len~interaction(supp,dose,sep=””),connect=list(c(1,3,5),c(2,4,6)),col=c(“red”,”darkgreen”),main=”interactionplot
with 95%CIs”,xlab=”treatment and dose combination”)

ca88亚洲城网站 41

用HH包中的interaction2wt()函数来可视化结果,图形对私下顺序的因子

规划的主效应和相互功用都会开始展览展现

> library(HH)

> interaction2wt(len~supp*dose)

ca88亚洲城网站 42

9.6 重复测量方差深入分析

含一个组间因子和叁个组内因子的双重度量方差剖判

w1b1<-subset(CO2,Treatment==”chilled”)

w1b1

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=”InteractionPlot
for Plant Type and Concentration”))

ca88亚洲城网站 43

boxplot(uptake~Type*conc,data=w1b1,col=c(“gold”,”green”),

main=”Chilled Quebec andMississippi Plants”,

ylab=”Carbon dioxide uptake rateumol/m^2 sec”)

ca88亚洲城网站 44

9.7 多元方差解析

当因变量(结果变量)不唯有一个时,可用多元方差深入分析(MANOVA)对它们同一时候举行辨析。

library(MASS)

head(UScereal)

attach(UScereal)

y<-cbind(calories,fat,sugars)

aggregate(y,by=list(shelf),FUN=mean)

cov(y)

fit<-manova(y~shelf)

summary(fit)

summary.aov(fit)

9.7.1 评估若是核查

单成分多元方差剖析有七个前提若是,贰个是体系正态性,贰个是方差—协方差矩阵同质性。

力排众议补充

若有多少个p*1的多级正态随机向量x,均值为μ,协方差矩阵为Σ,那么x与μ的马氏距离

的平方服从自由度为p的卡方布满。Q-Q图展现卡方布满的分位数,横纵坐标分别是样本量与

马氏距离平方值。倘诺点全部落在斜率为1、截距项为0的直线上,则评释数据遵循多元春态

分布。

查查多元日态性

> 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=”q-qplot assessing
multivariate normality”,ylab=”mahalanobis d2″)

> identify(coord$x,coord$y,labels=row.names(UScereal))

 

ca88亚洲城网站 45

9.7.2 稳健多元方差分析

假定多元春态性或许方差—协方差均值借使都不满意,又大概你怀恋多元离群点,那么可以

思索用肃穆或非参数版本的MANOVA查验。稳健单因素MANOVA可透过rrcov包中的

Wilks.test()函数完成。vegan包中的adonis()函数则提供了非参数MANOVA的等同形式。

> library(rrcov)

> Wilks.test(y,shelf,method=”mcd”)

9.8 用回归来做ANOVA

> library(multcomp)

> levels(cholesterol$trt)

[1] “1time” “2times” “4times””drugD”

[5] “drugE”

> fit.aov<-aov(response~trt,data=cholesterol)

> summary(fit.aov)

Df Sum Sq Mean Sq F value Pr(>F)

trt 4 1351.4 337.8 32.43 9.82e-13 ***

Residuals 45 468.8 10.4


Signif. codes:

0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’0.1 ‘ ’ 1

用lm()函数拟合同样的模子

> fit.lm<-lm(response~trt,data=cholesterol)

> summary(fit.lm)

Call:

lm(formula = response ~ trt, data =cholesterol)

 

Residuals:

Min 1Q Median 3Q Max

-6.5418 -1.9672 -0.0016 1.8901 6.6008

 

Coefficients:

Estimate Std. Error t valuePr(>|t|)

(Intercept) 5.782 1.021 5.665 9.78e-07 ***

trt2times 3.443 1.443 2.385 0.0213 *

trt4times 6.593 1.443 4.568 3.82e-05 ***

trtdrugD 9.579 1.443 6.637 3.53e-08 ***

trtdrugE 15.166 1.443 10.507 1.08e-13 ***


Signif. codes:

0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’0.1 ‘ ’ 1

 

Residual standard error: 3.227 on 45degrees of freedom

Multiple R-squared: 0.7425, AdjustedR-squared: 0.7196

F-statistic: 32.43 on 4 and 45DF, p-value: 9.819e-13

ca88亚洲城网站 46
ca88亚洲城网站 47

相关文章