縱向有序數(shù)據(jù)的臨床療效評(píng)價(jià)方法應(yīng)用研究*

劉 艷 李 揚(yáng),2 劉 罡 張育銘

1.模型形式

1.研究背景

本文利用縱向有序logistic模型對(duì)肩周炎理療床臨床療效評(píng)價(jià)進(jìn)行了研究，結(jié)果表明該模型能夠有效地處理縱向有序數(shù)據(jù)，其固定效應(yīng)解釋了總體水平上變量之間的相互影響程度，隨機(jī)效應(yīng)解釋了數(shù)據(jù)間的相關(guān)、過度離散、異質(zhì)性等問題。在模型參數(shù)估計(jì)過程中，Gauss-Hermite積分和Quasi-Newton迭代算法克服了由二分類擴(kuò)展至多分類導(dǎo)致的參數(shù)增加、似然函數(shù)復(fù)雜、計(jì)算量大等困難。因此縱向有序logistic模型可以準(zhǔn)確地評(píng)價(jià)治療方案的有效性，反映影響因素之間的關(guān)系并體現(xiàn)出個(gè)體之間的差異性，為臨床療效評(píng)價(jià)提供了科學(xué)的依據(jù)。

縱向有序數(shù)據(jù)的臨床療效評(píng)價(jià)方法應(yīng)用研究*

劉 艷1李 揚(yáng)1,2劉 罡3張育銘1

目的 探討臨床療效研究中縱向有序數(shù)據(jù)的評(píng)價(jià)方法。方法 應(yīng)用廣義線性混合效應(yīng)模型，固定效應(yīng)解釋總體水平上變量之間的相互影響程度，隨機(jī)效應(yīng)解釋數(shù)據(jù)間的相關(guān)、過度離散、異質(zhì)性等問題。結(jié)果 在模型參數(shù)估計(jì)過程中，Gauss-Hermite積分和Quasi-Newton迭代算法克服了由二分類擴(kuò)展至多分類導(dǎo)致的參數(shù)增加、似然函數(shù)復(fù)雜、計(jì)算量大等困難。結(jié)論 縱向有序logistic模型可以準(zhǔn)確地評(píng)價(jià)治療方案的有效性，反映影響因素之間的關(guān)系并體現(xiàn)出個(gè)體之間的差異性，為臨床療效評(píng)價(jià)提供了科學(xué)的依據(jù)。

縱向有序數(shù)據(jù) logistic模型 臨床療效評(píng)價(jià)

臨床療效評(píng)價(jià)研究中存在著大量不獨(dú)立的縱向數(shù)據(jù)，例如：有遺傳、環(huán)境效應(yīng)的家系資料，縱向數(shù)據(jù)的個(gè)體資料，多中心的臨床試驗(yàn)數(shù)據(jù)以及具有嵌套效應(yīng)的毒理學(xué)實(shí)驗(yàn)數(shù)據(jù)等。這些數(shù)據(jù)來(lái)源于擁有相同親代、或是處于相同環(huán)境的個(gè)體，它們彼此相關(guān)，即使使用隨機(jī)盲法也不一定能夠?qū)⒂绊懸蛩赝耆刂贫鴥H考察研究因素的效果。因此，為了提高研究效率，研究人員需要建立一套有效的評(píng)價(jià)方法，以分析研究中的縱向數(shù)據(jù)。

在上述生物醫(yī)學(xué)實(shí)驗(yàn)中，評(píng)價(jià)實(shí)驗(yàn)結(jié)果、醫(yī)療成效的響應(yīng)變量通常是有序的。例如臨床上通常使用綜合反映疼痛和功能障礙[1]兩大主癥的指標(biāo)“有效性”評(píng)價(jià)肩周炎療效，將之分為四個(gè)等級(jí)，即治愈、顯效、有效(好轉(zhuǎn))和無(wú)效。可見，響應(yīng)變量“有效性”是有序變量。

