前言?

生物統(tǒng)計(jì)學(xué)由于其中無法回避的數(shù)學(xué)內(nèi)容，容易讓非數(shù)理或統(tǒng)計(jì)學(xué)背景的制藥業(yè)研發(fā)人員望而生畏。生存分析有三大支柱：(1) Kaplan-Meier曲線；(2) 對(duì)數(shù)秩檢驗(yàn) (Log Rank Test)；(3) Cox回歸模型。其中方法(1)和(2)背后的概念和運(yùn)算除了與條件概率有關(guān)的定義很難科普之外 (但可以訴諸直覺而繞過)，其余內(nèi)容只需中學(xué)數(shù)學(xué)知識(shí)就能理解；方法(3)涉及的數(shù)學(xué)內(nèi)容相對(duì)復(fù)雜，可能更適合讀完本文后意猶未盡的讀者通過教科書來自學(xué)。幸運(yùn)的是，對(duì)于臨床試驗(yàn)中常見的雙樣本分析而言 (也即比較試驗(yàn)組和對(duì)照組的兩條生存曲線時(shí)不考慮其它協(xié)變量和分層分析)，對(duì)數(shù)秩方法既是給出差異顯著性概率值的最常用手段，也能用來近似估算出與Cox回歸法漸近等價(jià) (asymptotically equivalent) 的風(fēng)險(xiǎn)比數(shù)值。

近年來中文出版界出現(xiàn)了包含數(shù)學(xué)公式的專業(yè)科普未必曲高和寡的可喜現(xiàn)象。在這些成功范例的鼓舞下，本文將從最基本的概念開始，對(duì)“生存分析海洋”的兩大水域大膽進(jìn)行深度潛水式科普。

撰文?|?徐亦迅

Robert Weinberg教授的“The Biology of Cancer”是癌癥生物學(xué)領(lǐng)域不可多得的經(jīng)典教科書，自2007年第一版發(fā)行以來，無數(shù)生物專業(yè)學(xué)生以及癌癥研究人員從中獲益。從2011年開始，CTLA-4抗體和PD-1抗體先后成功，令癌癥免疫療法名聲大噪，在全球制藥界變得炙手可熱。Weinberg教授自然要在2013年的第二版中增加不少關(guān)于免疫哨卡抑制劑?(immune checkpoint inhibitors)?的內(nèi)容，在解釋FDA（美國(guó)食品藥品監(jiān)督管理局）為何批準(zhǔn)CTLA-4抗體ipilimumab?(商標(biāo)名：Yervoy)?用于治療已經(jīng)擴(kuò)散的黑色素瘤病人時(shí)，他寫道：“試驗(yàn)組病人接受ipilimumab治療后的中位總生存期是10個(gè)月，而對(duì)照組病人接受gp100多肽疫苗治療后的中位總生存期是6個(gè)月”[1]。言下之意，這4個(gè)月的中位總生存期差異是FDA最后綠燈放行的關(guān)鍵。在這里，Weinberg教授對(duì)臨床研究結(jié)果[2]的解讀可謂“撿了芝麻，丟了西瓜”：他只關(guān)注了“橫切”兩條生存曲線后得到的中位總生存期的差異，而不知原文括號(hào)里的風(fēng)險(xiǎn)比?(HR = 0.66)?和概率值（P = 0.003）才是FDA專家委員會(huì)決策的更重要依據(jù)！要想真正理解Weinberg教授產(chǎn)生誤讀的原因，我們首先需要掌握生存分析的一些基本概念。

生存分析中的基本概念

生存分析?(survival analysis)

生存分析是一系列統(tǒng)計(jì)學(xué)方法的總稱，用以分析研究者感興趣的事件發(fā)生前的持續(xù)時(shí)間數(shù)據(jù)。生存分析的應(yīng)用范圍非常廣泛，其中包括抗癌藥物臨床試驗(yàn)中的生存時(shí)間分析、工程學(xué)中的故障時(shí)間分析、計(jì)量經(jīng)濟(jì)學(xué)中的持續(xù)時(shí)間分析。

事件 (event)

事件的具體定義取決于新藥試驗(yàn)中事先指定的臨床終點(diǎn)?(clinical endpoint)?。如果終點(diǎn)是總體生存期?(overall survival, OS)，那么事件就是患者的死亡。如果終點(diǎn)是無進(jìn)展生存期?(progression-free survival, PFS)，那么事件就是患者病情的進(jìn)展?(例如固體瘤增大或者白血病的血液指標(biāo)惡化)?或者死亡。生存分析能在生物醫(yī)學(xué)以外的許多不同領(lǐng)域有用武之地[3]，其關(guān)鍵就在于事件這個(gè)概念在定義上的靈活性。

生存時(shí)間 (survival time)

指從病人被隨機(jī)分派進(jìn)入臨床試驗(yàn)的分組?(相對(duì)的時(shí)間零點(diǎn))?直到事件發(fā)生所經(jīng)歷的時(shí)間跨度。臨床試驗(yàn)的病人招募通常是個(gè)持續(xù)的過程，不同病人的試驗(yàn)一般始于日歷上不同的具體時(shí)間點(diǎn)，在數(shù)據(jù)分析時(shí)只有采用相對(duì)時(shí)間，才能有同樣的時(shí)間軸及零點(diǎn)。對(duì)于一個(gè)臨床試驗(yàn)的病人群體而言，個(gè)體病人的生存時(shí)間是一個(gè)隨機(jī)變量，用大寫的T表示。而生存曲線橫坐標(biāo)則對(duì)應(yīng)各病人事件發(fā)生的時(shí)間點(diǎn)，它不是隨機(jī)變量?(而用做函數(shù)的自變量)，用小寫的t表示。隨機(jī)變量T一般不遵從正態(tài)分布。

