分类目录归档:生物与医学统计

生物统计与医学统计

医学统计学系列

非常感谢统计之都盛情邀请在卫生统计方面撰写系列文章,虽然我已经在自己的博客“卫生统计空间”写了百余篇统计文章,不过那些都是兴之所至,随手即兴而做,并无什么系统性。这次既然专门写这方面,我想尽量写得系统一些,所以这一次是第一次写,先不写具体方法什么的,而是写点学习统计最重要的东西,也就是医学统计学的重要性。

很多卫生统计学老师上课从不讲医学统计学是干什么的,有什么重要性,以至于好多学生毕业后依然不知道学了卫生统计学到底有什么用。不少老师都是第一节课开场白就是“统计学是……的科学”,反正我当年第一感觉就是不知所云,脑子一片空白,就知道跟着学习书中的公式、条件等等,完全不知学了有什么用。时至今日,依然存在这样的老师,也依然存在像我当年一样晕的学生。我希望各位在看了这个文章后能够清醒,能够明确地知道统计学到底有什么用。
继续阅读医学统计学系列

分组检测方法和 binGroup 包

本文作者:张博安,University of Nebraska统计系在读博士

$\qquad$今天给大家介绍一下分组检测(group testing)方法和我们写的关于该方法的 R 包 binGroup。分组检测(又叫 pooled testing)主要用在样本检测当中,就是把一定数量的单个样本混合在一起,然后对混合样本(称为组;group)检测是否有某种特征。举一个例子,现在要检验 1000 个血液样本是否有艾滋病毒。如果对所有单个样本挨个检测(称为单体检测;individual testing),费时费力并且花费很大。如果我们把每四个单个样本混合,我们只需要对 250 个混合样本进行检测。对检验呈阳性的组,我们可以再对其中每个单个样本进行再检测。当具有所检测的特征样本比例比较小的时候,分组检测可以大幅减少检验的次数,从而节省时间和成本。因而分组检测在传染病监测(Gaydos, 2005),药物研发(Remlinger et al., 2006),基因分型(Chi et al., 2009)等各种有关样本检测的领域中都有成功和广泛的应用。

1.均一总体(HOMOGENEOUS POPULATION)

$\qquad$分组检测方法在很长一段时间里都只用在均一(homogeneous)的总体上。在这个总体里,我们假定所有个体都是独立的,并且具有所检测特征的概率均为 $p$。$p$ 称为该特征的总体流行率(overall prevalence)。把 $Y_{ik}$ 记作在第 $k$ 组中样本 $i$ 的值,($Y_{ik}=0$ 表明该样本没有所检测的特征;$Y_{ik}=1$ 表明该样本具有所检测的特征)。把 $I_k$ 记做第 $k$ 组中样本个数(称为组大小;group size),则对于 $ i = 1,\dots, I_k$,$k = 1,\dots, K$,$Y_{ik}$ 是独立同分布的 Bernoulli(p) 随机变量。类似的,把 $Z_k$ 记做第 $k$ 组的观测值($\theta_{k}=0$ 表明该组检测呈阴性,$\theta_{k}=1$ 表明该组检测呈阳性),则对于 $k = 1,\dots, K$,$Z_k$ 是 Bernoulli($\theta_k$) 随机变量。我们的目的是估计 $p = P(Y_{ik} = 1)$,但是因为 $Y_{ik}$ 无法观测到,我们需要找出 $p$ 和 $Z_k$ 间的关系。 假定所有检测是完美的,我们有 $\theta_k=0 \Leftrightarrow \sum^{I_k}_{i=1}Y_{ik}=0$ 和 $\theta_k=1 \Leftrightarrow \sum^{I_k}_{i=1}Y_{ik}>0$。则不难推出 $\theta_k=P(Z_k=1)=1-P(Z_k=0)=1-(1-p)^{I_k}$。又因为 $Z_k$ 是独立的 Bernoulli 随机变量,我们可以写出如下的似然函数:

$L(p|z_1,\dots,z_k)=\sum^{K}_{k=1}[1-(1-P)^{I_k}]^{Z_k}(1-p)^{I_k(1-Z_k)}$

 再由MLE方法求得 $\widehat{p}$,并通过Fisher信息矩阵算出 $\widehat{p}$ 方差。我们写了两个 R 函数计算 $p$ 的置信区间。bgtCI() 用来计算当所有组有相同大小时 $p$ 的置信区间。下面我们用一个例子说明如何使用 bgtCI() 。刘沛等(1997)这篇论文研究了在我国徐州地区丙型肝炎的流行率。实验者把 1875 个献血者的血液样本每五个混合(组大小为 5),再用 ELISA 检验试剂对 375 个组进行检测。结果有 37 个组检验呈阳性。则徐州地区丙型肝炎的总体流行率 $p$ 的 95% 置信区间可以如下计算:

> bgtCI(n = 375, y = 37, s = 5,
+  conf.level = 0.95,
+  alternative = "two.sided",
+  method = "AC") 

95 percent AC confidence interval:
 [ 0.01487, 0.02821 ]
Point estimate: 0.02056

这里我们使用了 Agresti-Coull 置信区间,因为它通常有很好的覆盖概率。还有一些其它的置信区间可供选择,譬如 Clopper-Pearson,Wilson 和 Wald 区间(覆盖概率最差)等。如果大家学过属性数据分析,相信对这些二项分布参数的置信区间 不陌生。在组大小相等的情况下,我们有 [latex\theta=P(Z_k=1)=1-(1-p)^l]$。我们先算出关于 $\theta$ 的置信区间,再通过以上变换求得 $p$ 的置信区间。关于这些区间在分组检测情况下的比较,请参看 Tebbs and Bilder (2004)。当各组组大小不等的时候,置信区间的计算就变得复杂许多。我们提供了 bgtvs() 函数计算这种情形下的精确置信区间,这里就不详细介绍了。

$\qquad$分组检测实验中很重要的一个步骤就是选择合适的组大小。如果组太小,很少的混合样本会检测呈阳性,我们就浪费了检测试剂;如果组太大,那么大部分的混合样本会检测呈阳性,对 $p$ 的估计就会很差。一个凭经验的方法是选一个组大小使得差不多一半的组检测呈阳性。更精确的方法是,当我们对 $p$ 有一个初步估计,选一个组大小使得 MSE(均方误差)最小化(Swallow, 1985)。在这个例子中,假定我们初步估计 $p$ 为 0.025,并设定组大小的上界为 100,则 estDesign() 函数可以计算出最佳的组大小是 61。

