8.1 单比例

对单比例的检验基于二项分布。在大样本的情况下，可以用一个均值为\(Np\)，方差为\(Np(1-p)\)的正态分布来很好的近似。有一个经验结论，在“成功”与“失败”的期望次数都大于5的时候能得到较好的近似。 令\(x\)表示“成功”的次数，则对假设\(p=p_0\)的检验可以基于下式： \[ u=\frac{x-Np_0}{\sqrt{Np_0(1-p_0)}} \] 这个统计量近似服从一个均值为0，标准差为1的正态分布，或者说\(u^2\)近似服从一个自由度为1的卡方分布。 下面开始讲例子。215名病人的39名被观测到患有哮喘，然后对“随机病人”患有哮喘的概率是0.15这个假设做检验。采用函数prop.test。

#这个函数的3个参数分别是，正观测数、总数以及想对其做检验的概率参数
#最后这个参数默认值为0.5
prop.test(39,215,0.15)

## 
##  1-sample proportions test with continuity correction
## 
## data:  39 out of 215, null probability 0.15
## X-squared = 1.425, df = 1, p-value = 0.2326
## alternative hypothesis: true p is not equal to 0.15
## 95 percent confidence interval:
##  0.1335937 0.2408799
## sample estimates:
##         p 
## 0.1813953

也可以用函数binom.test在二项分布下做检验，这时能得到精确地检验概率，所以一般它比上一个函数更受欢迎。不过除了单比例检验，函数prop.test还能做其他事情。为了得到这个\(p\)值，先计算出\(x\)取每一个可能值的点概率，然后再将观测到的小于等于\(x\)的概率加起来。

binom.test(39,215,.15)

## 
##  Exact binomial test
## 
## data:  39 and 215
## number of successes = 39, number of trials = 215, p-value = 0.2135
## alternative hypothesis: true probability of success is not equal to 0.15
## 95 percent confidence interval:
##  0.1322842 0.2395223
## sample estimates:
## probability of success 
##              0.1813953

8.2 两个独立的比例

函数prop.test也能用于比较两个或多个比例。为达到这个目的，参数应该是两个向量，其中第一个表示正观测数，第二个表示每组总数。 以下是例子：

lewitt.machin.success <- c(9,4)
lewitt.machin.total <- c(12,13)
prop.test(lewitt.machin.success,lewitt.machin.total)

## 
##  2-sample test for equality of proportions with continuity
##  correction
## 
## data:  lewitt.machin.success out of lewitt.machin.total
## X-squared = 3.2793, df = 1, p-value = 0.07016
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  0.01151032 0.87310506
## sample estimates:
##    prop 1    prop 2 
## 0.7500000 0.3076923

这里给出的是比例之差的置信区间。 也可以不使用Yates连续型修正来计算这个检验，即加入参数correct=F。连续性修正一定程度上让所得置信区间变得更宽，不过注意到这里仍然不包含0。而双边假设检验显示两组之间没有显著差异，所以这里置信区间与其相矛盾。这是因为使用了不同的近似方法导致的，这对于稀疏如此的表尤为重要。 如果希望确认至少\(p\)值是正确的，可以用Fisher精确检验。用上一节的数据来展示它。这个检验在给定行和列的边际值的情况下计算\(2\times 2\)表格的条件分布。这可能难以想象，不过可以这样来看：拿出13个白球和12个黑球（分别是成功与失败），然后不重复地从抽样中得到两个样本量分别为12和13的小组。第一组的白球数量显然足以确定整个表格，而它的分布可以完全转化为一个纯粹的组合问题。这个分布被称为超几何分布。 相关的检验是fisher.test，它要求输入的数据时矩阵形式的，如下：

matrix(c(9,4,3,9),2)

##      [,1] [,2]
## [1,]    9    3
## [2,]    4    9

lewitt.machin <- matrix(c(9,4,3,9),2)
fisher.test(lewitt.machin)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  lewitt.machin
## p-value = 0.04718
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   0.9006803 57.2549701
## sample estimates:
## odds ratio 
##   6.180528

注意表格的第二列应该是负结果的次数，不是观测值的总数。 还需注意的是，这里给出了比值比（odds ratio）的置信区间，即是\((p_1/(1-p_1))/(p_2/(1-p_2))\)。可以发现，若\(p_1\neq p_2\)，那么表格的条件分布只依赖于比值比，所以这是一个用于衡量Fisher检验中相关程度的自然指标。这个检验的精确分布在比值比不为1的时候可以精确地求出，不过这里与binom.test一样，双边95%置信区间是由两个单边97.5%置信区间并起来的。这导致它的结果与prop.test的不一致：检验的结果刚好显著，但是比值比的置信区间包含1在内。 和fisher.test一样，在chisq.test中的标准\(\chi^2\)检验需要矩阵类型的数据。对于一个\(2\times 2\)表格来说，这个检验与prop.test的结果是完全一样的。