刪失（censoring）

指由于事件沒有被觀測(cè)到或者無法觀測(cè)到，而導(dǎo)致生存時(shí)間無法精確記錄的情況。其中最為常見的情形稱為右刪失（right censoring，圖1），對(duì)這樣的病人我們只知道其生存時(shí)間要大于從試驗(yàn)開始到刪失發(fā)生的時(shí)間。有多種原因可以導(dǎo)致右刪失情況的出現(xiàn)，其中包括：(1)病人在某時(shí)間點(diǎn)上退出試驗(yàn)或失去隨訪信息；(2)病人在整個(gè)試驗(yàn)結(jié)束時(shí)事件還未發(fā)生；(3)病人由于毒性等原因停用被分派的藥物或換用其它藥物；(4)競(jìng)爭(zhēng)風(fēng)險(xiǎn)事件的發(fā)生。

圖1. 生存時(shí)間數(shù)據(jù)的右刪失現(xiàn)象。

生存時(shí)間T的非正態(tài)分布以及刪失情況的存在讓傳統(tǒng)的統(tǒng)計(jì)方法在分析這類數(shù)據(jù)時(shí)無用武之地，于是統(tǒng)計(jì)學(xué)家們殫精竭慮，直到1970年代才使生存分析這一方法體系趨于成熟。

生存函數(shù) (survival function)

又叫累積生存概率（cumulative survival probability），其定義公式為S(t) = P(T > t)。該函數(shù)的意義就是生存時(shí)間大于時(shí)間點(diǎn)t的概率。顯然當(dāng)t = 0時(shí)，S(t) = 1，每個(gè)病人在各自的試驗(yàn)開始時(shí)都還沒有事件發(fā)生。而隨著t的增加，S(t)逐漸向0方向遞減（嚴(yán)格來說是不增）。理論上的T和t一般是連續(xù)變量，相應(yīng)的生存函數(shù)S(t)的圖像就是一條光滑的曲線?(圖2左)。而在臨床試驗(yàn)的實(shí)踐中，T只能在有限個(gè)不同的時(shí)間點(diǎn)t被觀測(cè)到，因此我們常用不增的階梯曲線來描述估測(cè)出的生存函數(shù)（圖2右）。

圖2.理論與實(shí)踐中的生存函數(shù)。丨來源：Kleinbaum, D.G. & Klein, M. (2012) Survival Analysis: A Self-Learning Text, 3rd Edition。

圖3可以幫助我們進(jìn)一步理解光滑曲線和階梯曲線之間的區(qū)別：紅色數(shù)據(jù)點(diǎn)的縱坐標(biāo)代表根據(jù)試驗(yàn)結(jié)果估算出的累積生存概率值，我們?nèi)粝胗霉饣€來連接就需要對(duì)隨機(jī)變量T的分布做出假設(shè)的參數(shù)擬合法，而曲線一般不宜正好經(jīng)過所有的紅點(diǎn)?(那樣會(huì)導(dǎo)致過度擬合而使得統(tǒng)計(jì)模型沒有多大效用)?；若用非參數(shù)的階梯函數(shù)來連接，那么曲線簡(jiǎn)單而唯一確定！這一區(qū)別正是Kaplan-Meier非參數(shù)估計(jì)法能夠成為生存分析領(lǐng)域一大支柱的關(guān)鍵。

圖3. 階梯函數(shù)的使用可以改善光滑曲線估計(jì)生存函數(shù)的過擬合問題。

風(fēng)險(xiǎn)函數(shù) (hazard function) 和風(fēng)險(xiǎn)比 (hazard ratio)

一個(gè)函數(shù)將自變量定義為時(shí)間，目的就是要研究因變量隨時(shí)間的變化規(guī)律，也就離不開“變化率”這個(gè)重要概念，需要用到微分或差分的數(shù)學(xué)工具。生存函數(shù)S(t)描述的是組內(nèi)尚未發(fā)生事件的累積病人比例，由于刪失情況的存在，它與已發(fā)生事件的累積病人比例之間的數(shù)學(xué)關(guān)系不再簡(jiǎn)單。統(tǒng)計(jì)學(xué)家們通過實(shí)踐發(fā)現(xiàn)，研究與病人發(fā)生事件的概率有關(guān)的風(fēng)險(xiǎn)函數(shù)能提取出更多蘊(yùn)含在生存數(shù)據(jù)中的信息。

風(fēng)險(xiǎn)函數(shù)?(又稱為風(fēng)險(xiǎn)率，hazard rate)?的直接數(shù)學(xué)定義需要用到條件概率和極限兩種概念的組合，不是科普文章的最佳選擇。不過該函數(shù)也可以通過生存函數(shù)S(t)來間接定義：