> estDesign(n = 375, smax = 100, p.tr = 0.025)
group size s with minimal mse(p) = 61
$varp [1] 2.554086e-06 

$mse
[1] 2.560173e-06 

$bias
[1] 7.80204e-05 

$exp
[1] 0.02507802

当然,在实践中我们通常无法混合如此多个单个样本(阳性样本在被过度稀释后会检测呈阴性,称为稀释效应)。但是这个函数可以给我们选取组大小提供参考。对均一总体,我们的包中还包括以下函数:

  • bgtTest(): 计算对 p 假设检验的 p 值
  • bgtPower(): 计算对 p 假设检验的功效
  • nDeisgn(), sDesign(), plot.bgtDesign(): 当 n 或 s 变化时,计算假设检验的功效并作图。

这些函数可以帮助使用者设计自己的分组检测实验,这里我们就不一一介绍。

2.非均一总体(HETEROGENEOUS POPULATION)

$\qquad$在现实中均一总体的假设总显得不切实际。拿之前那个的例子来说,酗酒者,饮食卫生条件较差或卫生习惯不良者会比其他人得肝炎的几率较大,所以认为所有人有同样感染肝炎的风险是不合常理的。如果能搜集到每一个献血者个人信息的数据,我们希望建立一个关于肝炎的流行率的回归模型,使我们能通过各种相关的因素预测个人感染肝炎的概率。当然,这不是一个简单的 logistic 模型,因为在分组检测中我们无法收集到个人是否患病的数据。我们需要新的估计参数的方法。

$\qquad$我们仍然把 $Y_{ik}$ 记做在第 $k$ 组中样本 $i$ 的值,则 $p_{ik}=P(Y_{ik}=1)$ 是该样本具有所检测的特征的概率。把可能影响此概率的 $p-1$ 个变量记做 $x_{ik1}, x_{ik2},\dots, x_{ik,p-1}$,我们的回归模型是

$f(p_{ik})=\beta_{0}+\beta_{1}x_{1ik}+\dots+\beta_{p-1}x_{p-1,ik}$

 其中 $f$ 是链接函数。对于 $i = 1, \dots, I_k$, $k = 1,\dots, K$,$Y_{ik}$ 是独立的 Bernoulli($p_{ik}$) 随机变量。目前主要有两种方法估计 。第一种方法是由 Vansteelandt et al. (2000)提出的。假定所有检测都是完美的,则不难推出

$\theta_{k}=1-\sum^{I_k}_{i=1}(1-p_{ik})$

又因为 $Z_k$ 是独立的 Bernoulli 随机变量,我们可以写出如下的似然函数:

$L(\beta_0,\dots,\beta_{p-1}|Z_1,\dots,Z_k)$

$=\prod^{K}_{k=1}(1-\prod^{I_k}_{i=1}(1-p_{ik}))^{z_K}(\prod^{I_k}_{i=1}(1-p_{ik}))^{1-z_K}$

$=\prod^{K}_{k=1}[1-\prod^{I_k}_{i=1}(1-f^{-1}(\beta_0+\beta_1x_{1ik}+\dots+\beta_{p-1}x_{p-1,ik}))]^{Z_k}$

$\times [\prod^{I_k}_{i=1}(1-f^{-1}(\beta_0+\beta_1x_{1ik}+\dots+\beta_{p-1}x_{p-1,ik}))]^{1-Z_k}$

再由 MLE 估计 $\beta$。在包中我们是用 optim 函数实现 MLE 估计的。这个方法浅显易懂,但是缺点是不够灵活。当我们对阳性组中单个样本再检测的时候,这些再检验的结果无法被包括在参数估计中。另外一种更灵活的方法是由 Xie (2001)提出的。假设我们有每个单个样本的值,则似然函数具有简单的二项分布的形式:

$L(\beta_0,\dots,\beta_{p-1}|y_{11},\dots,y_{I_kK})=\prod^{K}_{k=1}\prod^{I_k}_{i=1}p^{y_{ik}}_{ik}(1-p_{ik})^{1-y_{ik}}$

 

因为单个样本无法观测到,我们将上式中的每一个 $y_{ik}$ 替换为 $E(Y_{ik}|Z_1=z_1,\dots,Z_k= z_k)$,再利用 EM 算法估计 $\beta$。EM 算法是在参数估计中非常重要的方法,在处理缺失数据和隐性变量模型等领域应用广泛。有关 EM 算法的简介,可以参看 Casella and Berger (2001, p. 326 – 329)。对我们的问题,EM 算法由下面给出:

1. E-step:计算$\widehat{w}_{ik}=E(Y_{ik}|Z_1=z_1,\dots,Z_k=z_k)$
2. M-Step:找出β使得以下似然函数的条件期望最大化(其中 $p_{ik}$ 是 $\beta$ 的函数)$E[log(L(\beta|y_{11},\dots,y_{I_kK}))|Z_1=z_1,\dots,Z_k=z_k]=\sum^{K}_{k=1}\sum^{I_k}_{i=1}\widehat{w}_{ik}log(p_{ik})+(1-\widehat{w}_{ik})log(1-p_{ik})$
3. 重复步骤1和2直到 $|(\widehat{\beta}^{(r)}_d-\widehat{\beta}^{(r-1)}_d)/\widehat{\beta}^{(r-1)}_d|<\epsilon$。其中$\widehat{\beta}^{(r)}$是$\beta$的第r次估计,$\epsilon > 0$

$\beta$ 的方差可以由 Louis’formula 得到。该方法的优点是非常灵活,可容纳对单个样本再检验,以及应用在如阵列分组(matrix or array-based pooling, Phatarfod and Sudbury, 1994; Kim et al., 2007)等更复杂的分组检测的场合。在我们的包中,gtreg 函数用来拟合分组检测回归模型。在使用时,我们可以指定检验的敏感性(sensitivity)和特异性(specificity)(在以上的介绍中,为简单起见,我们没有讨论检测有误差的情形),分组编号,参数估计的方法(Vansteelandt 或 Xie)等。如果我们要包含对单个样本再检验结果,则只能使用 Xie 方法;除此以外,两种方法应给出非常相似的结果。