chisq.test(lewitt.machin)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  lewitt.machin
## X-squared = 3.2793, df = 1, p-value = 0.07016

8.3 k 比例，检验趋势

有时候想要比较多于两个部分。但有时数据的分类可能是有序的，所以希望找到一个随着分组序号递增或递减的趋势。 这一节的例子是，一组女性是否使用剖腹产生育孩子，以及鞋子码数的数据。数据集如下：

library(ISwR)
caesar.shoe

##     <4  4 4.5  5 5.5  6+
## Yes  5  7   6  7   8  10
## No  17 28  36 41  46 140

#为在caesar.shoe这样的数据上使用prop.test函数
#需要将其转换为一个“成功”的向量（在这里和相反的很像）和一个“实验”的向量
caesar.shoe.yes <- caesar.shoe["Yes",]
caesar.shoe.total <- margin.table(caesar.shoe,2)
caesar.shoe.yes

##  <4   4 4.5   5 5.5  6+ 
##   5   7   6   7   8  10

caesar.shoe.total

##  <4   4 4.5   5 5.5  6+ 
##  22  35  42  48  54 150

#这样就很容易进行检验
prop.test(caesar.shoe.yes,caesar.shoe.total)

## Warning in prop.test(caesar.shoe.yes, caesar.shoe.total): Chi-squared
## approximation may be incorrect

## 
##  6-sample test for equality of proportions without continuity
##  correction
## 
## data:  caesar.shoe.yes out of caesar.shoe.total
## X-squared = 9.2874, df = 5, p-value = 0.09814
## alternative hypothesis: two.sided
## sample estimates:
##     prop 1     prop 2     prop 3     prop 4     prop 5     prop 6 
## 0.22727273 0.20000000 0.14285714 0.14583333 0.14814815 0.06666667

可以看出，这个检验的结果是不显著的。但是从细节上说，剖腹产一组的数字太少，让人感觉这有些巧妙的不合理，此外，注意那条警告，它是因为一些格子的期望数小于5。 可以用prop.trend.test来检测不同部分的趋势。这个检测的本质是一个用每组的分数对不同部分进行的加权线性回归，其中对领斜率进行检验，就成为了一个自由度为1的\(\chi^2\)检验。

#这个函数有三个参数，x，n，score
#前两个与prop.test中的一致
#最后一个是赋予每组的分数，默认是简单的1,2，……，k
prop.trend.test(caesar.shoe.yes,caesar.shoe.total)

## 
##  Chi-squared Test for Trend in Proportions
## 
## data:  caesar.shoe.yes out of caesar.shoe.total ,
##  using scores: 1 2 3 4 5 6
## X-squared = 8.0237, df = 1, p-value = 0.004617

所以，若假设鞋子码数是线性的，那么可以看到一个显著的区别。这样的假设并不是为了保证检验的合理所必需的。反之，它表示了这个检验所对应得备择假设。

8.4 r×c表格

为了分析两边都多于两个类的表格数据，可以用函数chisq.test或fisher.test，不过要注意后者在每一铬数字比较大而且超过两行或两列时的计算量非常大。在一个简单的例子里，已经介绍了chisq.test，不过对于更大的表而言，还有一些其他有趣的特征。 一个\(r\times c\)表格看起来像这样：

\[ \begin{array}{cccc|c} n_{11} & n_{12} & \cdots & n_{1c} & n_{1.}\\ n_{21} & n_{22} & \cdots & n_{2c} & n_{2.}\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ n_{r1} & n_{r2} & \cdots & n_{rc} & n_{r.}\\ \hline n_{.1} & n_{.2} & \cdots & n_{.c} & n_{..}\\ \end{array} \]

这样的一个表会导致几种不同的抽样方法，从而对应的表述“行和列之间没有关系”也有不同的含义了。每一行的总和可能是事先确定的，此时可能会检验每列的分布是否在每行上都一样。反之，如果每列的和是固定的，也可能是每一个个体随机地根据行或列来分组的。在后一种情况中，你应该会对检验统计独立性这个假设感兴趣，即一个个体掉入第\(ij\)个格子的概率是边际概率的乘积：\(p_{i.}p{.j}\)。不过，这些不同的情况最后都使用了同样的分析方法。 如果行和列之间没有关系，那么每一格的期望值应该是： \[ E_{ij}=\frac{n_{i.}\times n_{.j}}{n_{..}} \] 这可以理解为把每一行的总数按照每一列总数的比例（或者反过来）进行分布，或者是将整个表格的总数按照行和列的比例进行分布。 检验统计量 \[ X^2=\sum{\frac{(O-E)^2}{E}} \] 服从一个自由度为\((r-1)\times (c-1)\)近似的\(\chi^2\)分布。这里是对整个表格求和，然后下标\(ij\)被省略了。这里\(O\)表示观测值，而\(E\)表示前文所述的期望值。 下面用4.5节提到的婚姻状况与咖啡因消费情况的表格来计算这个\(\chi^2\)检验：