對(duì)于這種數(shù)據(jù)，采取傳統(tǒng)的處理截面連續(xù)變量的分析方法是不可取的：其一，順序變量損失了89%～99%的信息[2]，降低了分析的有效性；其二，忽略了有序變量的天花板和地板效應(yīng)(ceiling and floor effects)，因此會(huì)產(chǎn)生相關(guān)的殘差和回歸變量，違反高斯馬爾科夫假定；其三，處理截面數(shù)據(jù)的方法無(wú)法反映變量之間時(shí)間上的相關(guān)性。另外，采取經(jīng)典的時(shí)間序列分析方法也是存在問題的：將ARMA模型與logistic模型進(jìn)行結(jié)合[3]，必須保證每個(gè)變量每次觀測(cè)的時(shí)間都是相同的，很多的生物、醫(yī)藥觀測(cè)數(shù)據(jù)都難以滿足此限制條件。由于生物醫(yī)藥研究的特殊性，受試個(gè)體可能因?yàn)楦鞣N原因不能在同一時(shí)刻獲得某一指標(biāo)，甚至可能退出研究項(xiàng)目。這樣產(chǎn)生的縱向數(shù)據(jù)有別于經(jīng)典時(shí)間序列分析的數(shù)據(jù)：每個(gè)對(duì)象不只有一個(gè)觀測(cè)值，這些觀測(cè)值是在不同的時(shí)間記錄的，每個(gè)對(duì)象的觀測(cè)次數(shù)不一定，觀測(cè)間隔不一定相等[4-5]。

目前，針對(duì)這種數(shù)據(jù)主要有兩種建模思路：邊際模型和廣義線性混合效應(yīng)模型。邊際模型主要應(yīng)用于對(duì)總體的平均水平進(jìn)行建模，常用GEE方法[6]進(jìn)行模型參數(shù)的估計(jì)。但邊際模型無(wú)法反映觀測(cè)個(gè)體之間的差異性。廣義線性混合效應(yīng)模型(GLMM)，通過引入隨機(jī)效應(yīng)[7]來(lái)刻畫數(shù)據(jù)間的相依性，體現(xiàn)觀測(cè)單元與總體間的偏差，允許觀測(cè)值的協(xié)方差結(jié)構(gòu)有異方差性，解決了過度離散(over dispersion)、異質(zhì)性(heterogeneity)等問題。在廣義線性混合效應(yīng)模型中，針對(duì)二分類響應(yīng)變量的處理方法已經(jīng)相對(duì)成熟，但處理有序響應(yīng)變量的方法還在不斷改進(jìn)。在二分類模型的基礎(chǔ)上，多分類模型需要引入新的參數(shù)——閾值(threshold)，使模型參數(shù)的估計(jì)變得更加復(fù)雜，估計(jì)方法也更多樣。

本文在臨床療效評(píng)價(jià)方法研究中引入縱向有序logistic模型，屬于廣義線性混合效應(yīng)模型在有序響應(yīng)變量上的一種應(yīng)用。本文所述方法克服了過去研究中對(duì)數(shù)據(jù)記錄時(shí)間要求苛刻、無(wú)法同時(shí)反映個(gè)體狀態(tài)隨時(shí)間變化以及個(gè)體間差異性等問題，并對(duì)以有序變量為療效指標(biāo)的治療方案進(jìn)行應(yīng)用嘗試。

縱向有序logistic模型

1.模型形式

(2)

2.模型參數(shù)的估計(jì)

(3)

(4)

(5)

為了最大化似然函數(shù)H，本文使用Quasi-Newton算法[11]來(lái)完成優(yōu)化。算法的迭代公式是x(k+1)=x(k)-λkGkgk，x(k)為第k次迭代后x的值，λk為步長(zhǎng)因子，gk為目標(biāo)函數(shù)在x(k)處的梯度，Gk滿足Gk+1Δgk=Δxk。設(shè)已知的目標(biāo)函數(shù)f(x)及梯度g(x)，問題的維度為r，并給出終止迭代的精度ε1，ε2和ε3，算法如表1所示。

表1 Quasi-Newton算法

將Gauss-Hermite積分和Quasi-Newton迭代算法結(jié)合使用可以極大化邊際似然函數(shù)，并得到相應(yīng)的取極大值的條件，從而得到固定效應(yīng)β、隨機(jī)效應(yīng)b的協(xié)方差陣∑b以及閾值γc的最大邊際似然估計(jì)。

3.模型參數(shù)的檢驗(yàn)

