R语言单因素方差分析及两两比较

一、导语

两个样本均数的比较用t检验，那么多个样本均数的比较应该采用什么方法分析呢？就是接下来介绍的方差分析。方差分析由统计学家R.A.Fisher提出，又称为F检验。是通过对数据变异的分析来推断两个或者多个样本均数所代表总体均数是否有差别的一种统计学方法。

二、方差分析思想

方差分析的基本思想是把全部观察值间的变异按设计类型的不同，分解成两个或多个组成部分，然后将各个部分的变异与随机误差进行比较，以判断各部分的变异是否具有统计学意义。

三、单因素方差分析变异来源

完全随机化设计资料的方差分析，变异包括以下三种变异：总变异、组间变异、组内变异。

四、单因素方差分析的应用条件

各样本是相互独立的随机样本，均服从正态分布；
各样本的总体方差相等，即具有方差齐性；对方差齐性的判断通常采用方差齐性检验的方差，检验多个样本所代表的的总体方差是否不等，可以采用的方法有Bartlett卡方检验和Levene检验。

五、单因素方差分析步骤

1、各组进行正态性检验；如果不满足正态分布，考虑非参数检验。

2、进行方差齐性检验；如果不满足方差齐性，考虑非参数检验。

3、单因素方差分析；

4、多个样本均数的两两比较；多个样本均数的两两比较不可以采用两样本均数t检验进行分析，不然会增加I类错误的概率。多个样本均数多重比较的方法较多，常用的方法有SNK法和Dunnett法；

六、R代码实现

####单因素方差分析####
library(car)
library(multcomp)#加载包，里面有数据
library(agricolae)
attach(cholesterol)#正态性检验
shapiro.test(response[trt=="1time"])
shapiro.test(response[trt=="2times"])
shapiro.test(response[trt=="4times"])
shapiro.test(response[trt=="drugD"])
shapiro.test(response[trt=="drugE"])
#方差齐性检验
leveneTest(response,trt)#单因素方差分析
fit<-aov(response~trt)
summary(fit)#绘图
library(gplots)#绘制各组均值及其置信区间的图形
plotmeans(response~trt,xlab="treatment",ylab="Response",main="Mean Plot with 95%CI")
detach(cholesterol)####方差分析多重比较####
#Turkey法
TukeyHSD(fit)
par(las=1) #轴文字方向
par(mar=c(5,8,4,2))#图形边界距离
plot(TukeyHSD(fit))#Turkey法和Dunnett法
par(mar=c(5,4,8,2))
tuk<-glht(fit,linfct=mcp(trt="Tukey"),alternative="two.side")
summary(tuk)
plot(cld(tuk,level=.05),col="lightgrey")dunnett<-glht(fit,linfct=mcp(trt="Dunnett"),alternative="two.side")#
summary(dunnett)
#SNK法
library(agricolae)
out<-SNK.test(fit,"trt",console=T)
print(SNK.test(fit,"trt",group=FALSE))

结果：

