朱 玉梅 楊李 杰陳佰鋒姚應(yīng)水△

Cox比例風(fēng)險(xiǎn)Frailty模型簡(jiǎn)介與軟件實(shí)現(xiàn)*

朱 玉1梅 楊2李 杰1陳佰鋒1姚應(yīng)水1△

目的介紹Cox比例風(fēng)險(xiǎn)Frailty模型的原理及其在SAS 9.3軟件中的實(shí)現(xiàn)過程。方法利用具體數(shù)據(jù)的分析過程介紹Cox比例風(fēng)險(xiǎn)Frailty模型在SAS 9.3軟件中的實(shí)現(xiàn)，并比較Cox比例風(fēng)險(xiǎn)模型與Cox比例風(fēng)險(xiǎn)Frailty模型的分析效果。結(jié)果Frailty項(xiàng)對(duì)數(shù)變換后的方差估計(jì)值為0.831，與0比較差異有統(tǒng)計(jì)學(xué)意義，有必要在Cox比例風(fēng)險(xiǎn)模型中加入Frailty項(xiàng)。結(jié)論Cox比例風(fēng)險(xiǎn)Frailty模型能夠揭示資料的異質(zhì)性，準(zhǔn)確地分析因素對(duì)結(jié)局變量的影響，獲得更為客觀的分析結(jié)論。

Frailty模型 Cox比例風(fēng)險(xiǎn)模型 Cox比例風(fēng)險(xiǎn)Frailty模型

以時(shí)間到事件（Time-To-Event）為結(jié)局變量的生存分析方法眾多，其中經(jīng)典的分析方法是Cox比例風(fēng)險(xiǎn)模型。它的理論假設(shè)之一是研究對(duì)象間相互獨(dú)立，即研究對(duì)象間具有同質(zhì)性，這暗示每個(gè)研究對(duì)象有相同的基線風(fēng)險(xiǎn)。在實(shí)際研究中，這種假設(shè)不易達(dá)到，由于未知因素或設(shè)計(jì)等原因，導(dǎo)致研究對(duì)象個(gè)體間或組別間存在異質(zhì)性，研究對(duì)象表現(xiàn)出組內(nèi)相關(guān)的特性，經(jīng)典的生存分析不再適合分析這類型的數(shù)據(jù)。這時(shí)，需要對(duì)經(jīng)典的生存分析的方法進(jìn)行改進(jìn)，引入Frailty項(xiàng)（異質(zhì)性變量），發(fā)展為Frailty模型（異質(zhì)性模型）。本文介紹加入Frailty項(xiàng)的Cox比例風(fēng)險(xiǎn)模型——Cox比例風(fēng)險(xiǎn)Frailty模型和軟件實(shí)現(xiàn)。

基本原理

以多中心臨床隨機(jī)化對(duì)照試驗(yàn)為例，介紹Cox比例風(fēng)險(xiǎn)Frailty模型。中心i（i＝1，2，…，k），每個(gè)中心的研究對(duì)象j（j＝1，2，…，ni）。在介紹Cox比例風(fēng)險(xiǎn)Frailty模型之前，先回顧下Cox比例風(fēng)險(xiǎn)模型。

1.Cox比例風(fēng)險(xiǎn)模型

Cox提出比例風(fēng)險(xiǎn)模型［1］如下：

其中h0（t｜X）稱為基線風(fēng)險(xiǎn)函數(shù)，h（t｜X）表示在協(xié)變量為X條件下的風(fēng)險(xiǎn)函數(shù)，exp（β）為風(fēng)險(xiǎn)比（hazard ratio，HR）。Cox比例風(fēng)險(xiǎn)模型最突出的優(yōu)點(diǎn)在于，對(duì)參數(shù)β的估計(jì)并不依賴于基線風(fēng)險(xiǎn)函數(shù)的取值。但是，它要求研究對(duì)象相互獨(dú)立，即研究對(duì)象具有同質(zhì)性。

