2009-06-24

R语言：常用统计检验方法

R已经成为当前国际学术界最流行的统计和绘图软件之一，该语言较为简单易学，统计分析功能强大，且具有很强的绘图功能，能够绘制学术出版要求的多种图表.R语言在生物信息学，进化生物学、生态学与环境、经济学、语言学等领域有着极为广泛的应用。

R软件是跨平台的，可以在Linux， MacOs， Windows等多种系统上运行。针对每个研究方向，有大量的科研人员编写了相关的程序包，可以导入到基本的程序平台上运行。现有的程序包已经超过了1800个，并且还在增加中。

不仅如此，R是完全免费的，而且全部代码是公开的。

读者可以到 https://cran.r-project.org/ 下载windows版的R软件，安装程序仅为30M。

学习并掌握R语言，对于需要用到统计学的研究人员和学生都是非常必要的。

这里选取了R语言中若干操作实例，所有的命令行均可以在R中运行，并得到结果。

正态总体均值的假设检验

t检验

单个总体

例一
某种元件的寿命X（小时），服从正态分布，N（mu,sigma^2），其中mu,sigma^2均未知，16只元件的寿命如下：问是否有理由认为元件的平均寿命大于255小时。
命令：

1
2
3

X <- c(159, 280, 101, 212, 224, 379, 179, 264,
222, 362, 168, 250, 149, 260, 485, 170)
t.test(X, alternative = "greater", mu = 225)

两个总体

例二
X为旧炼钢炉出炉率，Y为新炼钢炉出炉率，问新的操作能否提高出炉率
命令：

1
2
3

X <- c(78.1,72.4,76.2,74.3,77.4,78.4,76.0,75.5,76.7,77.3)
Y <- c(79.1,81.0,77.3,79.1,80.0,79.1,79.1,77.3,80.2,82.1)
t.test(X, Y, var.equal=TRUE, alternative = "less")

配对t检验

例三
对每个高炉进行配对t检验
命令：

1
2
3

X <- c(78.1,72.4,76.2,74.3,77.4,78.4,76.0,75.5,76.7,77.3)
Y <- c(79.1,81.0,77.3,79.1,80.0,79.1,79.1,77.3,80.2,82.1)
t.test(X - Y, alternative = "less")

正态总体方差的假设检验

例四
从小学5年级男生中抽取20名，测量其身高（厘米）如下：
问，在0.05显著性水平下，平均值是否等于149，sigma^2 是否等于 75

命令：

1
2
3

X <- c(136, 144, 143, 157, 137, 159, 135, 158, 147, 165,
158, 142, 159, 150, 156, 152, 140, 149, 148, 155)
var.test(X,Y)

例五
对炼钢炉的数据进行分析，命令：

1
2
3

X <- c(78.1,72.4,76.2,74.3,77.4,78.4,76.0,75.5,76.7,77.3)
Y <- c(79.1,81.0,77.3,79.1,80.0,79.1,79.1,77.3,80.2,82.1)
var.test(X,Y)

二项分布的总体检验

例六有一批蔬菜种子的平均发芽率为P=0.85,现在随机抽取500粒，用种衣剂进行浸种处理，结果有445粒发芽，问种衣剂有无效果。
命令：

1	binom.test(445,500,p=0.85)

例七按照以往经验，新生儿染色体异常率一般为1%，某医院观察了当地400名新生儿，有一例染色体异常，问该地区新生儿染色体是否低于一般水平？
命令：

1	binom.test(1,400,p=0.01,alternative="less")

非参数检验

数据是否正态分布的Neyman-Pearson 拟合优度检验-chisq

例八 5种品牌啤酒爱好者的人数如下：

啤酒品牌	人数
A	210
B	312
C	170
D	85
E	223

问不同品牌啤酒爱好者人数之间有没有差异？

命令：

1 2	X <- c(210, 312, 170, 85, 223) chisq.test(X)

例九
检验学生成绩是否符合正态分布
命令：

X <- c(25,45,50,54,55,61,64,68,72,75,75,
78,79,81,83,84,84,84,85,86,86,86,
87,89,89,89,90,91,91,92,100)