縱向有序logistic模型的檢驗(yàn)主要由三方面構(gòu)成：模型整體的似然比檢驗(yàn)(LRT)、模型每個(gè)參數(shù)的t檢驗(yàn)以及模型之間比較的BIC準(zhǔn)則[12]。

對(duì)于模型整體的檢驗(yàn)，原假設(shè)H0：所有模型參數(shù)均為0，H1：模型參數(shù)中至少有一個(gè)不為0。似然比檢驗(yàn)評(píng)估了原假設(shè)模型與備擇假設(shè)模型哪個(gè)更適合當(dāng)前數(shù)據(jù)分析。記H0、H1下模型的邊際似然函數(shù)分別為L(zhǎng)0、L1，LR=2(lnL1-lnL0)服從χ2分布，自由度為備擇假設(shè)相對(duì)原假設(shè)增加的參數(shù)的數(shù)目。若模型未通過似然比檢驗(yàn)，說明可能存在模型誤設(shè)。

對(duì)于具體參數(shù)的檢驗(yàn)，采用t檢驗(yàn)判斷其顯著性。檢驗(yàn)統(tǒng)計(jì)量t=b/Se(b)，服從自由度為n-2的t分布，n是數(shù)據(jù)數(shù)量。若自變量未通過t檢驗(yàn)，需要重新考察變量選擇的合理性。

對(duì)于模型之間的比較，一般選擇BIC準(zhǔn)則進(jìn)行模型優(yōu)劣的判別。貝葉斯信息量BIC=-2 ln(L)+ln(n)×k，其中L是在該模型下的最大似然，n是數(shù)據(jù)數(shù)量，k是模型的變量個(gè)數(shù)。BIC準(zhǔn)則綜合考慮了以上因素，防止模型的過度擬合，全面評(píng)價(jià)模型的擬合狀況，其值越小越好。

CGM理療床療效評(píng)價(jià)研究

1.研究背景

為評(píng)價(jià)“CGM理療床”治療肩周炎的有效性及安全性，研究人員將72名受試者隨機(jī)分為試驗(yàn)組(使用CGM理療床)和對(duì)照組，分別記錄每位患者的初期疾癥程度。經(jīng)過4周的臨床觀察，研究人員分別在1、2、4周記錄每位患者的生理指標(biāo)以及肩周炎有效性以研究治療方案的持續(xù)療效。試驗(yàn)研究的所有變量及其符號(hào)如表2所示。影響臨床療效的各項(xiàng)指標(biāo)(分組、周數(shù)、初期疾癥)都有3個(gè)時(shí)期的觀測(cè)，這些觀測(cè)值之間是相關(guān)的。評(píng)價(jià)療效的指標(biāo)(有效性)是一個(gè)順序變量，受到其他指標(biāo)與患者個(gè)人因素的影響。個(gè)人因素是無(wú)法觀測(cè)的潛變量，假定它服從均值為0 的多元正態(tài)分布。在這些條件下，建立一個(gè)縱向有序logistic模型可以有效地解決療效評(píng)價(jià)的問題。

表2 試驗(yàn)研究變量表

2.模型的建立與解釋

為確定縱向有序logistic模型的基本形式，本文通過圖像初步探究各個(gè)變量的顯著性以及自變量間可能存在的交互效應(yīng)。圖1、圖2(對(duì)照組與試驗(yàn)組患者研究時(shí)點(diǎn)有效性分布的柱狀圖)表明使用理療床的患者四周后的康復(fù)程度顯著優(yōu)于未接受該治療方案的患者。圖3是所有受試者四周的有效性變化，圖4是對(duì)照組與試驗(yàn)組四周平均有效性變化折線圖。圖4中兩條折線近乎平行，表明時(shí)間與治療方案之間沒有交互效應(yīng)。由于缺乏治療方案療效與個(gè)人因素之間關(guān)系的信息，治療方案是否存在倍差估計(jì)量無(wú)法確定，因此本文設(shè)定了以下兩個(gè)可能的縱向有序logistic模型。隨機(jī)系數(shù)模型認(rèn)為治療方案會(huì)對(duì)不同的人產(chǎn)生不同的療效，這種差異包含了兩層含義：第一，在對(duì)照組和試驗(yàn)組，不同個(gè)體的療效都會(huì)受個(gè)體固有差異的影響。第二，在試驗(yàn)組，除了個(gè)體固有差異，療效還額外受治療方案附帶的個(gè)人因素影響。隨機(jī)截距模型則認(rèn)為無(wú)論是對(duì)照組還是試驗(yàn)組，療效都只受個(gè)體固有差異影響。

