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

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

1.模型形式

1.研究背景

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

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

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

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

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

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

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

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

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

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

縱向有序logistic模型

1.模型形式

(2)

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

(3)

(4)

(5)

為了最大化似然函數(shù)H，本文使用Quasi-Newton算法[11]來完成優(yōu)化。算法的迭代公式是x(k+1)=x(k)-λkGkgk，x(k)為第k次迭代后x的值，λk為步長因子，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的最大邊際似然估計。

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

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

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

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

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

CGM理療床療效評價研究

1.研究背景

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

表2 試驗研究變量表

2.模型的建立與解釋

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

圖1 對照組有效性隨時間變化

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

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

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

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

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

圖2 試照組有效性隨時間變化

圖3 時間折線圖

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

表3 參數(shù)估計

結(jié)論及展望

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

然而，實證研究中的數(shù)據(jù)分布不一定符合模型假設(shè)的均勻分布，使用logit變換作為連接函數(shù)略有不妥。隨著縱向數(shù)據(jù)的廣泛應(yīng)用，對它的研究也將更深入，針對一般的縱向數(shù)據(jù)要用更復(fù)雜的模型，更深厚廣泛的知識去解決。在連接函數(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)模型可以反映個體差異的變化，但其參數(shù)結(jié)果不利于解釋針對研究對象總體的整體療效。因此，在后續(xù)研究中可以考慮引入邊際化隨機(jī)效應(yīng)模型的研究思路，通過隨機(jī)效應(yīng)部分的選擇增強(qiáng)對數(shù)據(jù)相關(guān)性和異方差性的利用，同時通過邊際模型部分增強(qiáng)變量選擇結(jié)果的解釋性，提升理論方法的實證價值。

附錄：

/* 利用NLMIXED過程估計隨機(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ī)截距模型參數(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)寬,等.肩周炎臨床療效評定方法概況.江西中醫(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.北京：中國人民大學(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)計計算.北京：北京大學(xué)出版社，2012：375-379.

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

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

本文是中國人民大學(xué)科學(xué)研究基金(中央高?；究蒲袠I(yè)務(wù)費(fèi)專項資金資助)項目(15XNI011)

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

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

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