2.Cox比例風(fēng)險(xiǎn)Frailty模型

當(dāng)分析多中心臨床隨機(jī)化對(duì)照試驗(yàn)的數(shù)據(jù)時(shí)，每個(gè)中心的基線風(fēng)險(xiǎn)可能不一致，需要對(duì)Cox比例風(fēng)險(xiǎn)模型進(jìn)行如下改進(jìn)，發(fā)展為Cox比例風(fēng)險(xiǎn)Frailty模型［2-7］，如下：

其中γi為中心i（i＝1，2，…，k）的隨機(jī)效應(yīng)，其它同Cox比例風(fēng)險(xiǎn)模型。常對(duì)γi做如下假設(shè)：

上式改寫為下式：

其中μi被稱為Frailty項(xiàng)（異質(zhì)性變量），μi是相互獨(dú)立同分布的，其期望是1，方差為未知參數(shù)θ2，通常還假設(shè)其與X相互獨(dú)立，當(dāng)μi＝1時(shí)Cox比例風(fēng)險(xiǎn)Frailty模型蛻變?yōu)镃ox比例風(fēng)險(xiǎn)模型。此外，Elbers和Ridder［8］證明了在Frailty項(xiàng)的期望為1的條件下，模型是可估的。

當(dāng)μi＞1時(shí)說明中心i內(nèi)的研究對(duì)象傾向于加速失效（事件發(fā)生）；當(dāng)μi＜1時(shí)說明中心i內(nèi)的研究對(duì)象傾向于減速失效（事件發(fā)生）；當(dāng)μi＝1時(shí)說明中心i內(nèi)的研究對(duì)象的失效風(fēng)險(xiǎn)（事件發(fā)生）是正常風(fēng)險(xiǎn)。

為了估計(jì)未知參數(shù)，需要對(duì)Frailty項(xiàng)的分布做出假設(shè)。原則上期望為1和具有有限方差的任一正連續(xù)分布都可以作為Frailty項(xiàng)的分布。在應(yīng)用中常采用的分布有伽馬分布、逆高斯分布和對(duì)數(shù)正態(tài)分布等。方差θ2的估計(jì)方法有限制性最大似然估計(jì)（residual maximum likelihood，REML）和最大似然估計(jì)（maximum likelihood，ML）。

實(shí)例分析與軟件實(shí)現(xiàn)

下面結(jié)合SASPROC PHREG過程和SAS help中提供的數(shù)據(jù)介紹Cox比例風(fēng)險(xiǎn)Frailty模型的應(yīng)用和軟件實(shí)現(xiàn)［7］。

1.數(shù)據(jù)簡(jiǎn)介

該數(shù)據(jù)來源于一項(xiàng)關(guān)于糖尿病患者眼睛失明的研究，共包含197例糖尿病患者，他們的眼睛有很高的失明風(fēng)險(xiǎn)，每個(gè)研究對(duì)象的兩只眼睛中隨機(jī)選擇一只眼睛接受激光光凝治療，研究的目的主要是探索激光光凝治療有沒有延緩眼睛失明的進(jìn)程。因?yàn)榍嗌倌旰统赡晏悄虿∮胁煌陌l(fā)展過程，所以同時(shí)檢查發(fā)病年齡和眼睛失明時(shí)間有無關(guān)聯(lián)。每個(gè)患者是一個(gè)“群（clusters）”，類似于上文提到的“中心”，因?yàn)樽笥已劬]有生物學(xué)差異，我們假設(shè)它們有相同的基線風(fēng)險(xiǎn)。變量（id）表示患者編號(hào)，變量（time）是患者眼睛失明時(shí)間，變量（status）表示眼睛失明狀況（是二分類變量，“0”代表沒有失明，“1”代表失明），變量（age）表示患者的發(fā)病年齡（是二分類變量，“0”代表發(fā)病年齡≤20，“1”代表發(fā)病年齡＞20），變量（treatment）表示眼睛接受治療的方法（是二分類變量，“0”代表其他治療，“1”代表激光光凝治療）。用數(shù)據(jù)步建立數(shù)據(jù)集blind，見表1。