圖1 對(duì)照組有效性隨時(shí)間變化

模型1：隨機(jī)系數(shù)模型

λic=γc-(β1group+β2week+β3severity+b0i+b1igroup)

模型2：隨機(jī)截距模型

λic=γc-(β1group+β2week+β3severity+b0i)

模型參數(shù)可用SAS的NLMIXED過程進(jìn)行估計(jì)(SAS代碼見附錄)，兩個(gè)模型分別經(jīng)過20和17步迭代達(dá)到收斂條件，得到表3的參數(shù)估計(jì)結(jié)果。

圖2 試照組有效性隨時(shí)間變化

圖3 時(shí)間折線圖

圖4 交互效應(yīng)檢驗(yàn)

表3 參數(shù)估計(jì)

結(jié)論及展望

本文利用縱向有序logistic模型對(duì)肩周炎理療床臨床療效評(píng)價(jià)進(jìn)行了研究，結(jié)果表明該模型能夠有效地處理縱向有序數(shù)據(jù)，其固定效應(yīng)解釋了總體水平上變量之間的相互影響程度，隨機(jī)效應(yīng)解釋了數(shù)據(jù)間的相關(guān)、過度離散、異質(zhì)性等問題。在模型參數(shù)估計(jì)過程中，Gauss-Hermite積分和Quasi-Newton迭代算法克服了由二分類擴(kuò)展至多分類導(dǎo)致的參數(shù)增加、似然函數(shù)復(fù)雜、計(jì)算量大等困難。因此縱向有序logistic模型可以準(zhǔn)確地評(píng)價(jià)治療方案的有效性，反映影響因素之間的關(guān)系并體現(xiàn)出個(gè)體之間的差異性，為臨床療效評(píng)價(jià)提供了科學(xué)的依據(jù)。

然而，實(shí)證研究中的數(shù)據(jù)分布不一定符合模型假設(shè)的均勻分布，使用logit變換作為連接函數(shù)略有不妥。隨著縱向數(shù)據(jù)的廣泛應(yīng)用，對(duì)它的研究也將更深入，針對(duì)一般的縱向數(shù)據(jù)要用更復(fù)雜的模型，更深厚廣泛的知識(shí)去解決。在連接函數(shù)問題上，就有很多問題值得研究：在高層類別出現(xiàn)幾率較大的情況下應(yīng)使用Complementary log-log作為連接函數(shù)；在底層類別出現(xiàn)幾率較大的情況下應(yīng)使用Negative log-log作為連接函數(shù)；在兩端類別出現(xiàn)較大的情況下應(yīng)使用Cauchit變換等。另一方面，雖然隨機(jī)效應(yīng)模型可以反映個(gè)體差異的變化，但其參數(shù)結(jié)果不利于解釋針對(duì)研究對(duì)象總體的整體療效。因此，在后續(xù)研究中可以考慮引入邊際化隨機(jī)效應(yīng)模型的研究思路，通過隨機(jī)效應(yīng)部分的選擇增強(qiáng)對(duì)數(shù)據(jù)相關(guān)性和異方差性的利用，同時(shí)通過邊際模型部分增強(qiáng)變量選擇結(jié)果的解釋性，提升理論方法的實(shí)證價(jià)值。

附錄：

/* 利用NLMIXED過程估計(jì)隨機(jī)系數(shù)模型參數(shù)*/

PROC NLMIXED data=one QPOINTS=21;

PARMS b1=-1.38 b2=-0.719 b3=-0.36 v0=1 v1=1 v01=0 g1=-9.9 g2=-6.3 g3=-1.9;/*設(shè)定參數(shù)初始值*/

z=b1*group+b2*week+b3*severity+u0+u1*group;

IF(pain=1) THEN

p=CDF(′NORMAL′,g1-z);

