纵向有序数据的临床疗效评价方法应用研究*

刘 艳 李 扬,2 刘 罡 张育铭

1.模型形式

1.研究背景

本文利用纵向有序logistic模型对肩周炎理疗床临床疗效评价进行了研究，结果表明该模型能够有效地处理纵向有序数据，其固定效应解释了总体水平上变量之间的相互影响程度，随机效应解释了数据间的相关、过度离散、异质性等问题。在模型参数估计过程中，Gauss-Hermite积分和Quasi-Newton迭代算法克服了由二分类扩展至多分类导致的参数增加、似然函数复杂、计算量大等困难。因此纵向有序logistic模型可以准确地评价治疗方案的有效性，反映影响因素之间的关系并体现出个体之间的差异性，为临床疗效评价提供了科学的依据。

纵向有序数据的临床疗效评价方法应用研究*

刘 艳1李 扬1,2刘 罡3张育铭1

目的 探讨临床疗效研究中纵向有序数据的评价方法。方法 应用广义线性混合效应模型，固定效应解释总体水平上变量之间的相互影响程度，随机效应解释数据间的相关、过度离散、异质性等问题。结果 在模型参数估计过程中，Gauss-Hermite积分和Quasi-Newton迭代算法克服了由二分类扩展至多分类导致的参数增加、似然函数复杂、计算量大等困难。结论 纵向有序logistic模型可以准确地评价治疗方案的有效性，反映影响因素之间的关系并体现出个体之间的差异性，为临床疗效评价提供了科学的依据。

纵向有序数据 logistic模型 临床疗效评价

临床疗效评价研究中存在着大量不独立的纵向数据，例如：有遗传、环境效应的家系资料，纵向数据的个体资料，多中心的临床试验数据以及具有嵌套效应的毒理学实验数据等。这些数据来源于拥有相同亲代、或是处于相同环境的个体，它们彼此相关，即使使用随机盲法也不一定能够将影响因素完全控制而仅考察研究因素的效果。因此，为了提高研究效率，研究人员需要建立一套有效的评价方法，以分析研究中的纵向数据。

在上述生物医学实验中，评价实验结果、医疗成效的响应变量通常是有序的。例如临床上通常使用综合反映疼痛和功能障碍[1]两大主症的指标“有效性”评价肩周炎疗效，将之分为四个等级，即治愈、显效、有效(好转)和无效。可见，响应变量“有效性”是有序变量。

对于这种数据，采取传统的处理截面连续变量的分析方法是不可取的：其一，顺序变量损失了89%～99%的信息[2]，降低了分析的有效性；其二，忽略了有序变量的天花板和地板效应(ceiling and floor effects)，因此会产生相关的残差和回归变量，违反高斯马尔科夫假定；其三，处理截面数据的方法无法反映变量之间时间上的相关性。另外，采取经典的时间序列分析方法也是存在问题的：将ARMA模型与logistic模型进行结合[3]，必须保证每个变量每次观测的时间都是相同的，很多的生物、医药观测数据都难以满足此限制条件。由于生物医药研究的特殊性，受试个体可能因为各种原因不能在同一时刻获得某一指标，甚至可能退出研究项目。这样产生的纵向数据有别于经典时间序列分析的数据：每个对象不只有一个观测值，这些观测值是在不同的时间记录的，每个对象的观测次数不一定，观测间隔不一定相等[4-5]。

目前，针对这种数据主要有两种建模思路：边际模型和广义线性混合效应模型。边际模型主要应用于对总体的平均水平进行建模，常用GEE方法[6]进行模型参数的估计。但边际模型无法反映观测个体之间的差异性。广义线性混合效应模型(GLMM)，通过引入随机效应[7]来刻画数据间的相依性，体现观测单元与总体间的偏差，允许观测值的协方差结构有异方差性，解决了过度离散(over dispersion)、异质性(heterogeneity)等问题。在广义线性混合效应模型中，针对二分类响应变量的处理方法已经相对成熟，但处理有序响应变量的方法还在不断改进。在二分类模型的基础上，多分类模型需要引入新的参数——阈值(threshold)，使模型参数的估计变得更加复杂，估计方法也更多样。

本文在临床疗效评价方法研究中引入纵向有序logistic模型，属于广义线性混合效应模型在有序响应变量上的一种应用。本文所述方法克服了过去研究中对数据记录时间要求苛刻、无法同时反映个体状态随时间变化以及个体间差异性等问题，并对以有序变量为疗效指标的治疗方案进行应用尝试。

纵向有序logistic模型

1.模型形式

(2)

2.模型参数的估计

(3)

(4)

(5)

为了最大化似然函数H，本文使用Quasi-Newton算法[11]来完成优化。算法的迭代公式是x(k+1)=x(k)-λkGkgk，x(k)为第k次迭代后x的值，λk为步长因子，gk为目标函数在x(k)处的梯度，Gk满足Gk+1Δgk=Δxk。设已知的目标函数f(x)及梯度g(x)，问题的维度为r，并给出终止迭代的精度ε1，ε2和ε3，算法如表1所示。

表1 Quasi-Newton算法

将Gauss-Hermite积分和Quasi-Newton迭代算法结合使用可以极大化边际似然函数，并得到相应的取极大值的条件，从而得到固定效应β、随机效应b的协方差阵∑b以及阈值γc的最大边际似然估计。

3.模型参数的检验

纵向有序logistic模型的检验主要由三方面构成：模型整体的似然比检验(LRT)、模型每个参数的t检验以及模型之间比较的BIC准则[12]。