表1 建立數(shù)據(jù)集和數(shù)據(jù)分析的程序

2.語句介紹

在SAS 9.3版本的PROC PHREG過程中加入了分析Frailty模型的語句——Random語句。Random語句規(guī)定Frailty項(xiàng)服從對(duì)數(shù)正態(tài)分布，即γi～N（0，σ2），其中σ為未知參數(shù)，需要估計(jì)。Random語句的結(jié)構(gòu)和功能如下：

Vaiable是用來指定“群”，它必須是分類變量，即是class語句中的變量。options有如下：

Abspconv＝r，是規(guī)定Frailty項(xiàng)對(duì)數(shù)變換后方差估計(jì)時(shí)迭代的收斂準(zhǔn)則，如果迭代收斂。

Alpha＝value，是規(guī)定隨機(jī)效應(yīng)的（1-α）可信區(qū)間，默認(rèn)值是0.05。

Method＝reml｜ml，是規(guī)定方差參數(shù)估計(jì)的方法，默認(rèn)方法是rem l。

Noclprint，控制在輸出時(shí)不打印隨機(jī)效應(yīng)的分類信息表。

Pconv＝r，是規(guī)定Frailty項(xiàng)對(duì)數(shù)變換后方差估計(jì)時(shí)迭代的收斂準(zhǔn)則，如果迭代收斂。

Solution，規(guī)定顯示正態(tài)分布的隨機(jī)效應(yīng)估計(jì)值，同時(shí)也顯示對(duì)數(shù)正態(tài)分布的估計(jì)值，即Frailty的估計(jì)值。

Initialvariance｜initial＝value，規(guī)定Frailty項(xiàng)對(duì)數(shù)變換后方差估計(jì)的初始值，默認(rèn)為1。

3.分析與結(jié)果

首先不考慮Frailty項(xiàng)，建立含有變量treatment、age的主效應(yīng)和交互效應(yīng)的Cox比例風(fēng)險(xiǎn)模型，然后模型再加入Frailty項(xiàng)，分析程序見表1。為了方便對(duì)比，把主要結(jié)果列于表2、3。

表2 Cox比例風(fēng)險(xiǎn)模型與Frailty模型的參數(shù)估計(jì)

表3 兩個(gè)模型的激光光凝治療相對(duì)于其他治療的條件風(fēng)險(xiǎn)比估計(jì)

從表2可見，F(xiàn)railty項(xiàng)對(duì)數(shù)變換后的方差估計(jì)值為0.831，與0的差異有統(tǒng)計(jì)學(xué)意義，提示分析時(shí)需要考慮個(gè)體間的異質(zhì)性，有必要在Cox比例風(fēng)險(xiǎn)模型中加入Frailty項(xiàng)。兩個(gè)模型的參數(shù)估計(jì)值相近，Cox比例風(fēng)險(xiǎn)模型比Cox比例風(fēng)險(xiǎn)Frailty模型要保守，在檢驗(yàn)水準(zhǔn)為0.05的情況下，都提示治療方法與發(fā)病年齡間存在交互作用。表3給出了不同發(fā)病年齡下激光光凝治療相對(duì)于其他治療的風(fēng)險(xiǎn)比，從中可看出激光光凝治療方法可延緩眼睛失明進(jìn)程，并且成年組延緩眼睛失明進(jìn)程的效果要好于青少年組。

討 論