> ####单因素方差分析####
> library(car)
> library(multcomp)#加载包，里面有数据
> library(agricolae)
> attach(cholesterol)
The following objects are masked from cholesterol (pos = 5):response, trt> table(trt)
trt1time 2times 4times  drugD  drugE 10     10     10     10     10
> #正态性检验
> shapiro.test(response[trt=="1time"])Shapiro-Wilk normality testdata:  response[trt == "1time"]
W = 0.93063, p-value = 0.4541> shapiro.test(response[trt=="2times"])Shapiro-Wilk normality testdata:  response[trt == "2times"]
W = 0.9432, p-value = 0.5892> shapiro.test(response[trt=="4times"])Shapiro-Wilk normality testdata:  response[trt == "4times"]
W = 0.95487, p-value = 0.7262> shapiro.test(response[trt=="drugD"])Shapiro-Wilk normality testdata:  response[trt == "drugD"]
W = 0.93406, p-value = 0.489> shapiro.test(response[trt=="drugE"])Shapiro-Wilk normality testdata:  response[trt == "drugE"]
W = 0.98085, p-value = 0.9696经过正态性检验，五组p>0.05,均符合正态分布。> #方差齐性检验
> leveneTest(response,trt)
Levene's Test for Homogeneity of Variance (center = median)Df F value Pr(>F)
group  4  0.0755 0.989345
经过方差齐性检验，P=0.9893>0.05,方差齐。可进行方差分析。> fit<-aov(response~trt)
> summary(fit)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
> 方差分析结果显示，p<0.05,各组间具有统计学差异。> library(gplots)#绘制各组均值及其置信区间的图形
> plotmeans(response~trt,xlab="treatment",ylab="Response",main="Mean Plot with 95%CI")
> detach(cholesterol)> ####方差分析多重比较####
> TukeyHSD(fit)Tukey multiple comparisons of means95% family-wise confidence levelFit: aov(formula = response ~ trt)$trtdiff        lwr       upr     p adj
2times-1time   3.44300 -0.6582817  7.544282 0.1380949
4times-1time   6.59281  2.4915283 10.694092 0.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.251092 0.2050382
drugD-2times   6.13620  2.0349183 10.237482 0.0009611
drugE-2times  11.72255  7.6212683 15.823832 0.0000000
drugD-4times   2.98639 -1.1148917  7.087672 0.2512446
drugE-4times   8.57274  4.4714583 12.674022 0.0000037
drugE-drugD    5.58635  1.4850683  9.687632 0.0030633> par(las=1) #轴文字方向
> par(mar=c(5,8,4,2))#图形边界距离
> plot(TukeyHSD(fit))
>
> library(multcomp)
> par(mar=c(5,4,8,2))
> tuk<-glht(fit,linfct=mcp(trt="Tukey"),alternative="two.side")
> summary(tuk)Simultaneous Tests for General Linear HypothesesMultiple Comparisons of Means: Tukey ContrastsFit: aov(formula = response ~ trt)Linear Hypotheses:Estimate Std. Error t value Pr(>|t|)
2times - 1time == 0     3.443      1.443   2.385  0.13810
4times - 1time == 0     6.593      1.443   4.568  < 0.001 ***
drugD - 1time == 0      9.579      1.443   6.637  < 0.001 ***
drugE - 1time == 0     15.166      1.443  10.507  < 0.001 ***
4times - 2times == 0    3.150      1.443   2.182  0.20505
drugD - 2times == 0     6.136      1.443   4.251  < 0.001 ***
drugE - 2times == 0    11.723      1.443   8.122  < 0.001 ***
drugD - 4times == 0     2.986      1.443   2.069  0.25125
drugE - 4times == 0     8.573      1.443   5.939  < 0.001 ***
drugE - drugD == 0      5.586      1.443   3.870  0.00307 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)> plot(cld(tuk,level=.05),col="lightgrey")
>
> dunnett<-glht(fit,linfct=mcp(trt="Dunnett"),alternative="two.side")#
> summary(dunnett)Simultaneous Tests for General Linear HypothesesMultiple Comparisons of Means: Dunnett ContrastsFit: aov(formula = response ~ trt)Linear Hypotheses:Estimate Std. Error t value Pr(>|t|)
2times - 1time == 0    3.443      1.443   2.385     0.07 .
4times - 1time == 0    6.593      1.443   4.568   <0.001 ***
drugD - 1time == 0     9.579      1.443   6.637   <0.001 ***
drugE - 1time == 0    15.166      1.443  10.507   <0.001 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)> dunnettGeneral Linear HypothesesMultiple Comparisons of Means: Dunnett ContrastsLinear Hypotheses:Estimate
2times - 1time == 0    3.443
4times - 1time == 0    6.593
drugD - 1time == 0     9.579
drugE - 1time == 0    15.166>
> library(agricolae)
> out<-SNK.test(fit,"trt",console=T)Study: fit ~ "trt"Student Newman Keuls Test
for response Mean Square Error:  10.41668 trt,  meansresponse      std  r     Min     Max
1time   5.78197 2.878113 10  2.3039 10.3868
2times  9.22497 3.483054 10  3.5054 15.8258
4times 12.37478 2.923119 10  8.0534 18.1794
drugD  15.36117 3.454636 10  8.8194 19.9793
drugE  20.94752 3.345003 10 15.7707 27.2445Alpha: 0.05 ; DF Error: 45 Critical Range2        3        4        5
2.907109 3.498187 3.850494 4.101282 Means with the same letter are not significantly different.response groups
drugE  20.94752      a
drugD  15.36117      b
4times 12.37478      c
2times  9.22497      d
1time   5.78197      e> print(SNK.test(fit,"trt",group=FALSE))
$statisticsMSerror Df     Mean       CV10.41668 45 12.73808 25.33731$parameterstest name.t ntr alphaSNK    trt   5  0.05$snkTable CriticalRange
2 2.848372      2.907109
3 3.427507      3.498187
4 3.772697      3.850494
5 4.018417      4.101282$meansresponse      std  r     Min     Max       Q25      Q50       Q75
1time   5.78197 2.878113 10  2.3039 10.3868  3.260975  5.41510  7.672825
2times  9.22497 3.483054 10  3.5054 15.8258  7.913800  8.61465 10.106325
4times 12.37478 2.923119 10  8.0534 18.1794 10.554875 12.60540 13.949850
drugD  15.36117 3.454636 10  8.8194 19.9793 13.873475 16.21975 17.621125
drugE  20.94752 3.345003 10 15.7707 27.2445 18.825625 21.01590 22.561875$comparisondifference pvalue signif.        LCL          UCL
1time - 2times    -3.44300 0.0213       *  -6.350109  -0.53589081
1time - 4times    -6.59281 0.0001     *** -10.090997  -3.09462347
1time - drugD     -9.57920 0.0000     *** -13.429694  -5.72870568
1time - drugE    -15.16555 0.0000     *** -19.266832 -11.06426832
2times - 4times   -3.14981 0.0344       *  -6.056919  -0.24270081
2times - drugD    -6.13620 0.0003     ***  -9.634387  -2.63801347
2times - drugE   -11.72255 0.0000     *** -15.573044  -7.87205568
4times - drugD    -2.98639 0.0443       *  -5.893499  -0.07928081
4times - drugE    -8.57274 0.0000     *** -12.070927  -5.07455347
drugD - drugE     -5.58635 0.0003     ***  -8.493459  -2.67924081$groups
NULLattr(,"class")
[1] "group"