$\qquad$我们用 Vansteelandt et al. (2000)中的一个例子来说明如何使用 gtreg 函数。数据来自一项对肯尼亚偏远地区怀孕妇女 HIV 感染率监测的研究,研究者收集了所有怀孕妇女的个人信息,这里我们选取年龄和教育水平两项变量预测每位妇女感染 HIV 的概率。数据的最后几行如下:

> data(hivsurv)
> tail(hivsurv[,c(3,5,6:8)], n = 7)
    AGE EDUC. HIV gnum groupres
422  29     3   1   85        1
423  17     2   0   85        1
424  18     2   0   85        1
425  18     2   0   85        1
426  22     3   0   86        0
427  30     2   0   86        0
428  34     3   0   86        0

可以看到每一组中单个样本都有同样的分组编号(gnum)和)值(groupres).例如第 422 个样本是 HIV 阳性,所以在第 85 组中的所有样本都有 “1” 的组值。这组数据中也包括单体检测的数据(HIV),因为实验者希望证实分组检测在节约大量的成本的同时可以得到和单体检测几乎一样好的估计。现在我们就用 gtreg 函数拟合数据,并用包中的 summary 函数给出总结性的输出。

> fit1  summary(fit1)

Call:
gtreg(formula = groupres ~ AGE + EDUC., data = hivsurv, groupn = gnum,
    sens = 0.99, spec = 0.95, linkf = "logit", method = "Vansteelandt") 

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-1.1813  -0.9385  -0.8221   1.3297   1.6694   

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.99039    1.59911  -1.870   0.0615 .
AGE         -0.05163    0.06748  -0.765   0.4443
EDUC.        0.73621    0.43885   1.678   0.0934 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    Null deviance: 191.4  on 85  degrees of freedom
Residual deviance: 109.4  on 83  degrees of freedom
AIC: 115.4 

Number of iterations in optim(): 138

拟合的结果存放在 fit1 中,我们定义它为”gt” class。得到的模型可以写成

$logit(\widehat{p}_{ik})=-2.99-0.0516Age_{ik}+0.7362Educ_{ik}$

其中 $\widehat{p}_{ik}$ 是第 $k$ 组中第 $i$ 个样本 HIV 阳性的概率的估计。

$\qquad$在 binGroup 包中还包括一些其他函数。如 sim.gt 产生分组检测数据,predict.gt 根据参数的估计值预测 $p_{ik}$ 等。gtreg.mp 函数适用于在基因分型领域有非常重要应用的阵列分组。在这类实验中,所有单个样本被放在一个或多个阵列的方格中(例如 $n^2$ 个样本就分配到 $n \times n$ 个方格里), 然后将每一行的样本作为一组,每一列的样本作为一组检测。如果某一行和某一列检测呈阳性,则在它们交叉的位置的单个样本很可能具有所检测的特征。这种分组方式特别之处在于每个单个样本同时出现在所在行和所在列的组中,所以 Vansteelandt 方法无法应对这种分组方式。而对 Xie 方法,我们只需把所有行的观测值,列的观测值作为已知信息包含在 E-step 中的条件期望,就可通过 EM 算法求得参数估计。

$\qquad$值得一提的是,在很多分组检测的的场合(如对单个样本再检验,阵列分组等),E-step 中的条件期望很难求出,或者没有解析表达,这种情况下我们可以使用 Gibbs sampling 方法估计条件期望。Gibbs sampling(可以参看 Carlin and Louis, 2008, Section 3.4.1)方法是 Metropolis-Hastings 算法的一个特例。简单来说,当多个变量的联合概率分布不明确,而各个变量的条件分布已知的时候,Gibbs sampling 根据其他变量的当前值,依次对分布的每个变量生成一个样本,最后建立一个马尔可夫链,其平衡分布就是这多个变量的联合分布。为叙述方便,我们把所有单个样本的值重新标记为 $Y_1,\dots, Y_N$,把所有可观测到的变量(所有组值,再检验结果等)记作 $T$。因为 $Y_1,\dots, Y_N|T=t$ 的分布很难求出,而对每一个 $i$,条件分布 $Y_i|Y_s=y_s,s \ne i,T=t$ 则容易求得,所以我们可以生成 Gibbs 样本 $y^{*}_i \sim f(Y_i|Y_s=y_s,s \ne i,T=t)$。在 $N$ 个单个样本中循环 $K$ 次($K$ 通常很大)后,我们得到一个 Monte Carlo Markov Chain 并且 $(y^{*}_1,\dots,y^{*}_N)$ 的联合分布收敛到 $Y_1,\dots, Y_N|T=t$。现在我们可以用 $\sum y^{*}_i/K$ 估计 $\widehat{w}_i=E(Y_i|T=t)$ ,其中 $\sum$ 是对所有 $K$ 个 Gibbs 样本求和。这时的 EM 算法称作 Monte Carlo Expected Maximization(MCEM)算法。当然,因为在每个 E-step 中我们要产生大量的 Gibbs 样本,这种方法通常较慢。

$\qquad$以上向大家简单介绍了分组检测的一些方法及其在 R 中的实现。我们这里只讨论了分组检测中参数估计的问题。还有很多方法(参看 Bilder, Tebbs, and Chen, 2010)专注于如何通过再检验最快找出所有 $Y_{ik}=1$ 的样本,这些方法反过来也依赖于我们估计的 $\widehat{p}_{ik}$(在阳性组中优先再检验 $\widehat{p}_{ik}$ 大的样本)。如果你的工作中能用到这些方法,欢迎你使用我们的 binGroup 包以及向我咨询;如果和你的工作没有直接的联系,我们的模型和 R 程序也涉及了统计学中一些热门的方法,希望能对大家有所帮助和启发。

下载本文PDF文档: 分组检测方法和binGroup包

假设检验初步

准备尝试一下,用大白话叙述一遍统计推断中最基础的东西(假设检验、P值、……),算是把这段时间的阅读和思考做个梳理(东西不难,思考侧重在如何表述和展示)。这次打算用一种“迂回的”表达方式,比如,本文从我们的日常逻辑推理开始说起。