A <- table(cut(X, br=c(0,69,79,89,100)))
p <- pnorm(c(70,80,90,100), mean(X), sd(X))
p <- c(p[1], p[2]-p[1], p[3]-p[2], 1-p[3])
chisq.test(A,p=p)

说明：

cut 将变量区域划分为若干区间
table 计算因子合并后的个数
均值之间有无显著区别

大麦的杂交后代芒性状的比例无芒：长芒：短芒=9：3：4,而实际观测值为335：125：160 ,检验观测值是否符合理论假设？
命令：

1	chisq.test(c(335, 125, 160), p=c(9,3,4)/16)

例十现有42个数据，分别表示某一时间段内电话总机接到呼叫的次数，
接到呼叫的次数 0 1 2 3 4 5 6
出现的频率 7 10 12 8 3 2 0

问：某个时间段内接到的呼叫次数是否符合Possion分布？

命令：

x <- 0:6
y <- c(7,10,12,8,3,2,0)
mean <- mean(rep(x,y))
q <- ppois(x,mean)
n <- length(y)
p[1] <- q[1]
p[n] <- 1-q[n-1]
for(i in 2:(n-1)) {
    p <- q-q
}
chisq.test(y, p=p)
Z <- c(7, 10, 12, 8)
n <- length(Z); p <- p[1:n-1]; p[n] <- 1-q[n-1]
chisq.test(Z, p=p)

按：处理方法有点儿复杂，貌似不是最好的方法

理论分布依赖于若干未知参数时

Kolmogorov-Smirnov 检验 ks.test()

例一对一台设备进行寿命检验，记录十次无故障操作时间，并按从小到大的次序排列如下，用ks检验方法检验此设备无故障工作时间是否符合rambda=1/1500的指数分布

命令：

1 2	X <- c(420, 500, 920, 1380, 1510, 1650, 1760, 2100, 2300, 2350) ks.test(X, "pexp", 1/1500)

例二假设从分布函数F(x)和G(x)的总体中分别随机抽取25个和20个观察值样本，检验F(x)和G(x)是否相同。

命令

X <- c(
0.61,0.29,0.06,0.59,-1.73,-0.74,0.51,-0.56,0.39,
1.64,0.05,-0.06,0.64,-0.82,0.37,1.77,1.09,-1.28,
2.36,1.31,1.05,-0.32,-0.40,1.06,-2.47)

Y <- c(2.20,1.66,1.38,0.20,0.36,0.00,0.96,1.56,0.44,
1.50,-0.30,0.66,2.31,3.29,-0.27,-0.37,0.38,0.70,
0.52,-0.71)

ks.test(X, Y)

ks多样本检验的局限性，只用在理论分布为一维连续分布，且分布完全已知的情形。ks检验可用的情况下，功效一般优于Pearson chisq检验

列联表（contingerncy table）的独立性检验

Pearson chisquare 进行独立性检验

例三为了研究吸烟是否与肺癌有关，对63位患者及43名非肺癌患者调查了其中的吸烟人数，得到2*2列联表

数据	肺癌	健康	合计
吸烟	60	32	92
不吸烟	3	11	14
合计	63	43	106

命令

x <- c(60, 3, 32, 11)
dim(x) <- c(2,2)
chisq.test(x,correct = FALSE) # 不带连续校正的情况
chisq.test(x) # 带连续校正的情况

例四在一次社会调查中，以问卷方式调查了901人的年收入，及其对工作的满意程度，其中年收入A分为四档：小于6000元，6000-15000元，15000元至25000元，超过25000元。对工作的满意程度B 分为很不满意，较不满意，基本满意和很满意四档，结果如下

	很不满意	较不满意	基本满意	很满意	合计
< 6000	20	24	80	82	206
6000 ~15000	22	38	104	125	289
15000~25000	13	28	81	113	235
> 25000	7	18	54	92	171
合计	62	108	319	412	901

命令如下

x <- c(20,24,80,82,22,38,104,125,
13,28,81,113,7,18,54,92)

dim(x) <- c(4,4)
chisq.test(x)

Fisher 精确的独立检验

试用条件样本数小于4

例五某医师研究乙肝免疫球蛋白防止子宫内胎儿感染HBV的效果，将33例HBsAg阳性孕妇随机分为预防注射组和对照组，结果由下表所示，两组新生儿HBV总体感染率有无差别