從定義式可以看出，h(t)基本就是S(t)對(duì)時(shí)間的導(dǎo)數(shù)除以S(t)。由于S(t)一般隨時(shí)間遞減，加上一個(gè)負(fù)號(hào)是可以讓h(t)在取值上非負(fù)，從而在使用上更加方便。而h(t)在從零開始的一個(gè)時(shí)間段內(nèi)的積分，H(t)，則被稱為累積風(fēng)險(xiǎn)函數(shù)?(cumulative hazard function)：

以上兩個(gè)定義式為了簡(jiǎn)潔而假定時(shí)間t為連續(xù)變量，若t為離散變量時(shí)也有相對(duì)應(yīng)的定義，我們只需記得微分的離散對(duì)應(yīng)是差分，而積分的離散對(duì)應(yīng)是求和。

在新藥臨床試驗(yàn)中，兩條不同生存曲線的背后隱藏著兩個(gè)不同的風(fēng)險(xiǎn)函數(shù)：h1(t)和h2(t)。統(tǒng)計(jì)學(xué)家們通過研究許多生存數(shù)據(jù)后發(fā)現(xiàn)，如果試驗(yàn)組和對(duì)照組病人的隨機(jī)分派做得足夠好， 我們?cè)诤芏嗲闆r下就能用一個(gè)近似成立的比例風(fēng)險(xiǎn)假設(shè)?(proportional hazards)?來簡(jiǎn)化生存數(shù)據(jù)的分析：

在這里我們假設(shè)被比較的兩個(gè)風(fēng)險(xiǎn)函數(shù)之比率是一個(gè)不隨時(shí)間變化的常數(shù)，用HR?(hazard ratio，風(fēng)險(xiǎn)比)?來表示。制藥業(yè)的常用慣例是：試驗(yàn)組的風(fēng)險(xiǎn)函數(shù)作為分子，而對(duì)照組的風(fēng)險(xiǎn)函數(shù)作為分母?[注意：某些論文和教科書可能會(huì)反過來]。如果試驗(yàn)組生存曲線總體在對(duì)照組之上，那么HR就在0到1之間取值，數(shù)值越小，兩條曲線的差別就越大，也就意味著新藥的相對(duì)療效優(yōu)勢(shì)越大。需要指出的是，雖然生存分析中的常用方法無需比例風(fēng)險(xiǎn)假設(shè)嚴(yán)格成立，但是試驗(yàn)數(shù)據(jù)有時(shí)會(huì)有大幅度偏離該假設(shè)的情況，其統(tǒng)計(jì)分析就需要對(duì)本文所介紹的方法做一些在數(shù)學(xué)上比較復(fù)雜的改進(jìn)，這已經(jīng)超出了本文的討論范圍。

揭秘Kaplan-Meier曲線

可以毫不夸張地說，近幾十年來發(fā)表的臨床癌癥學(xué)會(huì)議報(bào)告和期刊論文中，幾乎每篇都至少要包含一幅Kaplan-Meier生存曲線圖。沒有統(tǒng)計(jì)學(xué)背景的新藥研發(fā)人員在專家的指引下通過多看雖能大致明白生存曲線的含義，但是對(duì)其背后的簡(jiǎn)單運(yùn)算往往有一種不必要的神秘感。而統(tǒng)計(jì)學(xué)家們賦予Kaplan-Meier法的專業(yè)術(shù)語“乘積極限估計(jì)”?(product limit estimator)?更是讓非專業(yè)人士望而卻步。

若想揭開籠罩在Kaplan-Meier曲線上的神秘面紗，最佳途徑似乎是“解剖麻雀法”：通過徹底展示一個(gè)非常簡(jiǎn)單實(shí)例的詳細(xì)計(jì)算過程，我們就可以讓非專業(yè)讀者也能用手畫出一條生存曲線。筆者在此改編了文獻(xiàn)[4]中的一組只有10位病人的生存數(shù)據(jù)，用字母依次給病人編號(hào)后列表如下：

表1. 臨床試驗(yàn)中的生存數(shù)據(jù)示例

其中第二列是時(shí)間點(diǎn)t，以月為單位。而第三列的病人狀態(tài)只有兩種取值：0代表右刪失，1代表死亡事件。整個(gè)試驗(yàn)持續(xù)的最長(zhǎng)相對(duì)時(shí)間是12個(gè)月，期間10人中有7人發(fā)生死亡事件。病人I在第7月出現(xiàn)刪失?(退出試驗(yàn)或失去隨訪)，我們對(duì)其生存時(shí)間的信息不完整，只知道大于7個(gè)月。而病人D在試驗(yàn)結(jié)束時(shí)依然存活。

要從表1估計(jì)出生存曲線，我們首先需將各行按照時(shí)間點(diǎn)t從小到大排序，記為ti，其中下標(biāo)i表示序號(hào)。由于發(fā)生死亡事件與刪失的離散時(shí)間點(diǎn)一般不是均勻分布，這些點(diǎn)把時(shí)間軸分割成了長(zhǎng)短不一的區(qū)間，每個(gè)ti值代表了區(qū)間的起點(diǎn)。然后我們把每個(gè)區(qū)間里發(fā)生的死亡事件累加，記為di。若同時(shí)出現(xiàn)了死亡和刪失，我們將用獨(dú)立的兩行來表示該區(qū)間。例如第7月就對(duì)應(yīng)兩個(gè)ti值：t3?= 7，t4?= 7+，其中的“+”號(hào)代表刪失。每個(gè)區(qū)間起點(diǎn)時(shí)的存活病人數(shù)記為ni，那么在該時(shí)間段的存活比例?(surviving proportion)?就是：pi?= (ni?– di) / ni?。這樣，我們通過整理和簡(jiǎn)單計(jì)算就得到了表2：