本質(zhì)上，F(xiàn)railty模型是在經(jīng)典的生存分析的基礎(chǔ)上加入了一個(gè)隨機(jī)效應(yīng)，以控制研究對(duì)象個(gè)體間或組別間存在的異質(zhì)性。雖然產(chǎn)生異質(zhì)性的變量也可以通過設(shè)置啞變量或用固定效應(yīng)分析，但固定效應(yīng)分析方法不是最優(yōu)的分析方法。按照異質(zhì)性產(chǎn)生的單位不同，把Frailty模型分為個(gè)體Frailty模型、共享Frailty模型兩種。個(gè)體Frailty模型假設(shè)研究對(duì)象個(gè)體之間具有異質(zhì)性，基線風(fēng)險(xiǎn)隨著個(gè)體而改變；共享Frailty模型假設(shè)組別間具有異質(zhì)性，基線風(fēng)險(xiǎn)隨著組別而改變，同一組別內(nèi)的個(gè)體享有相同的基線風(fēng)險(xiǎn)。研究對(duì)象個(gè)體間或組別間存在的異質(zhì)性是不能忽略，有文章也指出［9-10］：如果有異質(zhì)性存在，而沒有考慮，會(huì)導(dǎo)致效應(yīng)低估，這點(diǎn)在本文中也有體現(xiàn)——Cox比例風(fēng)險(xiǎn)模型比Cox比例風(fēng)險(xiǎn)Frailty模型要保守，但是兩者的意義是不同的，考慮異質(zhì)性時(shí)效應(yīng)是條件的效應(yīng)，不考慮異質(zhì)性時(shí)，效應(yīng)是邊際效應(yīng)，在預(yù)測(cè)時(shí)需要注意這點(diǎn)。

與在Cox比例風(fēng)險(xiǎn)模型中引入Frailty項(xiàng)類似，也可以在其它生存分析方法——Weibull回歸模型、加速失效模型中引入Frailty項(xiàng)，擴(kuò)展生存分析的分析方法。

關(guān)于“Frailty模型”的翻譯，有學(xué)者譯為“脆弱模型”，但筆者認(rèn)為有不妥之處。雖然“Frailty”有“脆弱、虛弱等”之意，但是“脆弱模型”所表達(dá)出來的意思不能反映“Frailty模型”的本質(zhì)。結(jié)合“Frailty模型”的產(chǎn)生原因與形成思想，筆者反復(fù)斟酌，認(rèn)為譯為“異質(zhì)性模型”更合適。在沒有規(guī)范的譯名情況下，使用“Frailty模型”是最恰當(dāng)?shù)摹?/p>

1.Cox DR.Regression models and life tables.JR Stat Soc Series B Stat Methodol，1972，34（2）：187-220.

2.M cGilchrist CA，Aisbett CW.Regression with frailty in survival analysis.Biometrics，1991，47：461-466.

3.Therneau TM，Grambsch PM.Modelling survival data：extending the cox model.Springer：New York，2000.

4.Ha ID，Lee Y，MacKenzie G.Model selection for multi-component frailty models.Stat Med.2007，26（26）：4790-807.

5.Andersen PK，Klein JP，Knudsen KM，et al.Estimation of variance in Cox′s regressionmodel with shared gamma frailties.Biometrics，1997，53（4）：1475-84.

6.Elbers C，Ridder G.True and spurious duration dependence：the identifiability of the proportional hazard model.Review of Econom ic Studies，1982，49（3）：403-40.

7.SAS Institute Inc.SAS/STAT?9.3 User′s Guide.Cary，NC：SAS Institute Inc，2011.

8.Elbers C，Ridder G.True and spurious duration dependence：the identifiability of the proportional hazard model.Review of Economic Studies，1982，49（3）：403-9.

9.Gail MH，W ieand S，Piantadosi S.Biased estimates of treatment effect in random ized experimentsw ith nonlinear regressions and om itted covariates.Biometrika，1984，71（3）：431-44.

10.Yashin AI，Iachine IA，Begun AZ，etal.Hidden frailty：myths and reality.Research report，2001，34，48 pages.

（責(zé)任編輯：郭海強(qiáng)）

安徽省自然科學(xué)基金（090413126）

1.皖南醫(yī)學(xué)院預(yù)防醫(yī)學(xué)系（241002）

2.南通出入境檢驗(yàn)檢疫局

△通信作者：姚應(yīng)水，E-mail：yingshuiyao＠163.com