2426 字

列联表分析那点事儿

列联表是一类常见表格,其行与列代表了两种分类变量,行列交织给出分别隶属两个分类的频数或者概率,其基本形式如下:

\(n_{11}\) $n_{12} $ \(n_{1+}\)
\(n_{21}\) $n_{22} $ \(n_{2+}\)
\(n_{+1}\) \(n_{+2}\) \(n_{ij}\)

简单说就是4个数据间的故事,当将其中组成部分的频数除以总数,我们可以得到概率列联表,所有表格概率综合为1。另外,固定某行(或列)求和,可以得到该行(或列)对应属性的总频数,除以总数可得某属性在样本总体中的概率。

控制实验

在控制实验中,如果列表示处理是否有效而行表示是否处理,其列联表常用来判定某一处理是否有效或某一因子是否起作用。处理间的对比用相对风险(relative risk, RR)表示,直观上理解为两个处理是否产生变化,也就是用\(\frac{\frac{n_{11}}{n_{1+}}}{\frac{n_{21}}{n_{2+}}}\)来进行估计。RR的分布为独立变量的除法不易计算,所以取对数来计算标准误(取对数后符合正态假设),这个标准误的对数估计为\(\sqrt{\frac{n_{12}}{n_{11}n_{1+}} + \frac{n_{22}}{n_{21}n_{2+}}}\),证明用到Delta Method,另一个前提为抽样比例的标准误估计为\(\sqrt{\frac{p \cdot (1-p)}{n}}\),据此可计算RR的置信区间。

另一个常用的对比描述为胜率比(odds risk, OR),有效比例较低时近似于RR,其公式为\(\frac{\frac{n_{11}}{n_{12}}}{\frac{n_{21}}{n_{22}}}\),同样,其标准误估计为\(\sqrt{\frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}}}\),可据此计算OR的置信区间。

小样本

在样本数较大时,卡方检验会得到较好的结果;但样本数较小时,由于正态近似不能达到,并不能保证错误率,可以使用Fisher exact test进行保证\(\alpha\)

关于Fisher exact test,有个很经典的故事叫女士品茶,后人以此为书名写了一本统计学史的科普书,十分精彩。精确检验的核心在于H0假设的模拟,也就是处理前后没区别,这时两独立变量也就是处理前后的概率服从超几何分布,这样可以精确计算出所有可能,进而得到精确的可能性,也就是p值,也可进行两尾检验。同样的,我们也可以用Monte Carlo模拟方法来对这个值进行估计,本质一样。以下为R中测试代码:

dat <- matrix(c(4, 1, 2, 3), 2)
fisher.test(dat, alternative = "greater")

大样本

样本数较大时,精确检验同样有效,但这时使用卡方检验在推断上功效可能更好。卡方统计量的定义为\(\sum \frac{(Observed − Expected)^2}{Expected}\),卡方统计量的本质为多个iid的平方在一定自由度下的分布,从公式上看可以理解为Expected在计数条件下既是期望也是方差,实际的证明可以参考这里。在2*2列联表里,出现对比的就两个期望,所以卡方统计量的自由度始终为1。在R中,进行卡方检验的代码为:

dat <- matrix(c(44, 77, 56, 43), 2)
chisq.test(dat)
chisq.test(dat, correct = FALSE)
chisq.test(x, simulate.p.value = TRUE)

其中,correct为连续性校正,默认开启。simulate.p.value表示用permutation test进行p值的精确计算。卡方检验的适用度很广,非常适合计数变量列联表的处理,特别是样本的独立性检验或符合度检验,但样本数小于10还是推荐精确检验来提高功效。

Simpson 悖论

如果用列联表进行统计推断,那么很难绕过Simpson悖论:当存在第三方因子变量的干扰时,分层概率与边际概率经常会得到相反的结论。其实从数学角度并不算悖论,满足下面关系的数字可以轻松构建一个Simpson悖论:

\(\frac{a}{b} < \frac{c}{d}\)

\(\frac{e}{f} < \frac{g}{h}\)

\(\frac{a + e}{b + f} > \frac{c + g}{d + h}\)

在统计推断上确实这个悖论不好解释,因为这涉及到实际的因果推断。统计上有时是会考虑分层抽样所造成的混杂变量影响,一种方法就是将分层造成的方差均一化加权,如果考虑计算的是胜率比,可以使用CMH test,其H0假设为所有分层间的胜率是一致为1的。在R中可用以下代码实现:

dat <- array(c(11, 10, 25, 27, 
16, 22, 4, 10,
14, 7, 5, 12,
2, 1, 14, 16,
6, 0, 11, 12,
1, 0, 10, 10,
1, 1, 4, 8,
4, 6, 2, 1),
c(2, 2, 8))
mantelhaen.test(dat, correct = FALSE)

同样,correct为连续性校正,也可使用exact = TRUE来进行精确计算。这样,我们可以得到某种分层策略会不会影响到胜率变化或者简单理解为胜率的异质性检验。最后,在推断上可以这样考虑:边际胜率比适合决策,条件胜率比适合分析影响。

重看控制实验

控制实验不同于观察实验,其中的一些设计方法可以对应更有针对性的列联表分析方法或者满足特殊假设下的分析要求。在研究疾病问题上,大面积采样并不容易,所以我们可以采用前瞻式实验设计或回顾式实验设计,前者固定因素数,后者固定疾病数。于此对应可以计算胜率比,其实数学上这个数没什么变化。但在计算标准误上可以根据更精确的方法来计算置信区间,这时就可以用到permutation test了,这样对于胜率比的置信区间估计相对功效高一些。

此外,当控制实验是面向同一组人时,我们会去考虑是否存在配对,这样上述分析的独立性假设可能就不满足了。这时考虑因子影响的话其实是对\(n_{12}\)\(n_{21}\)的差异进行卡方检验,这样的检验H0为配对前后无变化,称做McNemar’s test。在R中,通过下面的方法实现:

mcnemar.test(matrix(c(794, 86, 150, 570), 2),
correct = FALSE)

同样的,correct为连续性校正。从理解上,配对在某种程度上考虑了内在方差一致性问题,因此该检验与CMH test在某种程度上是相通的,同样可以进行精确检验。

非参方法

上面设计参数估计都会有假设的分布,非参方法虽不涉及参数分布的假设,但一般会考虑样本为iid。非参方法有的完全不考虑数值,有的则会有一定的加权,具体情况可根据实际需要进行。R中常用wilcox test来进行非参检验:

diff <- c(.07, .07, .00, -.04, ...)
wilcox.test(diff, exact = FALSE)

在处理列联表问题上,非参方法更不受异常值影响,但损失功效,所以可采用精确检验的思想来进行,比较适合基因组数据的处理。

小结

列联表问题的处理上,小样本用精确检验,大样本用卡方检验,混杂因素考虑CMH检验,如果数据是配对的则采用McNemar检验,感觉有异常值就wilcox检验,在计算成本不大的情况下,给出精确的p值可能更有说服力。