表2. Kaplan-Meier法估計(jì)生存函數(shù)的計(jì)算過程

為了使讀者概念更加清晰，筆者特地添加了第一行的時(shí)間零點(diǎn)。如果不存在刪失，表2中每行的生存概率顯然可以用一個(gè)簡(jiǎn)單的比值公式來估算：

但是刪失情況的出現(xiàn)會(huì)給這個(gè)思路帶來困難。例如第8月結(jié)束時(shí)，病人F的死亡讓組內(nèi)存活人數(shù)下降為4人，用4 / 10 = 0.4來估算該時(shí)間段的生存概率顯然不合理，它把第7月病人I的刪失計(jì)為死亡，因此這里的0.4在統(tǒng)計(jì)學(xué)上被稱為有偏估計(jì)量?(biased estimator)?。

要想找到生存概率的無偏估計(jì)量?(unbiased estimator)?并非易事。Paul Meier和Edward Kaplan這兩個(gè)普林斯頓大學(xué)數(shù)學(xué)系博士生在1940年代末跟隨同一位導(dǎo)師(統(tǒng)計(jì)學(xué)巨匠John Tukey)，而且同在1951年完成博士論文答辯，兩人在校期間居然互不相識(shí)[5]。畢業(yè)離校各奔東西后，他們又在互不知曉的情況下獨(dú)立研究了同樣的問題，卻殊途同歸地找到了類似答案。最后在Tukey教授和期刊編輯的協(xié)調(diào)下，花了4年左右的時(shí)間才將兩篇風(fēng)格迥異的手稿整合成一篇論文[6]，著名的Kaplan-Meier乘積極限法終于在難產(chǎn)后誕生。此文問世半個(gè)多世紀(jì)以來，累計(jì)引用次數(shù)已經(jīng)超過6萬，這個(gè)數(shù)字在任何一個(gè)自然科學(xué)領(lǐng)域都是驚人的！

Kaplan和Meier基于條件概率的概念，經(jīng)過一番并不簡(jiǎn)單的數(shù)學(xué)推導(dǎo)和探索，發(fā)現(xiàn)表2中pi的連乘積作為生存概率 S(t) 的非參數(shù)估計(jì)非常合適：

需要注意的是，出現(xiàn)刪失的時(shí)間點(diǎn)沒有對(duì)應(yīng)的pi值，這些項(xiàng)也就自動(dòng)被排除在算式之外?(本例中的p4, p7, p9都不存在)。這個(gè)無偏估計(jì)生存概率的連乘積公式，可以通過一個(gè)類比來獲得直覺上的理解：很多電子游戲都采用序列沖關(guān)的設(shè)計(jì)，游戲者只有先“活過”前面每一關(guān)，才能到達(dá)目前要闖的這一關(guān)。這其實(shí)就是條件概率連乘的思路，與物理學(xué)中串聯(lián)電路的乘法原理也一脈相通。因此所謂“乘積極限法”其實(shí)一點(diǎn)也不神秘?，F(xiàn)在我們可用該公式來重新估算第8月時(shí)間段的生存概率：

0.48這個(gè)數(shù)值顯然比把在前一時(shí)間段發(fā)生刪失的病人I計(jì)為死亡而得出的0.4更為合理?(也要比把病人I當(dāng)成肯定在第8月依然存活的有偏估計(jì)值0.5更合理)?。通過Kaplan-Meier法計(jì)算出的各時(shí)間段內(nèi)生存概率值放在表2的最右一列，然后將這些數(shù)值作為相應(yīng)區(qū)間的縱坐標(biāo)對(duì)時(shí)間作圖，就可得到下面的生存曲線：

圖4. 由表1數(shù)據(jù)計(jì)算得出的Kaplan-Meier生存曲線。

圖4只有6個(gè)數(shù)據(jù)點(diǎn)，有興趣的讀者完全可以用手工來完成。不過用專業(yè)軟件來畫?(例如本文使用的R語言)?可以更方便地在圖中添加注解性內(nèi)容，例如在縱坐標(biāo)為0.5處用一條水平線來就橫切該曲線，就發(fā)現(xiàn)本試驗(yàn)組的中位生存期為8個(gè)月。

通過研究這個(gè)只有10位病人的簡(jiǎn)單實(shí)例，沒有統(tǒng)計(jì)學(xué)背景的讀者跟隨筆者一路走來，或許也能對(duì)Kaplan-Meier曲線有“會(huì)心不遠(yuǎn)”的感覺。

統(tǒng)計(jì)學(xué)中的假設(shè)檢驗(yàn)