R语言单因素方差分析及两两比较相关推荐

R语言单因素方差分析与协方差分析
R语言单因素方差分析与协方差分析条件: 各个样本是相互独立的随机: 各个样本来自正态总体: 具有方差齐性: 用途: 检验两个或多样本均数间的差异有无统计学意义:注:本均数的比较可以采用 t检验或 F ...
R语言单因素方差分析（One-Way ANOVA）实战：探索性数据分析（EDA）、单因素方差分析模型结果解读（检查模型假设）、分析不同分组的差异TukeyHSD、单因素方差分析的结果总结
R语言单因素方差分析(One-Way ANOVA)实战:探索性数据分析(EDA).单因素方差分析模型结果解读(检查模型假设).分析不同分组的差异TukeyHSD.单因素方差分析的结果总结目录 R语言 ...
R语言单因素方差分析（附代码）
单因素方差分析概念不给大家赘述了,可自行查阅相关书籍. 直接上题目: 题目 pseudomonas aeruginosa菌在样本A中的丰度为:0.014,0.015,0.017:在样本B中的丰度为: ...
R语言效用分析 ( 效能分析、Power analysis)确定样本量、假设检验与两类错误、pwr包进行效用分析 ( 效能分析、Power analysis)的常用函数列表
R语言效用分析 ( 效能分析.Power analysis)确定样本量.假设检验与两类错误.pwr包进行效用分析 ( 效能分析.Power analysis)的常用函数列表目录
R语言生存分析寿命表（life table）实战案例：比较两种药物治疗感染患者的生存时间
R语言生存分析寿命表(life table)实战案例:比较两种药物治疗感染患者的生存时间目录
R语言生存分析COX回归分析实战:两种治疗方法发生肾功能损害的情况
R语言生存分析COX回归分析实战:两种治疗方法发生肾功能损害的情况目录
R语言ggplot2可视化左对齐两个可视化图像的画图区域（边缘）实战
R语言ggplot2可视化左对齐两个可视化图像的画图区域(边缘)实战目录 R语言ggplot2可视化左对齐两个可视化图像的画图区域(边缘)实战
R语言outer函数应用于两个数组实战
R语言outer函数应用于两个数组实战目录 R语言outer函数应用于两个数组实战 #基本语法 #outer函数应用于向量和单个值
R语言使用order函数按照两个数据列的值排序data.table数据（主变量降序排序、次变量升序排序）
R语言使用order函数按照两个数据列的值排序data.table数据(主变量降序排序.次变量升序排序) 目录

R语言单因素方差分析及两两比较

R语言单因素方差分析及两两比较相关推荐

最新文章

热门文章