caff.matital <- matrix(c(652,1537,598,242,36,46,28,21,218
                         ,327,106,67),
                       nrow = 3,byrow = T)
colnames(caff.matital) <- c("0","1-150","151-300",">300")
rownames(caff.matital) <- c("Married","Prev.married","Single")
caff.matital

##                0 1-150 151-300 >300
## Married      652  1537     598  242
## Prev.married  36    46      28   21
## Single       218   327     106   67

chisq.test(caff.matital)

## 
##  Pearson's Chi-squared test
## 
## data:  caff.matital
## X-squared = 47.386, df = 6, p-value = 1.567e-08

检验结果高度显著，所以可以放心地拒绝独立性的假设。不过，一般来说，也会想知道偏差的程度。为了这个目的，可以仔细查看函数chisq.test的一些额外的返回值。 注意，函数chisq.test（就像函数lm一样）的返回值实际上比显示出来的信息更丰富：

chisq.test(caff.matital)$expected

##                      0      1-150   151-300     >300
## Married      707.65188 1491.84889 571.74523 257.7540
## Prev.married  30.60495   64.52037  24.72718  11.1475
## Single       167.74317  353.63074 135.52759  61.0985

chisq.test(caff.matital)$observed

##                0 1-150 151-300 >300
## Married      652  1537     598  242
## Prev.married  36    46      28   21
## Single       218   327     106   67

#接下来便可以对这两个表格进行彻底检查，看看差别在哪里
#检查每个格子对整体χ2的贡献通常是有用的
#这样的表格不能直接被析取，但很容易计算出来
E <- chisq.test(caff.matital)$expected
O <- chisq.test(caff.matital)$observed
(O-E)^2/E

##                       0    1-150   151-300      >300
## Married       4.3766322 1.366507 1.2056296 0.9628887
## Prev.married  0.9510407 5.316215 0.4331815 8.7079428
## Single       15.0572411 2.005471 6.4332189 0.5700246

这里有一些格子有很大的贡献，特别是很多“戒除咖啡因”的单身者，同时曾经结过婚的分布方向则被转换为朝向更多的摄入——他们摄入的更多。不过，要在这些数据里找到一个偏离独立性的简单描述仍然是不容易的。 也可以对原始（没有成为列表）数据使用chisq.test，这里用4.5节提到的juul数据。

attach(juul)
chisq.test(tanner,sex)

## 
##  Pearson's Chi-squared test
## 
## data:  tanner and sex
## X-squared = 28.867, df = 4, p-value = 8.318e-06

对这两个变量检验独立性或许没什么意义，因为Tanner Stage的定义就是与性别有关的。

练习题

1. 再来考虑练习题3.3中的情况。其中有连续10个病人进行了手术都没有并发症，而其期望的概率是20%。通过二项分布计算相关的单边检验。需要多大的样本量（仍然完全没有并发症）来使得检验结果变成统计上显著的？
2. 在美国西部的洛杉矶斑疹热事件中，747例病患死亡了210人，而东部的661例病患中死亡了122人。这个差异在统计上是显著的吗？
3. 对两种治疗胃溃疡的药物进行比较，结果如下

计算\(\chi^2\)检验和Fisher精确检验，并讨论它们的不同点。给出治愈率之差的95%置信区间。

4. 从1968年9月20日至1969年2月1日，一位老师买了254个蛋。他每天都记录在蒸的过程中有多少个蛋破裂到蛋清流出来的程度，以及有多少个蛋破裂但是蛋清没有流出。此外，他还记录了这个蛋尺寸A还是尺寸B。从1969年2月4日至1969你那4月10日，他还买了130个蛋，但是这次他用一个开孔器在蛋上开了个小孔，以防破裂。结果如下：

分析一下这样做是否有效果。

5. 对于15次试验中有3个成功的现象，在成功概率是\(x\)的时候进行双边检验，让\(x\)从0到1以0.001为间隔进行变化，并且对检验的\(p\)值作图。解释为什么对双边置信区间的定义比较困难。