對(duì)于新藥的臨床試驗(yàn)而言，最重要的問題顯然是比較試驗(yàn)組和對(duì)照組這兩條生存曲線的差異：曲線分開的程度以及相應(yīng)的統(tǒng)計(jì)顯著性。衡量統(tǒng)計(jì)顯著性的常用方法是對(duì)數(shù)秩檢驗(yàn)，這個(gè)命名起源于非?；逎牡葍r(jià)數(shù)學(xué)推導(dǎo)[7]，我們?cè)趯?shí)際應(yīng)用中根本不會(huì)感覺到“對(duì)數(shù)”和“秩”的存在。一個(gè)更好的名稱也許是Mantel檢驗(yàn)，可惜文獻(xiàn)中的習(xí)慣一旦形成就很難更改。

用概率值來度量顯著性是基礎(chǔ)統(tǒng)計(jì)學(xué)的核心支柱之一，其標(biāo)準(zhǔn)名稱為假設(shè)檢驗(yàn)(hypothesis testing)，這是一個(gè)幾乎人人有所耳聞卻又經(jīng)常被誤解的論題。統(tǒng)計(jì)學(xué)界內(nèi)部都沒能統(tǒng)一思想，在假設(shè)檢驗(yàn)如何操作上分為兩大流派：Fisher流派是大多數(shù)教科書中的主流，但是Neyman-Pearson流派的影響也不容忽視。如何正確理解假設(shè)檢驗(yàn)是大多數(shù)非統(tǒng)計(jì)學(xué)背景的讀者在自學(xué)中需要跨越的路障之一。2016年3月，英國(guó)著名的《Nature》雜志發(fā)表了一篇點(diǎn)擊率很高的新聞報(bào)道，其中傳達(dá)的重要信息發(fā)人深?。好绹?guó)統(tǒng)計(jì)學(xué)協(xié)會(huì)對(duì)多年來科學(xué)文獻(xiàn)中存在的大量對(duì)檢驗(yàn)概率值的誤用和濫用深表擔(dān)憂[8]。

一言以蔽之，假設(shè)檢驗(yàn)就是一個(gè)統(tǒng)計(jì)推斷過程，它能基于樣本數(shù)據(jù)?(sample)?而在兩個(gè)相互對(duì)立的描述總體?(population)?特性的論斷中做出合理選擇。其中一個(gè)我們稱為零假設(shè)?(null hypothesis)，而與它對(duì)立的就稱為備擇假設(shè)?(alternative hypothesis)。零假設(shè)的設(shè)立原則可以通過與一個(gè)現(xiàn)實(shí)生活實(shí)例的類比來說明：刑事訴訟中的無罪推定原則。刑事法庭的“零假設(shè)”就是被告無罪，只有在公訴方取得足夠的可靠證據(jù)?(“樣本數(shù)據(jù)”)?可以超越合理懷疑?(“統(tǒng)計(jì)顯著性閾值”)?之后，審判方才能采納被告有罪這一“備擇假設(shè)”。兩個(gè)對(duì)立的假設(shè)確定后，我們還需找到一個(gè)可從觀測(cè)數(shù)據(jù)推算的統(tǒng)計(jì)量（test statistic），其零分布（就是在零假設(shè)成立前提下的概率分布）是已知的，從而計(jì)算出得到該統(tǒng)計(jì)量比從來自樣本的觀測(cè)值更極端的概率值。如果概率值小于事先約定的顯著性閾值（比如常用的0.05），我們就拒絕零假設(shè)而接受備擇假設(shè)；而若概率值大于顯著性閾值，我們就無法拒絕零假設(shè)。值得注意的是，當(dāng)零假設(shè)不能被拒絕時(shí)，我們并不能證明零假設(shè)成立！另外概率值越小雖然說明統(tǒng)計(jì)顯著性越強(qiáng)，但它并不能用來描述兩組之間的差異大小，這是由于概率值還依賴于樣本量。若對(duì)這兩個(gè)概念要點(diǎn)沒有深刻領(lǐng)悟，非統(tǒng)計(jì)學(xué)專業(yè)背景的研究人員在陳述實(shí)驗(yàn)結(jié)果時(shí)就容易出錯(cuò)。

孟德爾遺傳學(xué)中的擬合優(yōu)度檢驗(yàn)

生存分析中的對(duì)數(shù)秩檢驗(yàn)屬于經(jīng)典統(tǒng)計(jì)學(xué)中著名的卡方擬合優(yōu)度檢驗(yàn)?(chi-squared goodness-of-fit test)?范疇，而擬合優(yōu)度檢驗(yàn)在生物學(xué)領(lǐng)域最著名的例子就是統(tǒng)計(jì)學(xué)泰斗Ronald Fisher將其用于分析孟德爾?(Gregor Mendel)?的遺傳學(xué)實(shí)驗(yàn)數(shù)據(jù)[9]。

圖5. 孟德爾根據(jù)豌豆雜交實(shí)驗(yàn)提出的性狀分離理論[9]