对于模型整体的检验，原假设H0：所有模型参数均为0，H1：模型参数中至少有一个不为0。似然比检验评估了原假设模型与备择假设模型哪个更适合当前数据分析。记H0、H1下模型的边际似然函数分别为L0、L1，LR=2(lnL1-lnL0)服从χ2分布，自由度为备择假设相对原假设增加的参数的数目。若模型未通过似然比检验，说明可能存在模型误设。

对于具体参数的检验，采用t检验判断其显著性。检验统计量t=b/Se(b)，服从自由度为n-2的t分布，n是数据数量。若自变量未通过t检验，需要重新考察变量选择的合理性。

对于模型之间的比较，一般选择BIC准则进行模型优劣的判别。贝叶斯信息量BIC=-2 ln(L)+ln(n)×k，其中L是在该模型下的最大似然，n是数据数量，k是模型的变量个数。BIC准则综合考虑了以上因素，防止模型的过度拟合，全面评价模型的拟合状况，其值越小越好。

CGM理疗床疗效评价研究

1.研究背景

为评价“CGM理疗床”治疗肩周炎的有效性及安全性，研究人员将72名受试者随机分为试验组(使用CGM理疗床)和对照组，分别记录每位患者的初期疾症程度。经过4周的临床观察，研究人员分别在1、2、4周记录每位患者的生理指标以及肩周炎有效性以研究治疗方案的持续疗效。试验研究的所有变量及其符号如表2所示。影响临床疗效的各项指标(分组、周数、初期疾症)都有3个时期的观测，这些观测值之间是相关的。评价疗效的指标(有效性)是一个顺序变量，受到其他指标与患者个人因素的影响。个人因素是无法观测的潜变量，假定它服从均值为0 的多元正态分布。在这些条件下，建立一个纵向有序logistic模型可以有效地解决疗效评价的问题。

表2 试验研究变量表

2.模型的建立与解释

为确定纵向有序logistic模型的基本形式，本文通过图像初步探究各个变量的显著性以及自变量间可能存在的交互效应。图1、图2(对照组与试验组患者研究时点有效性分布的柱状图)表明使用理疗床的患者四周后的康复程度显著优于未接受该治疗方案的患者。图3是所有受试者四周的有效性变化，图4是对照组与试验组四周平均有效性变化折线图。图4中两条折线近乎平行，表明时间与治疗方案之间没有交互效应。由于缺乏治疗方案疗效与个人因素之间关系的信息，治疗方案是否存在倍差估计量无法确定，因此本文设定了以下两个可能的纵向有序logistic模型。随机系数模型认为治疗方案会对不同的人产生不同的疗效，这种差异包含了两层含义：第一，在对照组和试验组，不同个体的疗效都会受个体固有差异的影响。第二，在试验组，除了个体固有差异，疗效还额外受治疗方案附带的个人因素影响。随机截距模型则认为无论是对照组还是试验组，疗效都只受个体固有差异影响。

图1 对照组有效性随时间变化

模型1：随机系数模型

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

模型2：随机截距模型

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

模型参数可用SAS的NLMIXED过程进行估计(SAS代码见附录)，两个模型分别经过20和17步迭代达到收敛条件，得到表3的参数估计结果。

图2 试照组有效性随时间变化

图3 时间折线图

图4 交互效应检验

表3 参数估计

结论及展望

本文利用纵向有序logistic模型对肩周炎理疗床临床疗效评价进行了研究，结果表明该模型能够有效地处理纵向有序数据，其固定效应解释了总体水平上变量之间的相互影响程度，随机效应解释了数据间的相关、过度离散、异质性等问题。在模型参数估计过程中，Gauss-Hermite积分和Quasi-Newton迭代算法克服了由二分类扩展至多分类导致的参数增加、似然函数复杂、计算量大等困难。因此纵向有序logistic模型可以准确地评价治疗方案的有效性，反映影响因素之间的关系并体现出个体之间的差异性，为临床疗效评价提供了科学的依据。

然而，实证研究中的数据分布不一定符合模型假设的均匀分布，使用logit变换作为连接函数略有不妥。随着纵向数据的广泛应用，对它的研究也将更深入，针对一般的纵向数据要用更复杂的模型，更深厚广泛的知识去解决。在连接函数问题上，就有很多问题值得研究：在高层类别出现几率较大的情况下应使用Complementary log-log作为连接函数；在底层类别出现几率较大的情况下应使用Negative log-log作为连接函数；在两端类别出现较大的情况下应使用Cauchit变换等。另一方面，虽然随机效应模型可以反映个体差异的变化，但其参数结果不利于解释针对研究对象总体的整体疗效。因此，在后续研究中可以考虑引入边际化随机效应模型的研究思路，通过随机效应部分的选择增强对数据相关性和异方差性的利用，同时通过边际模型部分增强变量选择结果的解释性，提升理论方法的实证价值。

附录：

/* 利用NLMIXED过程估计随机系数模型参数*/

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;/*设定参数初始值*/

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过程估计随机截距模型参数 */

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;/*设定参数初始值*/

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;

注意：执行以上程序前已将数据存放在包含week、group、severity、pain变量的数据集one中。

[1]胡幼平,刁骧,杨运宽,等.肩周炎临床疗效评定方法概况.江西中医药,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.北京：中国人民大学出版社.

[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]高惠璇.统计计算.北京：北京大学出版社，2012：375-379.

[12]吴喜之.复杂数据统计方法——基于R的应用.北京：中国人民大学出版社，2012：92-98.

(责任编辑：刘 壮)

本文是中国人民大学科学研究基金(中央高校基本科研业务费专项资金资助)项目(15XNI011)

1.中国人民大学统计学院(100872)

2.中国人民大学应用统计科学研究中心

3.哈佛大学生物统计系