0.普通逻辑

复习一下普通逻辑的基本思路。假设以下陈述为真:

你打了某种疫苗P,就不会得某种流行病Q。

我们把这个先决条件表述如下:

如果P 则非Q

其中,

P表示打了疫苗P,

Q表示得流行病Q

或者,更形式化一点:

if P then NOT Q

然后,如果观察到你得了流行病Q,那么就可以推出你没有打疫苗P——这个推断只不过是上述前提条件的逆反命题而已。我们把以上推理过程表述如下:

if P then NOT Q (先决条件)

Q (前提)

———————–

then NOT P (结论)

还有,如果你没有得流行病Q,就能推断出你打了疫苗P吗?显然不能。打疫苗P是不得流行病Q的充分条件,但非必要条件:你没有得流行病Q,可能是因为打了疫苗P,也可以是因为其他任何原因。即,if P then NOT Q,不能够推出if NOT Q then P。

到此为止没有任何令人惊奇的地方。下面将表明,假设检验背后的统计推断规则也只不过是我们以上日常逻辑推理的一个衍生而已。这只需要思维的一次小小的“跳跃”。

1.假设检验

在统计推断中,我们不说“你打了疫苗P,就不会得流行病Q”,而是说,比如,“你打了疫苗P,就有95%的把握不会得流行病Q”,即if P then probably NOT Q。把上面的逻辑推理规则改写成统计推断规则:

if P then probably NOT Q    (先决条件)

Q                                                     (前提)

———————–

then probably NOT P         (结论)

回到以前“万能”的硬币实验,我们做实验来考察一枚硬币是不是均匀的。改写成现在我们熟悉的形式:

P:硬币是均匀的。

Q:在100次投掷中,得到90次正面,10次反面。

我们说,如果是一个均匀的硬币,就不太可能发生这样的情形:投100次,出现90次正面,10次反面(if P then probably NOT Q)。现在如果在100次投掷实验中,观察到出现90次正面,10次反面(Q),那就可以有把握地说,这个硬币不是均匀的(NOT P)。这个推理可以写成与上面一致的统计推断的形式,其中,P是原假设H0,NOT P是备择假设Ha:

H0:硬币是均匀的  (P

Ha:硬币是有偏的 (NOT P

如果原假设为真,即硬币是均匀的,就不太可能发生这样极端的事情,比如:在100次投掷实验中,观察到出现90次正面,10次反面(Q)。如果真的观察到这样极端的事情,你就有把握认为硬币不是均匀的,即拒绝原假设(P),接受备择假设(NOT P)。

另外,如果在100次投掷实验中,观察到60个正面,40个反面(NOT Q)。这时你就不好下结论了,因为一个均匀的硬币可能投出这样的结果,一个有偏的硬币也可能投出这样的结果。最后,你只能说,如果实验结果是这样的,那就没有把握拒绝原假设。这枚硬币是否有偏,需要更多的证据来证明(这通常意味着更多的实验,比如,再投1000次)。

总结一下。在搜集数据之前,我们把想证明的结论写成备择假设,把想拒绝的结论写成原假设。之所以写成这个形式,因为从上面不厌其烦的讨论中得知,这是方便逻辑/统计推断的形式:当我们难以拒绝原假设时,只能得到结论,原假设也许是真的,现在还不能拒绝它;而当我们能够拒绝原假设时,结论是:它就很有把握是不真的。注意,在看到数据之前,我们不知道自己想证明的结论是否能够被证据所支持。

在确定假设检验的形式的同时,我们对之前一直随意说的“把握”、“可能”也做一个限定,即指定一个显著性水平α(significance level),也叫犯第一类错误的概率(type I error,在上面的硬币实验中,就是否定一个均匀硬币的错误,也叫“弃真”错误)。

根据某些保守或稳健的原则(比如,我们认为,把一个无辜的人判决为有罪,比放掉一个有罪的人,后果更为严重),我们要尽量把犯“弃真”错误的概率控制在一个很小的水平里。通常α=0.05,这时候就是说,如果拒绝了原假设,你就有95%的把握说原假设是不真的。这里,95%(=1-α)就是置信水平(confidence level)。

又,放掉一个有罪的人,即把一个有罪的人判为无罪,这犯的是第二类错误β(type II error,在硬币实验中,就是把一个有偏的硬币当成均匀硬币的错误,也叫“取伪”错误)。关于第一类和第二类错误之间的权衡取舍(trade off),详见《决策与风险》。在我们的假设检验里,我们认为犯一类错误的后果比犯第二类错误的后果更为严重。

需要注意的是,在这里,我强调的是先提出需要检验的假设,然后再搜集收据。这是统计推断的原则之一。如果看到了数据之后再提出假设,你几乎可以得到所有你想要的结果,这是不好的机会主义的倾向。强调这些,是因为在学校里,我们大多是看了别人搜集好的数据之后再做统计练习。

事先确定好你想拒绝/证明的假设,在看到数据之前,你不知道结果如何。

2.P值(P Value)

上面提到“极端”事件,比如,在100次硬币投掷实验中,观察到出现90次正面,10次反面(Q)。怎么样的事件才是“极端的”?简单地说,一个事件很极端,那么少比它本身“更极端”的事件就非常少(比如,只有“91次正面,9次反面”、“91次反面,9次正面”等情况才比它更极端)。

但这个Q只是从一次实验中得出的。我们可以重复做这个实验,比如100次,每次都投掷100次,记录下的正面数X,它构成一个二项分布,X~B(n,p),其中,n=100,p=0.5。根据某个中心极限定理,正态分布是二项分布的极限分布,上面的二项分布可以由均值为np=50,方差为np(1-p)=25的正态分布来近似。我们在这个近似的正态分布的两端来考察所谓“更极端”的事件,那就是正面数大于90或者小于10。

重复一遍,“P值就是当原假设为真时,所得到的样本观察结果更极端的结果出现的概率”。如果P值很小,就表明,在原假设为真的情况下出现的那个分布里面,只有很小的部分,比出现的这个事件(比如,Q)更为极端。没多少事件比Q更极端,那就很有把握说原假设不对了。

在上述近似的正态分布中,P值就等于X>90 或 X<10的概率值(记做,P{X>90 or X<10})。根据对称性,这个概率值等于2*P{X<10}=1.2442E-15。

上面我们的确求出了一个非常小的P值,但如何不含糊地确定它就是很“极端”呢? 事先确定的显著性水平α,本身就是一个判定法则。只要P值小于显著性水平α,我们就认为,在认为原假设为真的情况下出现的事件Q,是如此地极端,以至于我们不再相信原假设本身。一句话,我们的判定法则是:

P值小于显著性水平α,拒绝原假设。

3.一个手算示例

用一个双侧的单样本T检验做例子。假设我们想知道,螃蟹的平均温度,跟空气的温度(24.3)有没有统计差别(α=0.05)。事先确定的假设检验的形式表达如下:

零假设(H0):   μ=24.3°C

备择假设(Ha):  μ≠24.3°C

以下是25只螃蟹在温度为24.3°C下的体温(单位:°C):

25.8    24.6    26.1    22.9    25.1
27.3    24        24.5    23.9    26.2
24.3    24.6    23.3    25.5    28.1
24.8    23.5    26.3    25.4    25.5
23.9    27        24.8    22.9    25.4

一些基本的算术结果:

样本均值:$\bar{X}=25.3$

样本量:n=25

样本方差:$s^2$=1.8

样本均值的标准误差:$s(\bar{X})=\sqrt{s^2/n}=0.27$

这里T检验的思路如下:

  1. 我们先假设H0为真,即认为螃蟹的平均温度跟空气温度没有差异(P),  μ=24.3°C。有一个极端事件Q,如果原假设H0成立,Q就不成立(if H0 then probably NOT Q);但如果在原假设为真的情况下,出现了这么一个Q,那我们就有把握拒绝原假设。
  2. 样本均值:$\bar{X}$是总体均值μ的最好的估计,在本例中,$\bar{X}=25.03$。
  3. 这个样本均值只是一个估计值。它只是从总体的一个随机样本中得到的(样本是上述25只螃蟹)。我们不知道这次实验结果是不是“极端”事件。而判断一个事件是不是极端事件,根据第二节的讨论,我们可以重复做上述实验,比如100次,每次都抓25只螃蟹,都在空气温度24.3的状态下测量其体温,然后也各自求出一个样本均值来。
  4. 容易得出,这种实验出来样本均值,辅以适当的数学形式,就服从一个自由度为24(=25-1)的t分布,即$(\bar{X}-\mu)/s(\bar{X})\sim t(24)$。
  5. 样本均值$\bar{X}=25.03$,在这个自由度为24的t分布下,有一个对应的t值,t=25.03-24.3/0.27=2.704。现在我们可以在整个分布里考察这个t值。在这个自由度为24的t分布里,我们看 t=2.704是不是一个“极端”事件Q。根据对称性,比Q更极端的是那些大于2.704或者小于-2.704的点。

t

从上图可以看到,在这个t分布里,比t=2.704更“极端”的点占整个分布的0.0124。这个0.0124就是我们要求的P值。这个P值小于我们事先选定的显著性水平α=0.05,因此我们可以拒绝原假设,认为这批螃蟹的平均体温不等于空气温度。

这个双侧P值可以手算如下:

在SAS里,P=2*(1-probt(t,df))=2*(1-probt(2.704,24))=0.012392

在R里,     P=2*(1-pt(t,df))=2*(1-pt(2.704,24))=0.012392

———-

以上是用P值作为判定条件。一个等价的做法是用临界值来判断。我们事先给定的显著性水平α=0.05,在这个自由度为24的t分布里,就对应着一个临界t值2.064。下图的阴影部分,也称作拒绝区域。上面求出的跟样本均值$\bar{X}=25.03$对应的t值=2.704,处在这个拒绝区域内(2.704>2.064),于是我们一样拒绝原假设。

t2

又,上述临界值可以手算(或查表)如下:

在SAS里,tCritic=tinv(1-alpha/tail,df)=2.06390

其中,alpha=0.05,tail=2表示双侧检验,df=24.

在R里,tCritic=qt(1-alpha/tail,df)=2.063899

4.注

本文是对近期阅读做的一个笔记。作为一个非统计科班出身的程序员,我一直在思考,如何来理解统计概念,以及如何把自己的理解向同行传达。关于用日常逻辑推理来理解假设检验的思路,来自

Common Statistical Methods for Clinical Research with SAS Examples(2nd edition, SAS Inc., 2002, by Glenn A. Walker)

关于决策与风险的讨论,参考了

维恩堡《数理统计初级教程》(常学将等译,太原:山西人民出版社,1986,Statistics: An Intuitive Approach By George H. Weinberg and John Abraham Schumaker)

第三节示例的数据,来自

Biostatistical Analysis (5th Edition) by Jerrold H. Zar, Prentice Hall, 2009

第三节的t分布图,来自一个在线的t分布生成器(很好用):

http://onlinestatbook.com/analysis_lab/t_dist.html

附录: 用SAS来计算

上面的文字尽量做到“平台无关”。这里附出SAS例子,是想把以上的手算结果跟机器结果做个对照,让读者更有信心一些。 欢迎读者贴出自己趁手的工具得出的结果。

/*data*/
data body;
input temp @@;
h0=24.3;
diff=temp-h0;
datalines;
25.8    24.6    26.1    22.9    25.1
27.3    24      24.5    23.9    26.2
24.3    24.6    23.3    25.5    28.1
24.8    23.5    26.3    25.4    25.5
23.9    27      24.8    22.9    25.4
;

/*method 1: use proc means*/
proc means data=body T PRT;
var diff ;
run;

结果是:

t Value    Pr > |t|
——————-
2.71      0.0121
——————-

上面的t Value 就是计算出来的t值,Pr > |t| 就是P值(这里的|t|就是上面计算出来的t值2.704,Pr > |t|求的是比t值更极端的概率,即P值)。proc means没有提供临界t值(即通常说的查表得出的t值),下同。

/*method 2 (prefered): use proc ttest*/
proc ttest data=body h0=24.3 alpha=0.05;
var temp;
run;

proc ttest的结果更为丰富:

N      Mean     Std Dev  Std Err    Minimum   Maximum

25     25.0280      1.3418      0.2684 22.9000        28.1000

Mean     95% CL     Mean       Std Dev     95% CL   Std Dev

25.0280 24.4741  25.5819            1.3418             1.0477   1.8667

DF    t Value    Pr > |t|

24       2.71           0.0121

北京数据管理与生物统计论坛(BBF)第三次聚会见闻录

9月4号下午,周六,去北大医学部参加了北京数据管理与生物统计论坛(Beijing Biometrics Forum, BBF)的第三次聚会,这次活动由SAS China北京大学临床研究所赞助。这里写些会议见闻和一些零散的感想,不算是会议的正式“纪要”。东西贴这,大致想给“统计之都”(COS)的朋友交流下北京SAS技术社区的氛围、工作市场情况以及一些相关技术评论等等。

1. 话题

西安杨森的薛富波博士做的开幕致辞,他也是BBF的召集人。SAS China的市场部总监罗威先生代表SAS公司做了欢迎致辞。下面是六位主题发言人的演讲简介以及我的一些杂感。

继续阅读北京数据管理与生物统计论坛(BBF)第三次聚会见闻录

泊松低方差计数数据建模问题

本文作者为中国人民大学统计学院饶燕芳同学,由COS编辑部审核发表,略有修改。点击此处下载/阅读本文PDF版本

一、问题的引出

在数据分析和数据建模的过程中,我们通常需要假定数据变量服从某种分布,以便于建立与分布参数有关的模型或方程,之后利用观测值对参数进行估计,从而达到研究和分析的目的。由于变量是随机的,我们无法确定变量在某个情况下的具体取值,因此通过假定它服从某个分布,然而感兴趣的只是它们的平均水平,即各变量之间的关系都建立在均值的基础上,方差则用于计算估计的精度和假设检验。而大多数情况下,一旦分布的假定确定,随之确定的也就是数据必须符合该分布的均值和方差特征。对于许多单参数分布,均值和方差均有一一对应的关系,如果均值确定,方差由于是均值的函数也就自然地确定下来,例如伯努利分布具有单参数$p$,均值$\mu=p$,方差$v=p(1-p)$,即有$v=v(\mu)=\mu(1-\mu)$。在这种单参数的情况下,如果观测数据的均值符合假定(即认为$p\approx\bar{Y}$),数据的方差和均值就必须满足一定条件(即例如假定$Y$服从两点分布,认为$p\approx\bar{Y}$,则方差应该有$Var(Y)=p(1-p)\approx\bar{Y}(1-\bar{Y})$),此时若观测到方差系统地大于分布假设下(此时常被观测均值决定)的方差,就出现了所谓的“超散布性”(overdispersion),类似地,若出现方差偏小的情况,也就相应出现了“超聚集性”(underdispersion)。

具体到本文需要讨论的泊松分布:现实中常常出现方差不满足假定的情况。由于参数为 的泊松分布具有均值和方差相等的特点,如果假定服从泊松分布的数据的样本方差大于模型估计的方差——即样本均值,就出现了“超散布性”,本文称之为泊松分布高方差(extra-Poisson variation),而当样本方差低于样本均值时,称此时的“超聚集性”为泊松分布低方差,后文出现的泊松低方差都符合该定义。

正如之前所说,通常建立模型如线性回归都基于均值,因此方差违反假定分布并不影响参数估计效率,但在区间估计和假设检验时就会出现问题。当“超聚集性”出现时,真实的方差会被低估,这将会错误的表现出数据中原本不显著的差异,相反地,“超聚集性”出现时,真实的方差会被高估,这样可能无法检验出组间分布的真实差异,参数的置信区间也会给得过大。因此对于方差超扩散或超聚集的数据,方差问题的处理显得尤为重要,针对此的模型建立是该类问题分析的关键。

泊松分布的超散布性数据在现实中较为常见,简单的序列正相关和非齐次性都可能引起超散布性的出现。泊松低方差的情况则较为少见,但在医学和社会领域中却经常出现。本文的目标就在于探讨针对泊松低方差数据的分布模型。

二 两种泊松低方差问题的处理方法

泊松分布为模拟计数数据提供了良好的模型,但均值和方差相等的要求在现实中却显得太为苛刻。因此处理泊松低方差的方法探究就集中在合适的修正分布的寻找上。能够描述计数数据且具有泊松低方差特点(即均值大于方差)的分布选择并不太多,其中包括两种典型的泊松低方差模型:加权泊松分布模型和 CBR 分布模型。

1. 加权泊松分布

Rao CR (1965)提出,若随机变量Y服从加权泊松分布,其密度函数为

$
P(Y=k)=\frac{e^{-\lambda}\lambda^k\omega_k }{Wk!}, k=0,1, \ldots ; \lambda>0
$

$
W=\sum^{\infty}_{k=1}\frac{e^{-\lambda}\lambda^k\omega_k }{k!}
$                                                     (a)

它是保证求和为1的标准化因子。

本文讨论一种较为简单的权重(Martin S R and P Besbeas,2004),即

$\omega_k=\left( \begin{array}{ll} e^{\beta_1(\lambda-k)} & ,k\leq\lambda\\e^{\beta_2(\lambda-k)} & ,k>\lambda\\\end{array} \right)$

对于$\beta_1,\beta_2>0$,它的分布类似将概率密度图线向均值“挤压”(见图1),分布更加集中,相对于标准的泊松分布就有更小的方差,称该分布为三参数指数加权泊松分布,记为 EWP3 。特殊地,当$\beta_1=\beta_2=\beta$时,称为两参数指数加权泊松分布,记为 EWP2 分布,当$\beta=0$时退化为标准泊松分布。对于EWP2 和 EWP3 ,它们拥有更高的峰值,标准化因子 W 可以由式 (a) 导出。尽管矩的表达没有显式,但可以确定分布的方差随着$\beta_1,\beta_2$或$\beta$的增大而降低。

图1 lambda=10,beta1=0.1,beta2=0.2的EWP分布

对于加权泊松分布的探究使得权重还有一些别的形式,如 Castillo and Pe´rez-Casany (1998) 提出的幂加权泊松分布(PLWP),Cameron 和 Johansson (1997) 使用的泊松多元加权分布 PPp 等。

2. 纯生过程模型(CBR)

不得不提的是,在处理泊松低方差数据的问题中还有一类较为有效的方法。由Faddy(1997)在随机过程的基础上提出这种变出生概率(CBR)分布。这个分布是建立在广义泊松分布的基础上:Faddy 认为,任何关于 {0,1,2,…}的离散分布都有广义泊松特性,即纯生过程。考虑一个 Markov 计数过程,$X(t)$为$(0,t)$内的事件发生数,在$(t,t+\sigma t)$内有转移概率:

$
P\left(X(t+\sigma)=n+1|X(t)=n\right)=\lambda_n\sigma+o(\sigma)
$

其中$\lambda_n$为事件数为n时的事件发生率,我们感兴趣的只是某一时刻$X(t)$的分布,这里t可以不失一般性地取1,在此模型中,时刻1时的事件数$X$的分布具有如下形式:

$(p_1(1),p_2(1),\ldots,p_n(1))=(1,0,0,\cdots,0)exp(Q)$

其中

$p_i(1)=P\left(X(1)=i\right)$

$exp(Q)=e^Q=E+\frac{Q}{1}+\frac{Q^2}{2!}+\frac{Q^3}{3!}+\cdots$

$
Q=\left( \begin{array}{ccc} -\lambda_1 & \lambda_1 & 0 & \ldots & 0\\0 &-\lambda_2 & \lambda_2 & \ldots & 0\\0 & 0 & -\lambda_3 & \ldots & \vdots\\ \vdots & \vdots & \vdots & \ddots & \lambda_{n-1}\\ 0 & 0 & 0 & \ldots & -\lambda_n\\\end{array} \right)
$

而上式的成立可由 Kolmogorov 微分方程简单推导而来:

$p'(t)=Qp(t); p(t)=ce^{Qt}; p(1)=ce^Q;$

$p(0)=c; p(1)=p_0e^Q; p(0)=(1,0,0,\ldots,0)$

这里认为初始时刻的事件数是从1开始的。因此,CBR 分布是由一系列不同的事件发生率参数$\left(\lambda_1,\lambda_2,\ldots,\lambda_k,\ldots\right)$决定的。通常可以认为$\lambda_k$是k的函数。Faddy 在1997年已经证明,对于递增的$\left(\lambda_1<\lambda_2<\cdots<\lambda_n<\cdots\right)$,$X(t)$将表现出泊松高方差特征,而当$\lambda_1>\lambda_2>\cdots>\lambda_n>\cdots$递减时,也就表现出泊松低方差特征。

三、参数估计

上述两种分布的参数估计都可通过极大似然法求出。记$x_i$为第i个样本的事件发生数,观测数据中中事件数k的频数$f_k$(k=1,2,3,…),则 EWP2 和 EWP3 分布的负对数似然方程为(已去除与参数无关的项$lnk!$):

$-LnL(\lambda,\beta_1,\beta_2)=n[\lambda-\bar{x}ln\lambda+lnW]+\beta_1\sum^{[\lambda]}_{k=x_{min}}(\lambda-k)f_k+
\beta_2\sum^{x_{max}}_{k=[\lambda]+1}(k-\lambda)f_k$ (b)

通过求使(b)式达到最小值的$\hat{\lambda},\hat{\beta}_1,\hat{\beta}_2$得到估计参数。

对于纯生过程模型,概率分布向量$(p_1(1),p_2(1),\ldots,p_N(1))$就是矩阵$exp(Q)$的第一行,若$N=x_{max}$,其负对数似然函数为:

$-LnL(\lambda_1,\lambda_2,\ldots,\lambda_N)=\sum^N_{k=1}f_klnp_k(1)$

通过最小化上式即可得到$(\lambda_1,\lambda_2,\ldots,\lambda_N)$的极大似然估计。而参数估计的方差可以通过数值计算时产生的 Hessian 矩阵得到。

四、数据实例

我们引用 Faddy(2001)的小鼠胚胎数据,作者已对该数据用 CBR 方法做了较好的分析。在产生该数据的实验中,对已经怀孕的小鼠用药(除草剂2,4,5-T),同时记录小鼠子宫上的胚胎着床数。数据给出了7种剂量水平下胚胎着床数的频率分布。在交配后的6-14天内,往怀孕的雌鼠体内注射某种剂量水平的2,4,5-T。在分娩之前,将雌鼠杀死,确定其体内的活胎数目。0剂量组的频数分布便具有泊松低方差特征。0剂量组数据的均值为11.55,方差为9.92,方差均值比为0.859。

图2

使用泊松分布对参数进行估计,参数 的最大似然估计为样本均值11.55,图2显示了估计的泊松分布和样本观测频率的差异。尽管拥有相同的均值,但由于数据具有泊松低方差特点,其经验分布的比泊松分布集中得多,可以说此时使用泊松分布模型的效果是差强人意的,显然不是一个合适的模型。

1. 加权泊松分布

分别使用 EWP2 和 EWP3 分布对0剂量组数据进行参数的最大似然估计。其中EWP2分布中,$\hat{\lambda}= 11.79$, $\hat{\beta}=0.11$,EWP3中,$\hat{\lambda}=14.56$,$\hat{\beta}_1=-0.15$,$\hat{\beta}_2=0.68$,估计的 EWP2 、EWP3 分布如图3。由于分布具有了泊松低方差特征,通过权重参数$\beta(\beta_1,\beta_2)$将分布向中部“挤压”,分布更加集中且峰值更高,估计的效果相比之下比泊松分布好很多。EWP2 由于在 $\lambda$ 的左右都赋予相等的权重,因此对称的“压缩”模式不如 EWP3 的非对称“压缩”有弹性,从观测数据经验分布的轻微左偏也暗示了允许$\beta_1,\beta_2$取不同值的 EWP3 分布在分布的拟合上更具优势。

图3 EWP2 、EWP3和CBR的分布拟合图

2. CBR分布

CBR参数的极大似然估计$\hat{\lambda}_1,\hat{\lambda}_2,\ldots,\hat{\lambda}_{21}$= {4.66, 7.90, 22.97, 15.02, 18.23, 19.18, 16.49, 14.40, 13.15, 10.01, 9.33, 6.90, 6.32, 6.26, 5.40, 6.28, 3.27, 0.00, 2.00, 2.00, 1.00},图3中显示的 CBR 拟合效果非常好,几乎与经验分布重合。但模型中含有过多的参数使得估计精度大大降低,其中$\hat{\lambda}_3$的估计标准误为12.59,置信区间宽近50。虽然对于拥有多个事件发生率$\hat{\lambda}_k$的 CBR 分布能够灵活地刻画事件发生数之间的任何概率变化,但参数时过多模型的过度拟合是没有太大意义的,也不易于控制和分析。如前所述,可以利用该模型的$\lambda$建立关于k的函数,从而减少模型中的参数个数,在该例子中事件数分类较多时这种做法就显得十分必要。Faddy(2001) 就从众多函数形式中找到了一种能很好地模拟事件数发生概率在最初增长较快而后缓慢下降特点的四参数模型:

$\lambda_k=a(k^be^{-ck}+d), k\geq1$

对于0剂量组数据,$\left(a,b,c,d\right)$的估计值为 {1.360, 3.507, 0.648, 2.953},估计的分布与经验分布的对比如图4,尽管模型中减少了17个参数,但由于函数形式的合理,该分布仍旧保持较好的拟合效果,拟合优度$\chi^2(13)=6.755$,p值为0.914。

图4

针对 2-4-5T 对小鼠胚胎着床影响的实验数据的泊松低方差特性,相对于标准泊松,加权泊松 EWP2、EWP3 和 CBR 在分布上的修正都使模型有许多改进,能够更好地表现出数据中蕴含的现实含义,其生物学意义的解释也更加清晰:在发育毒性的研究中,药物会影响早期的繁殖过程,阻止胚胎的形成。对于多胎的动物(如老鼠),参数起初的增长表示了受精胚胎着床过程可以进一步刺激排卵,之后的下降说明受精卵着床过程的减慢,这段平缓的过程也就使得方差有所降低。

3. 组别分析

之所以选择一个足够合适的分布的重要的意义不仅在于它能较合适地刻画观测数据,更在于它能够精确地刻画不同组别之间的差异。Faddy(2001)表明75剂量组和90剂量组的频数分布相似,并不存在显著性的差异,文章这部分就将比较标准泊松、EWP2、EWP3 和 CBR 这四种分布在检验0剂量组和75/90剂量组之间差异的能力。过程中使用似然比检验。

似然比检验在大样本时具有渐进性。似然比统计量为

$\Lambda(x)=\frac{sup{L(\theta|x):\theta\in\Theta_o}}{sup{L(\theta|x):\theta\in\Theta}}$

当样本量n趋于无穷,$-2log(\Lambda)$将渐进服从$\chi^2(r)$分布,r为参数空间$\Theta$ and $\Theta_o$的维数之差。

表1 四种分布检验效果比较
$-2log(\Lambda)$ 自由度r p值
标准泊松 2.690 1 0.1000
EWP2 4.789 2 0.0912
EWP3 11.980 3 0.0075
CBR 10.764 4 0.0294

标准泊松分布的0剂量组和75/90剂量组的极大对数似然函数值分别为-1837.763和-318.618 ,即负两倍似然比为2.691(自由度为1),p值为0.10(实际值大于0.1),即使在10%的显著性水平上都无法认为0剂量和75/90剂量对小鼠胚胎着床的影响是显著的。EWP2 的负两倍似然比为4.789,p值比标准泊松略小,为0.0912,在10%的显著性水平下可以认为0剂量组和75/90剂量组小鼠胚胎着床的显著差异,但如果显著性水平在5%则无法拒绝原假设。相比之下,CBR和EWP3的负两倍似然比统计量的p值都小得多,在通常5%的显著性水平下能够有力地表明0剂量组和75/90剂量组之间的差异是显著的,且其中 EWP3 的检验效率甚至明显高于 CBR,p值0.0075达到高度显著。

以上检验至少说明在0剂量组和75/90剂量组的比较上 EWP2、EWP3 和 CBR 都优于标准泊松,能够有效地检测出不同组别之间的分布差异,从而证明了本文之初的观点:标准泊松无法准确刻画该实验数据具有泊松低方差的特性,因此将高估剂量组内的方差,在检验上无法有效地识别组间真实存在的差异。如果轻易地使用泊松分布进行分析,将得出0剂量组和75/90剂量组无显著差异的错误结论。而加权泊松分布和 CBR 都在某种程度上克服了标准泊松的缺点,其中 EWP3 和 CBR 则“灵敏”地发现了组间的显著性不同。且 EWP3 能够表现地比 EWP3 出色,还因为剂量组下的频数分布略微左偏,2个加权参数容许 EWP3 更贴切地拟合原始数据的真实分布。

五、结论

当数据出现“超散布性”和“超聚集性”时可能出现的问题,分布假定的错误将分别低估和高估真实数据的方差,从而影响模型的合理性,有时甚至导致得出错误的结论。本文着眼于一类典型的“超聚集性”问题——泊松低方差特性,并针对该类问题的解决归纳了两种方法:泊松加权分布模型和纯生过程分布模型。前者通过对标准泊松分布进行加权修正,克服了泊松分布均值和方差必须相等的局限性,其中 EWP2 和 EWP3 具有形式简单且适用性强的特点,而 EWP3 在很多情况下会优于 EWP2,多一个参数能够较好地模拟较普遍的不对称的单峰经验分布。而纯生过程分布模型在思路上则有很大不同,它基于随机过程中的事件发生机制,对于分类的事件计数数据在理论上有很强的适用性。CBR 能够用足够多的参数$\lambda_k$模拟不同事件数间频率的变化特征,通过建立$\lambda_k$与k的合适的函数形式,可以构造出任何离散分布,尤其适和分析分类较多的数据。此外,$\lambda_k$的函数变化还有利于结合数据在真实世界中的内在的形成原理,对基于不同事件数时刻的时间发生率$\lambda_k$有比较完整的描述,能够赋予参数合理的解释。但 CBR 不适用于分类较少的计数数据,会由于缺少事件产生过程信息它无法表现出其优势。