圖5展示了孟德爾根據(jù)豌豆雜交實(shí)驗(yàn)結(jié)果提出的性狀分離假設(shè)。他在該實(shí)驗(yàn)中用了兩個(gè)不同性狀的親代植株?(P1)：父本是高豌豆株?(基因型記為TT)，母本是矮豌豆株?(基因型記為tt)，其中高性狀的等位基因T對(duì)于矮性狀的等位基因t顯性?；诿系聽柤僭O(shè)的理論推測(cè)是：子一代?(F1)?自花受粉后得到的子二代?(F2)?表現(xiàn)型將呈現(xiàn)3:1的分離。但是生命現(xiàn)象內(nèi)在的隨機(jī)性需要我們用概率論的視角來看問題，實(shí)驗(yàn)結(jié)果一般會(huì)和理論預(yù)期的3:1有一定的差異。孟德爾有一次實(shí)驗(yàn)后得到的子二代中有787株為高，277株為矮，觀測(cè)到的比值為2.84:1。我們可以通過擬合優(yōu)度檢驗(yàn)來判斷實(shí)測(cè)的高株數(shù)和矮株數(shù)是否與孟德爾假設(shè)相容?(注意：卡方擬合優(yōu)度檢驗(yàn)只能用于計(jì)數(shù)數(shù)據(jù)，而不能用于非整數(shù)的比值數(shù)據(jù))?，這里兩條對(duì)立的假設(shè)分別是：

• 零假設(shè)(H0)：觀測(cè)到的高矮株計(jì)數(shù)符合理論預(yù)期的3:1比例；

• 備擇假設(shè)(H1)：觀測(cè)到的高矮株計(jì)數(shù)不符合理論預(yù)期3:1比例。

要定義合適的統(tǒng)計(jì)量絕非易事，近代統(tǒng)計(jì)學(xué)先驅(qū)Karl Pearson經(jīng)過一番探索，率先在1900年找到了著名的擬合優(yōu)度統(tǒng)計(jì)量，其定義公式為：

其中Oi代表觀測(cè)到的高和矮兩類?(n=2)?性狀計(jì)數(shù)，Ei代表理論預(yù)期的兩類性狀計(jì)數(shù)。該統(tǒng)計(jì)量在零假設(shè)前提下的精確概率分布很難確定，好在樣本足夠大時(shí)非常接近著名的卡方分布?(chi-squared distribution)??？ǚ椒植记€的形狀還取決于自由度?(degrees of freedom)?，這個(gè)重要概念在很多教科書中都沒有透徹的解釋。本文也只能做簡(jiǎn)要講解，有興趣深入鉆研的讀者可以參閱文獻(xiàn)[10]。概率論告訴我們，k個(gè)標(biāo)準(zhǔn)正態(tài)分布的獨(dú)立隨機(jī)變量之平方和遵循自由度為k的卡方分布?？ǚ浇y(tǒng)計(jì)量是一個(gè)由多個(gè)平方和累加而得的非負(fù)實(shí)數(shù)，總類數(shù)n越大，相應(yīng)的項(xiàng)數(shù)也越多而導(dǎo)致統(tǒng)計(jì)量數(shù)值增大。由類數(shù)變多引起的統(tǒng)計(jì)量增值顯然與假設(shè)檢驗(yàn)的顯著性無關(guān)，因此我們有必要對(duì)卡方統(tǒng)計(jì)量的分布曲線按照n的大小細(xì)分。在本文涉及的擬合優(yōu)度檢驗(yàn)和對(duì)數(shù)秩檢驗(yàn)中，自由度都是k = n – 1。孟德爾實(shí)驗(yàn)中的n為2，因此其檢驗(yàn)統(tǒng)計(jì)量的近似分布就是圖6中自由度為1的棕黃色曲線：

圖6. 不同自由度下卡方分布的概率密度曲線。丨來源：Wikipedia

根據(jù)上述公式，我們可以基于孟德爾的實(shí)驗(yàn)數(shù)據(jù)通過下表來計(jì)算擬合優(yōu)度檢驗(yàn)所需的統(tǒng)計(jì)量：

表3. 根據(jù)孟德爾豌豆雜交實(shí)驗(yàn)數(shù)據(jù)推算的統(tǒng)計(jì)量

表3中的798和266分別是總數(shù)1064的3/4和1/4，也就是理論預(yù)期值?(Ei)。最后我們得到的統(tǒng)計(jì)量數(shù)值為0.60。通過查閱自由度為1時(shí)的卡方分布累積概率表，我們知道對(duì)應(yīng)于概率值p=0.05的統(tǒng)計(jì)量數(shù)值為3.841。0.60要比3.841小很多，這意味著擬合優(yōu)度檢驗(yàn)的概率值要遠(yuǎn)大于0.05?(用一條簡(jiǎn)單的R語言代碼可以算出p=0.439)，因此我們無法拒絕零假設(shè)，而認(rèn)為孟德爾的實(shí)驗(yàn)數(shù)據(jù)與他理論預(yù)期的3:1分離比例相容。需要再次強(qiáng)調(diào)的是，實(shí)驗(yàn)數(shù)據(jù)與零假設(shè)的相容性并不足以證明零假設(shè)本身。

比較兩條生存曲線