ELSE IF(pain=2) THEN

p=CDF(′NORMAL′,g2-z)-CDF(′NORMAL′,g1-z);

ELSE IF(pain=3) THEN

p=CDF(′NORMAL′,g3-z)-CDF(′NORMAL′,g2-z);

ELSE IF(pain=4) THEN

p=1-CDF(′NORMAL′,g3-z);

loglik=LOG(p);

MODEL pain ～ GENERAL(loglik);

RANDOM u0 u1 ～ NORMAL([0,0],[v0,v01,v1]) SUBJECT=id out=ebest2b;

ESTIMATE′recorr′v01/SQRT(v0*v1);

RUN;

/* 利用NLMIXED過程估計(jì)隨機(jī)截距模型參數(shù) */

PROC NLMIXED data=one QPOINTS=21;

PARMS b1=-1.38 b2=-0.719 b3=-0.36 v0=1 g1=-9.9 g2=-6.3 g3=-1.9;/*設(shè)定參數(shù)初始值*/

z=b1*group+b2*week+b3*severity+u0 ;

IF(pain=1) THEN

p=CDF(′NORMAL′,g1-z);ELSE IF(pain=2) THEN

p=CDF(′NORMAL′,g2-z)-CDF(′NORMAL′,g1-z);

ELSE IF(pain=3) THEN

p=CDF(′NORMAL′,g3-z)-CDF(′NORMAL′,g2-z);

ELSE IF(pain=4) THEN

p=1-CDF(′NORMAL′,g3-z);

loglik=LOG(p);

MODEL pain ～ GENERAL(loglik);

RANDOM u0 ～ NORMAL(0,v0) SUBJECT=id out=ebest2b;

RUN;

注意：執(zhí)行以上程序前已將數(shù)據(jù)存放在包含week、group、severity、pain變量的數(shù)據(jù)集one中。

[1]胡幼平,刁驤,楊運(yùn)寬,等.肩周炎臨床療效評(píng)定方法概況.江西中醫(yī)藥,2007，38(9):63-66.

[2]Armstrong BG,Sloan M.Ordinal regression models for epidemiologic data.American Journal of Epidemiology,1989,129(1):191-204.

[3]Wooldridge J.Introductory Econometrics.Fourth Edition.北京：中國(guó)人民大學(xué)出版社.

[4]Varin C,Czado C.A mixed autoregressive probit model for ordinal longitudinal data.Biostatistics,2010,11(1):127-138.

[5]Hartzel J,Agresti A,Caffo B.Multinomial logit random effects models.Statistical Modelling 2001,1(2):81-102

[6]Liang KY,Zeger SL.Longitudinal data analysis using generalized linear models.Biometrika,1986,73(1):12-22.

[7]Hedeker D.A mixed-effects multinomial logistic regression model.Statistics in medicine,2003,22(9):1433-1446.

[8]Liu LC,Hedeker D.A Mixed-Effects Regression Model for Longitudinal Multivariate Ordinal Data.Biometrics,2006，62(1):261-268.

[9]Gardiner JC,Luo Z,Roman LA.Fixed effects,random effects and GEE:What are the differences.Statistics in medicine,2009,28(2):221-239.

[10]Hedeker D,Gibbons RD.A Random-Effects Ordinal Regression Model for Multilevel Analysis.Biometrics,1994，50(4):933-944.

[11]高惠璇.統(tǒng)計(jì)計(jì)算.北京：北京大學(xué)出版社，2012：375-379.

[12]吳喜之.復(fù)雜數(shù)據(jù)統(tǒng)計(jì)方法——基于R的應(yīng)用.北京：中國(guó)人民大學(xué)出版社，2012：92-98.

(責(zé)任編輯：劉 壯)

本文是中國(guó)人民大學(xué)科學(xué)研究基金(中央高校基本科研業(yè)務(wù)費(fèi)專項(xiàng)資金資助)項(xiàng)目(15XNI011)

1.中國(guó)人民大學(xué)統(tǒng)計(jì)學(xué)院(100872)

2.中國(guó)人民大學(xué)應(yīng)用統(tǒng)計(jì)科學(xué)研究中心

3.哈佛大學(xué)生物統(tǒng)計(jì)系