组别	阳性	阴性	合计	感染率
预防注射组	4	18	22	18.8
对照组	5	6	11	45.5

命令如下

1
2
3

x <- c(4,5,18,6)
dim(x) <- c(2,2)
fisher.test(x)

对前面提到的肺癌进行检验

1
2
3

x <- c(60, 3, 32, 11)
dim(x) <- c(2,2)
fisher.test(x)

McNemar检验

McNemar检验不是独立性检验，但是是关于列连表的检验

例六甲、乙两种方法检测细菌的结果如下，问：

	甲	乙	合计
+	49	25	74
-	21	107	128
合计	70	132	202

命令

1
2
3

X  <-  c(49, 21, 25, 107);
dim(X)  <-  c(2,2)
mcnemar.test(X,correct=FALSE)

符号检验

1 假设一个样本是否来自某个总体

例七联合国人员在世界上66个大城市的生活花费指数（以纽约1996年12月为100），按照从小到大的次序排列如下，其中北京的指数为99。假设这个样本是从世界大城市中随机抽样得到的。用符号检验分析，北京是在中位数之上，还是中位数之下。

X <- c(
    66,75,78,80,81,81,82,83,83,83,83,
    84,85,85,86,86,86,86,87,87,88,88,
    88,88,88,89,89,89,89,90,90,91,91,
    91,91,92,93,93,96,96,96,97,99,100,
    101,102,103,103,104,104,104,105,106,109,109,
    110,110,110,111,113,115,116,117,118,155,192)

binom.test(sum(X>99),length(X), al="l")

2 用成对样本检验两总体间是否有差异

例八两种不同饲料，对猪增重情况如下，分析两种饲料养猪有无差异

命令

1
2
3

x <- c(25,30,28,23,27,35,30,28,32,29,30,30,31,16)
y <- c(19,32,21,19,25,31,31,26,30,25,28,31,25,25)
binom.test(sum(x<y), length(x))

例九某饮料店为调查了顾客对饮料的爱好情况，某日随机调查了13为顾客，喜欢奶茶超过咖啡用-表示，喜欢咖啡超过奶茶用+表示，两者都喜欢用0表示，结果如下，分析顾客是更喜欢咖啡开始奶茶。

顾客编号	1	2	3	4	5	6	7	8	9	10	11	12	13
喜欢咖啡	1		1	1	1	0	1		1	1	1		1
喜欢奶茶		1						1				1

1	binom.test(3,12,p=1/2, al="l", conf.level = 0.90)

Spearman秩相关检验

例十一项有六人参加表演的竞赛，有两人进行评定，评定结果用如表所示，试用Spearman秩相关检验方法检验这两个评定员对于等级评定有无相关性

选手编号	1	2	3	4	5	6
甲的打分	4	2	2	4	5	6
乙的打分	5	3	4	3	2	5

1
2
3

x <- c(4,2,2,4,5,6)
y <- c(5,3,4,3,2,5)
cor.test(x, y, method = "spearman")

Kendall相关检验

例十一某幼儿园对9对双胞胎的智力进行测验，并按照百分制打分，试用Kendall相关检验方法检验双胞胎的智力是否相关。

1	2	3	4	5	6	7	8	9
86	77	68	91	70	71	85	87	63
88	76	64	96	65	80	81	72	60

1
2
3

X <- c(86, 77, 68, 91, 70, 71, 85, 87, 63)
Y <- c(88, 76, 64, 96, 65, 80, 81, 72, 60)
cor.test(X, Y, method = "kendall")

Wilcoxon秩检验—— 考虑了样本观察值月总体中位数的差。

1 对于来自同一个总体样本的检验

例十二某电池厂称其生产的某种电池，中位数为140安培小时，现随机从其新生产的电池中抽取20个，检验其寿命：137.0, 140.0, 138.3, 139.0, 144.3, 139.1, 141.7, 137.3, 133.5, 138.2, 141.1, 139.2, 136.5, 136.5, 135.6, 138.0, 140.9, 140.6, 136.3, 134.1

用Wilcoxon符号检验分析该厂生产的电池是否符合标准