以上對(duì)孟德爾遺傳學(xué)中的擬合優(yōu)度檢驗(yàn)之回顧將有助于我們深度理解對(duì)數(shù)秩檢驗(yàn)何以能確定兩條生存曲線差異的統(tǒng)計(jì)顯著性。我們將通過一個(gè)來自文獻(xiàn)[11]的實(shí)例來闡明對(duì)數(shù)秩統(tǒng)計(jì)量的計(jì)算過程。這是一個(gè)總共54位病人的急性淋巴性白血病?(acute lymphoblastic leukemia，ALL)?臨床試驗(yàn)，目的是比較異體?(allogenic)?骨髓移植和自體?(autologous)?骨髓移植這兩種療法的差異。對(duì)ALL的一線治療需要高強(qiáng)度的化療和放療，但這同時(shí)也會(huì)損傷成年患者的骨髓造血系統(tǒng)。很多ALL病人在一線治療后還需要通過骨髓移植來修復(fù)造血功能和免疫系統(tǒng)。最理想的骨髓來自配型成功的病人親屬，這種異體骨髓移植的預(yù)后效果最佳。若病人沒有合適的異體骨髓來源，就只能退而求其次，將病人自身的骨髓抽出體外進(jìn)行藥物處理，在殺滅ALL惡性細(xì)胞后再移植回病人體內(nèi)。臨床實(shí)踐表明，在自體骨髓移植過程中，惡性細(xì)胞很難被消滅干凈，因此其預(yù)后效果一般不如異體移植。我們對(duì)該臨床試驗(yàn)的結(jié)果用上一節(jié)介紹的時(shí)間點(diǎn)排序后，整理成下面的表格：

表4. 急性淋巴性白血病臨床試驗(yàn)的生存數(shù)據(jù)[11, 4]

我們對(duì)這些數(shù)據(jù)進(jìn)行前文已經(jīng)介紹過的Kaplan-Meier乘積極限運(yùn)算，就可畫出這樣兩條生存曲線：

圖7. 通過對(duì)表4數(shù)據(jù)計(jì)算得出的Kaplan-Meier生存曲線。

從圖7可以看出，接受異體骨髓移植的ALL病人生存狀況要明顯好于自體骨髓移植。但是這個(gè)試驗(yàn)的樣本不大，隨機(jī)取樣引起的波動(dòng)依然可以讓我們對(duì)兩條生存曲線差異的顯著性存疑。對(duì)數(shù)秩檢驗(yàn)可以幫助我們來定量統(tǒng)計(jì)顯著性，它的零假設(shè)是：兩條生存曲線代表的病人群體在任意時(shí)間點(diǎn)上的生存概率沒有差別。與零假設(shè)對(duì)立的備擇假設(shè)無需贅述。需要特別指出的是：與Cox回歸模型一樣，對(duì)數(shù)秩檢驗(yàn)也依賴于比例風(fēng)險(xiǎn)假設(shè)的近似成立。

對(duì)數(shù)秩檢驗(yàn)的統(tǒng)計(jì)量與前述卡方擬合優(yōu)度統(tǒng)計(jì)量在概念上頗有相似之處，只是由于刪失數(shù)據(jù)點(diǎn)的存在，其計(jì)算過程要相對(duì)復(fù)雜一些。與Kaplan-Meier曲線的計(jì)算相似，第一步也是將兩組病人發(fā)生死亡事件的各時(shí)間點(diǎn)t從小到大排序，每個(gè)ti值代表時(shí)段的起點(diǎn)。我們用A代表自體骨髓移植，用B代表異體骨髓移植，并將表4重新整理以清晰呈現(xiàn)詳細(xì)的計(jì)算過程：

表5. 對(duì)表4生存數(shù)據(jù)進(jìn)行對(duì)數(shù)秩檢驗(yàn)的計(jì)算過程

表5中的變量多數(shù)都有兩個(gè)下標(biāo)：第一個(gè)下標(biāo)代表療法的分組，A或者B。下標(biāo)符號(hào)為T時(shí)表示同一行中A組和B組計(jì)數(shù)相加；第二個(gè)下標(biāo)i還是時(shí)間點(diǎn)序號(hào)。變量O代表時(shí)段內(nèi)發(fā)生死亡事件的人數(shù)，變量r代表時(shí)段起點(diǎn)的存活人數(shù)，變量f代表時(shí)段內(nèi)兩組累計(jì)的死亡比例。r與f的乘積就是時(shí)段內(nèi)的預(yù)期死亡人數(shù)，用E變量來表示，這和卡方擬合優(yōu)度檢驗(yàn)中的E變量是同一個(gè)概念。重要的是，發(fā)生刪失的時(shí)間點(diǎn)由于對(duì)E變量的計(jì)算沒有直接貢獻(xiàn)而無需列入表5，但是相應(yīng)的刪失病人數(shù)需要從下一時(shí)段的r值中扣除。

對(duì)于熟練的臨床生統(tǒng)專家而言，瑣碎的計(jì)算過程早就被隱藏在簡(jiǎn)單的幾行代碼背后而由電腦來完成。但對(duì)于生存分析的初學(xué)者而言，若能在閱讀本文時(shí)拿起紙筆對(duì)表5中行間與列間的聯(lián)系進(jìn)行簡(jiǎn)單驗(yàn)算，將會(huì)有助于徹底掌握對(duì)數(shù)秩檢驗(yàn)的神髓。算出各時(shí)段內(nèi)E變量的數(shù)值后，我們將表5中的四列數(shù)字累加，就得到了與前文孟德爾豌豆雜交實(shí)驗(yàn)結(jié)果的卡方擬合優(yōu)度檢驗(yàn)非常相似的四個(gè)數(shù)值：OA，OB，EA，EB。對(duì)數(shù)秩統(tǒng)計(jì)量最簡(jiǎn)明的計(jì)算方法可以采用下面的卡方公式：

