文章目录
示例数据:
MASS
包中的
birthwt
数据集。
首先将数据集中的分类变量因子化,具体参考
这里。
独立样本t检验
方差齐性检验 (两组):var.test()
语法:var.test(连续变量名~分组变量名,data=数据框名)
> var.test(bwt~smoke,data=birthwt)
F test to compare two variances
data: bwt by smoke
F = 1.3019, num df = 114, denom df = 73, p-value = 0.2254 # P值>0.05,表示可认为方差齐性
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.8486407 1.9589574
sample estimates:
ratio of variances
1.301927
独立样本t检验:t.test()
语法:t.test(连续变量名,分组变量名,var.equal=[TRUE/FALSE],data=数据框名)
根据方差齐性检验结果,声明var.equal
值 (默认为FALSE)
> t.test(bwt~smoke,data=birthwt,var.equal=TRUE)
Two Sample t-test
data: bwt by smoke
t = 2.6529, df = 187, p-value = 0.008667 # 读取P值
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
72.75612 494.79735 # 组0减去组1的差的95%CI
sample estimates:
mean in group 0 mean in group 1 # 组0和组1所在总体均值的估计值
3055.696 2771.919
若要返回非95%CI (如99%CI),需要额外声明参数conf.level
。
非独立样本t检验
即配对样本t检验。语法类似于独立样本t检验,仅需声明参数paired=TRUE
即可:
t.test(连续变量名,分组变量名,var.equal=[TRUE/FALSE],data=数据框名,paired=TRUE)
单因素方差分析 (ANOVA)
首先,复习单因素ANOVA应用条件:
- 需比较的组≥3个
- 可认定各组数据是从正态总体中独立抽样得到的 (正态性检验)
- 各组数据方差齐性
正态性检验:tapply(,shapiro.test)
语法:tapply(数据框名$连续变量名,数据框名$分组变量名,shapiro.test)
> tapply(birthwt$bwt,birthwt$race,shapiro.test)
$white
Shapiro-Wilk normality test
data: X[[i]]
W = 0.98727, p-value = 0.4861
$black
Shapiro-Wilk normality test
data: X[[i]]
W = 0.97696, p-value = 0.8038
$others
Shapiro-Wilk normality test
data: X[[i]]
W = 0.97537, p-value = 0.2046
方差齐性检验 (多组)
Bartlett检验 | Levene检验 |
---|---|
参数检验 | 非参数检验 |
对数据的正态性敏感 | 适用范围更广 |
Bartlett检验:bartlett.test()
语法:bartlett.test(连续变量名~分组变量名,data=数据框名)
> bartlett.test(bwt~race,data=birthwt)
Bartlett test of homogeneity of variances
data: bwt by race
Bartlett's K-squared = 0.65952, df = 2, p-value = 0.7191
Levene检验:leveneTest()
该函数位于car
包下。
语法:leveneTest(连续变量名~分组变量名,data=数据框名)
> leveneTest(bwt~race,data=birthwt)
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 2 0.4684 0.6267
186
建立单因素ANOVA模型:aov()
语法:aov模型名<-aov(连续变量名~分组变量名,data=数据框名)
使用summary()
函数查看单因素ANOVA结果。
race.aov<-aov(birthwt$bwt~birthwt$race)
> summary(race.aov)
Df Sum Sq Mean Sq F value Pr(>F)
birthwt$race 2 5015725 2507863 4.913 0.00834 **
Residuals 186 94953931 510505
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
事后检验
方差分析得出的P<0.05意味着不能认为几组变量均值全等。此后还要进一步探索变量分组两两之间的均值关系,即事后检验 (post-hoc test)。
基于不同的对第I类错误的控制方法,有多种事后检验的方法可以选择。
TukeyHSD检验:TukeyHSD()
语法:TukeyHSD(aov模型名)
> TukeyHSD(race.aov)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = birthwt$bwt ~ birthwt$race)
$`birthwt$race`
diff lwr upr p adj
black-white -383.02644 -756.2363 -9.816581 0.0428037
others-white -297.43517 -566.1652 -28.705095 0.0260124
others-black 85.59127 -304.4521 475.634630 0.8624372
其他校正方法:pairwise.t.test()
语法:pairwise.t.test(数据框名$连续变量名,数据框名$分组变量名,p.adjust.method=声明的校正方法)
可选用的校正方法包括:“holm”, “hochberg”, “hommel”, “bonferroni”, “BH”, “BY”, “fdr”, “none”
。以Bonferroni校正为例举例如下:
> pairwise.t.test(birthwt$bwt,birthwt$race,p.adjust.method = "bonferroni")
Pairwise comparisons using t tests with pooled SD
data: birthwt$bwt and birthwt$race
white black
black 0.049 -
others 0.029 1.000
P value adjustment method: bonferroni
组间差异非参数检验
类似于参数检验,不再赘述。
统计目标 | 非参数检验 | 语法 |
---|---|---|
两组独立样本 | Mann-Whitney-Wilcoxin检验 | wilcox.test(连续变量名,分组变量名,data=数据框名) |
两组配对样本 | Wilcoxin符号秩和检验 | wilcox.test(变量1,变量2,paired=TRUE) |
多组独立样本 | Kruskal-Wallis检验 | kruskal.test(连续变量名,分组变量名,data=数据框名) |
两组对比
独立:wilcox.test()
> wilcox.test(bwt~smoke,data=birthwt)
Wilcoxon rank sum test with continuity correction
data: bwt by smoke
W = 5249.5, p-value = 0.006768
alternative hypothesis: true location shift is not equal to 0
配对:wilcox.test(,paired=TRUE)
> x<-c(10,12,15,18,20,25,18,22)
> y<-c(20,15,20,18,22,23,17,22)
> wilcox.test(x,y,paired = TRUE)
Wilcoxon signed rank test with continuity correction
data: x and y
V = 3.5, p-value = 0.1718
alternative hypothesis: true location shift is not equal to 0
Warning messages: # 当配对样本中存在相同的结果时,会出现如此的警告信息
1: In wilcox.test.default(x, y, paired = TRUE) : 无法精確計算带连结的p值
2: In wilcox.test.default(x, y, paired = TRUE) : 有0时无法計算精確的p值
多组:kruskal.test()
> kruskal.test(bwt~race,data=birthwt)
Kruskal-Wallis rank sum test
data: bwt by race
Kruskal-Wallis chi-squared = 8.5199, df = 2, p-value =
0.01412
事后检验
可以使用Wilcoxin符号秩和检验 (控制第1类错误概率)。
也可以使用PMCMRplus
包中的bwsAllPairsTest()
进行BWS检验。
语法:bwsAllPairsTest(连续变量名,分组变量名,data=数据框名)
> bwsAllPairsTest(bwt~race,data=birthwt)
Pairwise comparisons using BWS All-Pairs Test
data: bwt by race
white black
black 0.027 -
others 0.030 0.535
P value adjustment method: holm
alternative hypothesis: two.sided