X <- c(137.0, 140.0, 138.3, 139.0, 144.3, 139.1, 141.7, 137.3, 133.5,
       138.2, 141.1, 139.2, 136.5, 136.5, 135.6, 138.0, 140.9, 140.6, 136.3, 134.1)

wilcox.test(X,mu=140, alternative="less",exact=FALSE, correct=FALSE, conf.int=TRUE)

该方法也可用于成对样本的检验

例十三为检验某种新肥料，将现有麦地分为十块，再将每一块分为两部分，一半施普通肥料，一半儿施新肥料，用Wilcoxon符号检验法检验新复合肥能否显著提高小麦产量。

1	2	3	4	5	6	7	8	9	10
459	367	303	392	310	342	421	446	430	412
414	306	321	443	281	301	353	391	405	390

x <- c(459, 367, 303, 392, 310, 342, 421, 446, 430, 412)
y <- c(414, 306, 321, 443, 281, 301, 353, 391, 405, 390)
wilcox.test(x, y, alternative = "greater", paired = TRUE)
wilcox.test(x-y, alternative = "greater")
binom.test(sum(x>y), length(x), alternative = "greater")

非成对样本的秩次和检验

Wilcoxon-Mann-Whitney 统计量 U

例十四测量了10名不同作业组的工人血铅含量，分析两组之间是否有差别。

非铅作业组 24 26 29 34 43 58 63 72 87 101

含铅作业组 82 87 97 121 164 208 213

x <- c(24, 26, 29, 34, 43, 58, 63, 72, 87, 101)
y <- c(82, 87, 97, 121, 164, 208, 213)
wilcox.test(x,y,alternative="less",exact=FALSE,correct=FALSE)
wilcox.test(x, y, alternative="less", exact=FALSE)

例十五学生数学能力排序: 新方法 3 5 7 9 10, 原方法 1 2 4 6 8。问新旧方法之间是否有差别。

1
2
3

x <- c(3, 5, 7, 9, 10);
y <- c(1, 2, 4, 6, 8)
wilcox.test(x, y, alternative="greater")

例十六检验一种药物对于慢性支气管炎有没有效果，抽取了216个病例，治疗效果。分析该药物对两种慢性支气管炎的治疗效果是否相同。

   控制 显效 进步 无效

单纯型 62 41 14 11
喘息型 20 37 16 15

1
2
3

x <- rep(1:4, c(62, 41, 14,11));
y <- rep(1:4, c(20, 37, 16, 15))
wilcox.test(x, y, exact=FALSE)

参考资料

薛毅陈立萍《统计建模与Ｒ软件》清华大学出版社 2006

1	2	3	4	5	6	7	8	9	10
459	367	303	392	310	342	421	446	430	412
414	306	321	443	281	301	353	391	405	390

1	2	3	4	5	6	7	8	9	10
459	367	303	392	310	342	421	446	430	412
414	306	321	443	281	301	353	391	405	390

正态总体均值的假设检验

t检验

单个总体

两个总体

配对t检验

正态总体方差的假设检验

二项分布的总体检验

数据是否正态分布的Neyman-Pearson 拟合优度检验-chisq

问：某个时间段内接到的呼叫次数是否符合Possion分布？

理论分布依赖于若干未知参数时

Kolmogorov-Smirnov 检验 ks.test()

命令

列联表（contingerncy table）的独立性检验

Pearson chisquare 进行独立性检验

命令

命令如下

Fisher 精确的独立检验

试用条件 样本数小于4

对前面提到的肺癌进行检验

McNemar检验

符号检验

1 假设一个样本是否来自某个总体

2 用成对样本检验两总体间是否有差异

命令

Spearman秩相关检验

Kendall相关检验

Wilcoxon秩检验—— 考虑了样本观察值月总体中位数的差。

1 对于来自同一个总体样本的检验

用Wilcoxon符号检验分析该厂生产的电池是否符合标准

该方法也可用于成对样本的检验

非成对样本的秩次和检验

Wilcoxon-Mann-Whitney 统计量 U

参考资料

试用条件样本数小于4

1	2	3	4	5	6	7	8	9	10
459	367	303	392	310	342	421	446	430	412
414	306	321	443	281	301	353	391	405	390