該統(tǒng)計(jì)量在零假設(shè)前提下也近似遵循自由度為1的卡方分布。把相應(yīng)的數(shù)值代入公式，我們得到的結(jié)果是：2.22 + 2.96 = 5.18。運(yùn)用統(tǒng)計(jì)學(xué)軟件或者查卡方分布的累積概率表，可以算出當(dāng)統(tǒng)計(jì)量為5.18時(shí)的對(duì)數(shù)秩檢驗(yàn)概率值：p = 0.023 < 0.05。因此試驗(yàn)數(shù)據(jù)提供了足夠的證據(jù)讓我們拒絕零假設(shè)而采納備擇假設(shè)：異體移植組和自體移植組的兩條生存曲線有統(tǒng)計(jì)顯著性差異。有興趣深究對(duì)數(shù)秩檢驗(yàn)的讀者，不妨先驗(yàn)證OA?+ OB?= EA?+ EB。隨后可以研讀生存分析的常用教科書?(比如圖2來源的那本)，掌握另一種運(yùn)用方差的對(duì)數(shù)秩統(tǒng)計(jì)量估算法。由此通過筆算或?qū)I(yè)軟件得到的卡方值為5.48，通常與本文采用的算法結(jié)果相似。

僅僅確立兩種療法差異的統(tǒng)計(jì)顯著性?(statistical significance)?還不能讓臨床醫(yī)生們完全滿意，他們還想知道這個(gè)差異的幅度有多大。如果我們對(duì)圖7中的兩條生存曲線“橫切”與“縱切”，雖然能提取出反映臨床顯著性?(clinical significance)?的一些數(shù)值，但是它們只是生存曲線“電影”的幾個(gè)“截屏”。例如本例中自體移植組的中位生存期是17個(gè)月，而異體移植組的中位生存期尚未達(dá)到?(not yet reached, NYR)，該“截屏”說明異體骨髓移植療法在臨床上有優(yōu)勢(shì)，但卻不能在整體上說出這個(gè)優(yōu)勢(shì)有多大。

其實(shí)我們?cè)诨靖拍钜还?jié)中定義的風(fēng)險(xiǎn)比?(HR)?才是衡量臨床顯著性的更重要指標(biāo)，它的數(shù)值體現(xiàn)了兩條曲線的整體差異度。雖然Cox回歸分析是臨床生統(tǒng)實(shí)踐中更為常用的估算HR方法，但是對(duì)數(shù)秩方法也能用一個(gè)簡(jiǎn)單公式來近似估算HR?[12]，它與用Cox模型算出的HR在樣本數(shù)足夠大的漸近意義上等價(jià)：

將前文算出的四個(gè)數(shù)值代入該公式，我們得到的HR數(shù)值為0.41，也就是說異體移植組病人的死亡風(fēng)險(xiǎn)比自體移植組病人下降了59%左右，這無疑體現(xiàn)了異體骨髓移植的巨大臨床優(yōu)勢(shì)。作為比較，用R語言的Cox回歸分析算出的HR數(shù)值是0.39，與對(duì)數(shù)秩方法的估算非常接近。HR數(shù)值代表的相對(duì)死亡風(fēng)險(xiǎn)下降程度雖然更加重要，但也不宜把它看成是衡量臨床顯著性的唯一指標(biāo)，對(duì)生存曲線的多點(diǎn)橫切與縱切可以起到輔助作用，從而為臨床醫(yī)學(xué)工作者提供一幅完整的畫面。

致謝：在修訂本文的2016年初稿過程中，北京天壇醫(yī)院谷鴻秋博士、美國(guó)Karyopharm Therapeutics公司湯詩杰博士、美國(guó)斯坦福大學(xué)統(tǒng)計(jì)系田魯教授先后提供了有益的建議，筆者在此特別感謝。

參考文獻(xiàn)

[1]?Weinberg, R.A. (2013) Chapter 15, The Biology of Cancer, 2nd ed. Garland Science, Taylor & Francis Group, LLC.

[2]?Hodi, F.S. et al. (2010) New England Journal of Medicine 363: 711-723.

[3]?Linoff, G.S. & Berry, M.J.A. (2011) Chapter 10, Data Mining Techniques: For Marketing, Sales, and Customer Relationship Management, 3rd ed. Wiley Publishing, Inc.

[4]?Glantz, S.A. (2012) Chapter 11, Primer of Biostatistics, 7th ed. The McGraw-Hill Companies, Inc.

[5]?Marks, H.M. (2004) Clinical Trials 1: 131-138.

[6]?Kaplan, E.L. & Meier, P. (1958) Journal of the American Statistical Association 53: 457-481.

[7]?Peto, R. & Pike, M.C. (1973) Biometrics 29: 579-584.[8]?Baker, M. (2016) Nature 531: 151.

[9]?Tamarin, R.H. (2001) Chapter 4, Principles of Genetics, 7th ed. The McGraw-Hill Companies, Inc.

[10]?Dallal, G.E. (2012) The Little Handbook of Statistical Practice ?(http://www.jerrydallal.com/LHSP/dof.htm)

[11]?Vey, D. et al. (1994) Bone Marrow Transplantation 14: 383-388.

[12]?Machin, D. & Gardner, M.J. (1988) British Medical Journal 296: 1369-1371.