２０１１年中国卫生统计学年会论文集
多变量随机效应模型及其ＳＡＳ实现
刘文华，吴家利，杨杨，宋婷婷，章顺悦，陈远方，尹平△
华中科技犬学尉济医学院公共卫生学院流行病与卫生统计学系（４３００３０）
摘要目的介绍多变量随机效应模型及其参数估计，促进其在多个阈值诊断试验ｍｅｔａ分析中的应用。方法对双变量随机效应模型进行扩展和延伸，构建多变量随机效应模型，采月｛｝最大似然法，并利用ＳＡＳ中的ＰＲＯＣＮＬＭＩＸＥＤ过程进行统计实现。结果采用多变量随机效应模型可得到每一闽值的合并灵敏度、特异度、诊断比值比，及曲线下面积ＡＵＣ等评价指标。亦可绘制ＳＲＯＣ曲线及其置信带、预测带等。结论多变量随机效应模型能够充分利用各研究的原始数据信息，得到稳定可靠的结果，有着较强的应用价值。
关键词多变量随机效应模型闽值诊断试验ｍｅｔａ分析
Ｍｕｌｔｉｖａｒｉａｔｅｒａｎｄｏｍ
ＬｉｕｅｆｆｅｃｔｓｍｏｄｅｌａｎｄｉｔｓｓＡＳｐｒｏｇｒａｍｓＷｅｎｈｕａ，ＷｕＪｉａｌｉ，ＹａｎｇＹａｎｇ，ＳｏｎｇＴｉｎｇｔｉｎｇ，ＺｈａｎｇＳｈｕｎｙｕｅ，ＣｈｅｎＹｕａｎｇｆａｎｇ，ＹｉｎＰｉｎｇ
ＤｅｐａｒｔｍｅｎｔｏｆＥｐｉｄｅｍｉｏｌｏｇｙａｎｄＨｅａｌｔｈＳｔａｔｉｓｔｉｃｓ，ＳｃｈｏｏｌｏｆＰｕｂｌｉｃ
ｏｆＳｃｉｅｎｃｅａｎｄＨｅａｔｈ，ＴｏｎｇｉｉＭｅｄｉｃａｌＣｏｌｌｅｇｅｏｆＨｕａｚｈｏｎｇＵｎｉｖｅｒｓｉｔｙＴｅｃｈｎｏｌｏｇｙ｛４３００３０），Ｗｕｈａｎ
ｔｏＡｂｓｔｒａｃｔ０Ｍｅｃｔｉｖｅ：Ｔｈｅｓｔｕｄｙｗａｓｃｏｎｄｕｃｔｅｄｉｎｔｒｏｄｕｃｅｔｈｅｍｕｌｔｉｖａｒｉａｔｅｒａｎｄｏｍｅｆｆｅｃｔｓｍｏｄｅｌａｎｄｉｔｓ
ｐａｒａｍｅｔｅｒｅｓｔｉｍａｔｉｏｎ，ａｎｄｔｏｆａｃｉｌｉｔｙｉｔｓａｐｐｌｉｃａｔｉｏｎｓｉｎｍｅｔａ－ａｎａｌｙｓｉｓｏｆｄｉａｇｎｏｓｔｉｃｔｅｓｔｓｗｉｔｈｍｕｌｔｉｐｌｅｔｈｒｅｓｈｏｌｄｓ．Ｍｅｔｈｏｄｓ：Ｔｈｉｓｐａｐｅｒｅｘｐａｎｄｅｄａｎｄｅｘｔｅｎｄｅｄｔｈｅｂｉｖａｒｉａｔｅｒａｎｄｏｍｅｆｆｅｃｔｓｍｏｄｅｌｔｏｄｅｖｅｌｏｐｔｈｅｍｕｌｔｉｖａｒｉａｔｅｒａｎｄｏｍｅｆｆｅｃｔｓｍｏｄｅｌ，ａｎｄｉｍｐｌｅｍｅｎｔｅｄｉｔｉｎｔｈｅＳＡＳＰｒｏｃＮＬＭＩＸＥＤｐｒｏｃｅｄｕｒｅｗｉｔｈｍａｘｉｍｕｍｌｉｋｅｌｉｈｏｏｄ．Ｒｅｓｕｌｔｓ：Ｔｈｅｍｕｌｔｉｖａｒｉａｔｅｒａｎｄｏｍｅｆｆｅｃｔｓｍｏｄｅｌａｌｌｏｗｓｃａｌｃｕｌａｔｉｎｇｔｈｅｐｏｏｌｅｄ
ａｒｅａｓｅｎｓｉｔｉｖｉｔｙ，ｓｐｅｃｉｆｉｃｉｔｙ，ｄｉａｇｎｏｓｔｉｃｏｄｄｓｒａｔｉｏｆｏｒｅａｃｈｔｈｒｅｓｈｏｌｄａｎｄｔｈｅ
ｔｈｅＳＲＯＣｃｕｒｖｅｕｎｄｅｒｔｈｅｃｕｒｖｅ，ｗｈｉｌｅｄｒａｗｉｎｇｗｉｔｈｉｔｓｃｏｎｆｉｄｅｎｃｅｂａｎｄａｎｄｐｒｅｄｉｃｔｉｏｎｂａｎｄ．Ｃｏｎｃｌｕｓｉｏｎｓ：Ｔｈｅｍｕｌｔｉｖａｒｉａｔｅｒａｎｄｏｍ
ｆｕｌｌｙｅｘｐｌｏｉｔｔｈｅｄａｔａ，ｏｂｔａｉｎｓｔａｂｌｅａｎｄｒｅｌｉａｂｌｅｒｅｓｕｌｔｓ，ａｎｄｈａｖｅａｅｆｆｅｃｔｓｍｏｄｅｌｃａｎｇｏｏｄａｐｐｌｉｃａｔｉｏｎ．
Ｋｅｙｗｏｒｄｓ：ｍｕｌｔｉｖａｒｉａｔｅｒａｎｄｏｍｅｆｆｅｃｔｓｍｏｄｅｌ；ｔｈｒｅｓｈｏｌｄ；ｄｉａｇｎｏｓｔｉｃｔｅｓｔｓ；ｍｅｔａ－ａｎａｌｙｓｉｓ
近年来，诊断试验ｍｅｔａ分析在医学研究中应用越来越广泛，其统计学方法主要针对每个研究只有一个诊断阈值，提供一对灵敏度和特异度或一个四格表的资料，一般采用分别合并灵敏度、特异度，或采用固定或随机效应模型联合灵敏度和特异度进行分析【１，２１。但在实际工作中会经常遇到检测结果有多个诊断阈值的情形，针对这类问题，以往ｍｅｔａ分析常用的方法是将其转换为二分类，即每个研究根据一个闽值得到一对灵敏度和特异度，再应用ＳＲＯＣ方法或双变量随机效应模型进行分析。而这种“降维”的方法会损失数据信息，降低数据的利用度。解决这类多阂值问题的有效方法是多变量随机效应模型，其可以充分利用原始数据信息，获得更为可靠的结论。国际上已进行了一些研究
【３Ｊ’但国内的相关研究报道甚少。为此，本文对多变量随机效应模型进行介绍，以期能更好地促进国内该方法在多个阂值诊断试验ｍｅｔａ分析中的应用。△通讯作者：尹平
２０１１年中国卫生统计学年会论文集
一、模型与参数估计
多变量随机效应模型可看作双变量随机效应模型的扩展和延伸．在双变量随机效应模型中，灵敏度和１．特异度经Ｌｏｇｉｔ转换后的真实值服从双变量正态分布，灵敏度和１一特异度的真实值分别服从二项分布。而在多变量随机效应模型中，针对多个阈值，可得到相应的多对灵敏度和１一特异度，经Ｌｏｇｉｔ转换后其真实值服从多变量正态分布，而检测试验取每个结果的频数在疾病人群和非疾病人群中均服从多项分布。根据乘积多项分布理论，可获得其联合概率，从而构造出似然函数，并采用最大似然法进行参数估计。
假设经纳入标准和排除标准筛选，１个研究，每个研究有ｇ－１个阈值，则检测指标可分为ｇ个组段或ｇ个有序变量。对于研究ｆ，按检测结果从高到低，依次列出各组段或等级的观测数，得到一个２ｘｇ列联表，如表１所示。
表１第ｆ研究金标准与检测试验结果的２×ｇ频数分布
列联表中每个格子的观测数目为随机变量，在病例人群中服从参数为ｎｌｉ，刀，∥…，兀岵的多项分布，在非病例人群中服从参数为，＇ｌＯｉ，巧０ｆｌ，…，石ｏｆｇ的多项分布。
以测定值大于阈值判断为检测试验阳性，则对于研究；，阈值，，其真阳性率研碾；，；翌，假。
ｎｌｉ
眦率脚。；擎，将珊玎，肿。进弛ｇ；ｔ转换，参数ｙ俩郴ｕ表示撇盯，脚。转换后的总体真实值，则殉，却，问线性关系为：
ＹⅡ＝口ｆ＋肛目
而ｘｄ，可用如下模型表示：
ｘ盯＝ｘ，＋６ｆ＋Ｐ口（１）（２）
式中，口；表示各研究常数项参数；斜率卢为对称性参数，各研究斜率固定，卢＝１表示ＳＲＯＣ曲线围绕“灵敏度：特异度，，直线对称；；，表示对于诊断阈值．『，‰的平均值；６ｉ表示系统误差，：白表示残差。该模型还作如下假设：Ｏｔｉ。Ⅳ五，仃。２），６，一Ｎ（Ｏ，仃ｄ２），ｅ．｝：－Ｎ（Ｏ，仃。２），仃。＝０，仃＆＝０ａ因此模型也可写成如下矩阵形式：
——————————————————————————————————————————————，——————————————————一一．一！！！！兰！望墨竺篁盐堂兰全笙壅篓一——
…
ａｄｒ口ｄｄｒ甜仃神
ｉ
２２…
ｘｉｌｏ口＋０ｅ仃６２仃６２
：
●
～Ｎｆ３）
Ｘ每仃６‘．Ｄｄ
●
：骥二
Ｚｉ，ｇ一１ｚｇ—ｌ』Ｉ口面…２
仃６口６…口ｄ２＋仃Ｐ２
采用最人似然法进行参数估计，根据乘积多项分布原理，其联合概率密度函数为：厂（，ｚ。玎，，ｚ・玎ｋ。ｒ，九，ｚ，万叫，万ｔ玎）２ｌ＝ｌｎ‘ｉ丽ｎｏｉ！．ｎｌｉ！导ｇｊ万叫疗叫万－玎”ｌ口）ｃ４，
ｅＸＯ，＿『＝１
１４－ｅＸ，Ｊ
其中，石叫＝‰（月孵）＝三一』兰，＿『．２，…，ｇ－１，。（５）。
１４－ｅｂ１＋ｅｘｉ＂卜１
１一』三，＿『；ｇ。。。
１＋ｅＸ，，Ｊ一１
ｇｙｇ
ａｊ；１
１＋ｅ’口
刀埘＝置牙Ｏｌ帮）＝』二一』二，－『：２，…，ｇ一１，。。。～（６）
１＋ｅｙ９１．＋ｅｙｉ．卜１
１一』兰，’『；ｇ。。。
１＋ｅｙｅ，卜１
骞鼽ｎｇ
Ｉ，ＪＪｃｌ－１，善荟～“。Ｊ—ＪＪ。Ｉ
因此，似然函数为：Ｌ＝Ｆ（臼卜ｏｉｊ＇ｎｌＯ）＝，Ｑ嘶，ｎ坷ｎｏｉ，ｎｉｌ∥Ｏｑ＇万坷）。同时，为得到非偏估计。采用迭代法不断修正最大似然估计值。
当然，也可将斜率参数芦设为随机变量，同对，为校正研究间协变量（ｚ）的影响，模型可写成：ｙ可＝绞ｆ＋屈ｘ口＋弦ａ进行参数估计时需增加仃芦２、仃肚２、仃彤２、），等４个参数，但较多的参数在一定程度上会影响模型的收敛性。
利用ＳＡＳ的ＰＲＯＣＮＬＭＩＸＥＤ过程可以完成多变量随机效应模型的拟合，得到每一阂值的合并灵敏度、特异度、诊断比值比，曲线下面积ＡＵＣ等评价指标。还可绘制ＳＲＯＣ曲线及其置信带、预测带，以及每个研究的ＲＯＣ曲线。１８９
二、ＳＡＳ实现
２００６年Ｋｏｂａｙａｓｈｉ等【４】提出了最新的川崎病并发冠脉损伤的高危评分体系，总积分为１１分，Ｏ～３分为低危组，４分以上为高危人群，７分以上为极高危人群。本文模拟Ｋｏｂａｙａｓｈｉ评分系统检测川崎病患者是否并发冠脉损伤的数据，将评分结果分为＝３分，４－６分，＝７分三组，按照川崎病和冠脉损伤诊断标准将川崎病ＩＶＩＧ无反应者分为ＣＡＬ（冠脉损伤）组和ＮＣＡＬ（无冠脉损伤）组，评价该评分体系的诊断准确性。ＳＡＳ过程如下。
ｄａｔａｄａｔａ；
ｌｅｎｇｔｈａｕｔｈｅｒａｒｅａ￥１０：
ｉｎｐｕｔｓｔｕｄｙａｕｔｈｅｒ￥ｙｅａｒａｒｅａ￥ｈｅａ＿ｈｉｈｅａ—ｍｉｄｈｅａ一１０Ｗｄｉｓ—ｈｉｄｉｓ—ｍｉｄｄｉｓ＿ｌｏｗ；
ｈｅａ＿＿ｎ＝ｈｅａ＿ｈｉ＋ｈｅａｍｉｄ＋ｈｅａｌｏｗ；
ｄｉｓ＿ｎ＝ｄｉｓ＿ｈｉ＋ｄｉｓ＿ｍｉｄ＋ｄｉｓ＿ｌｏｗ；
ｄａｔａｌｉｎｅｓ；
１Ｋｏｂａｙａｓｈｉ２００６ＵＳ６３９４８２１００１１６７８
２Ｙａｓｕｎｏｒｉ２００７ＵＳ９５９６０３２９３８８１
３Ｉｎｏｕｅ２００８ＵＳ１０８５７００３０５０９５
４Ｔａｋｅｕｃｈｉ２００８ＵＳ２８１５６２９９１８３０４
５Ａｅｒｔｇｅｅｒｔｓ２００９ＵＳ２１３１０２２０３８５２
６Ｈｕａｎｇ２００９Ｃｈｉｎａ１１８２２８２３２１１９
７０ＧＡＴＡ２００９Ｊａｐａｎ３４１３２２６１２８３９９
８Ｃｈｅｎ２０１０Ｃｈｉｎａ１８１２５３２２９１５
９Ｈｗａｎｇ２０１０Ｋｏｒｅａ２９９５２７６１２２６５
１０Ｋｅｓｔｅｒ２０１１ＵＫ３５１３４３１５２９３８１１
ｒｕｎ；
ｄａｔａｍｕｌｔｉｖａｒ；ｓｅｔｄａｔａ；
ｉｆａｒｅａ＝”ＵＳ”ｏｒａｒｅａ＝”ＵＫ”ｔｈｅｎｚ＝ｌ；ｅｌｓｅｚ＝０；
ｒｅｎａｍｅｈｅａ—ｈｉ＝ｎ０１ｈｅａｍｉｄ＝ｎ０２ｈｅａ—ｌｏｗ＝ｎ０３ｈｅａ＿ｎ＝ｎ０ｄｉｓ—ｈｉ＝ｎｌｌｄｉｓ＿ｍｉｄ＝ｈｉ２ｄｉｓ—ｌｏｗ＝ｎ１３ｄｉｓ＿ｎ＝ｎｌ；
ｒｕｎ；
ｐｒｏｃｎｌｍｉｘｅｄｄａｔａ＝ｍｕｌｔｉｖａｒｔｅｃｈ＝ｎｅｗｒａｐ；
ｐａｒｒｎｓｍａ＝２．８ｍｘｌ＝－３．６７ｒｅｘ２＝．１．２７ｂ＝ｌｃ＝０ｖａ＝０．３ｃｏｖａｄ－一０．２ｖｄ＝０．６ｖｅ＝０．０５；
ｂｏｕｎｄｓｍｘｌ—ｒｅｘ２＜＝０；
ｙｌ＝ａ＋ｂ＋ｘｌ＋ｃ４ｚ；
ｙ２＝ａ＋ｂ害ｘ２＋ｃ８ｚ：
ｐＯｌ＝ｌ／（１＋ｅｘｐ（一（ｘ１）））；
ｐ０２＝ｌ／（１＋ｅｘｐ（一（ｘ２）））－１／（１＋ｅｘｐ（一（ｘ１）））；
ｐ０３＝ｌ－ｌ／（１＋ｅｘｐ（－（Ｘ２）））；
ｐｌｌ＝ｌ／（１＋ｅｘｐ（－（ｙ１）））；
ｐ１２＝ｌ／（１＋ｅｘｐ（一（ｙ２）））－１／（１＋ｅｘｐ（－（ｙ１）））；
ｐ１３－－１－１／（１＋ｅｘｐ（一（ｙ２）”；
ｉｆ（ｐＯｌ“＝Ｏａｎｄｐ０２“＝Ｏａｎｄｐ０３“＝Ｏａｎｄｐ１１“＝０ａｎｄｐ１２“＝Ｏａｎｄｐ１３“＝ｏ）ｔｈｅｎｌｏｇｌ＝ｎ０１＋ｌｏｇ（ｐ０１）＋ｎ０２４ｌｏｇ（ｐ０２）＋ｎ０３＋ｌｏｇ（ｐ０３）＋ｎｌｌ４ｌｏｇ（ｐ１１）＋ｎ１２４ｌｏｇ（ｐ１２）＋ｎ１３４ｌｏｇ（ｐ１３）；
ｅｌｓｅｌｏｇｌ＝１４ｅｘｐ（－３）；
ｍｏｄｅｌｎ０１～ｇｅｎｅｒａｌＯ０９１）；
ｒａｎｄｏｍａｘｌｘ２一ｎｏｒｍａｌ（［ｍａ，ｍｘｌ，ｍｘ２］，
【ｖａ，
ｃｏｖａｄ，ｖｄ＋ｖｅ，
ｃｏｖａｄ，ｖｄ，ｖｄ＋ｖｅ］）
ｓｕｂｊｅｃｔ＝ｓｔｕｄｙｏｕｔ＝ｔｏｒｅｒｏ；
ｅｓｔｉｍａｔｅ”Ｓｅｎｌ”ｌ／（１＋ｅｘｐ（一（ｍａ＋０．９２＋ｍｘｌ）））；
ｅｓｔｉｍａｔｅ”Ｓｐｅｌ”１／（１＋ｅｘｐ（ｍｘｌ））；
ｅｓｔｉｍａｔｅ”ＤＯＲｌ”ｅｘｐ（ｍａ＋０．９２＋ｍｘｌ—ｍｘｌ）；
ｅｓｔｉｍａｔｅ”Ｓｅｎ２”ｌ／（１＋ｅｘｐ（一（ｍａ＋０．９２４ｍｘ２）”；
ｅｓｔｉｍａｔｅ”Ｓｐｅ２”ｌ／（１＋ｅｘｐ（ｍｘ２））；
ｅｓｔｉｍａｔｅ”ＤＯＲ２”ｅｘｐ（ｍａ＋０．９２＋ｍｘ２－ｍｘ２）；
ｒｕｎ；
结果显示，最大似然法估计随机变量ａ均值为２．６６（０．３２），参数卢为ｏ．８７（０．０６），ｙ为・０．６２（０．３１），若以Ｋｏｂａｙａｓｈｉ评分＝７分为阳性判定标准，则其综合灵敏度为０．３８，９５％ＣＩ为（０．２０，０．５７）；特异度为Ｏ．９７，９５％ＣＩ为（０．９６，０．９９）；诊断比值比为２２．８，９５％Ｃ１为（５．Ｏ，４０．５）。若以Ｋｏｂａｙａｓｈｉ评分＝４分为阳性判定标准，则其综合灵敏度为Ｏ．７８，９５％ＣＩ为（Ｏ．６７，０．９０）；特异度为０．８３，０９５％ＣＩ（０．７４，０．９２）；诊断比值比为１７．５，９５％ＣＩ为（３．７，３１．３）。曲线下面积采用梯形面积法计算，每研究的ＲＯＣ曲线，ＳＲＯＣ曲线及其置信带、预测带可根据估计的非线性混合模型转换而成，ＳＡＳ程序较复杂，这里就不详细介绍。
三、讨论
诊断试验中较常见的是一个闽值，获得的是一对灵敏度和特异度。对于每个研究只有一个阈值的ｍｅｔａ分析，采用双变量随机效应模型结合传统的ＳＲＯＣ法，可以得到较准确而可靠的分析结果。但每个研究如果有多个阂值，则可获得多对灵敏度和特异度。此时，若通过降维的方法采用双变量随机效应模型进行评价，则会造成数据信息的损失。为此，国内外的研究者们提出了针对有序检测结果的潜在正态分布模型法、采用分层结合有序回归的ｍｅｔａ分析法及针对癌症分级的多变量随机效应模型等一系列方法［３１。但这些方法并不是双变量随机效应模型的直接扩展，应用存在一定的局限性。本文所介绍的多变量随机效应模型在原理和方法上可看作是双变量随机效应模型的扩展和延伸，可通过ＳＡＳ的ＰＲＯＣＮＬＭＩＸＥＤ过程实现，有着广泛的应用前景。
本模型针对每个研究有相同数目的阈值，对于阈值数目不同的情况，虽然模型允许缺失数据的存在，但是否可以无需校正直接应用该模型，还有待进行进一步探讨和验证。
参考文献
１．ＬｉｔｔｅｎｂｅｒｇＢ，ＭｏｓｅｓＬＥ．Ｅｓｔｉｍａｔｉｎｇｄｉａｇｎｏｓｔｉｃ・ａｃｃｕｒａｃｙｆｒｏｍｍｕｌｔｉｐｌｅｃｏｎｆｌｉｃｔｉｎｇｒｅｐｏｒｔｓ－ａｎｅｗｍｅｔａ－ａｎａｌｙｔｉｃｍｅｔｈｏｄ．ＭｅｄｉｃａｌＤｅｃｉｓｉｏｎＭａｋｉｎｇ．１９９３，１３：３１３－３２１．
２．ＡｒｅｎｄｓＬＲ，ＨａｍｚａＴＨ，ｖａｌｌＨｏｕｗｅｌｉｎｇｅｎＪＣ，ｅｌａ１．Ｂｉｖａｒｉａｔｅｒａｎｄｏｍｅｆｆｅｃｔｓｍｅｔａ－ａｎａｌｙｓｉｓ
Ｄｅｃｉｓｉｏｎｏｆ—ＲＯＣｃ—ｕｒｖｅｓ．ＭｅｄｉｃａｌＭａｋｉｎｇ，２００８，２８（５）：６２１・６３８．
３．ＡｒｅｎｄｓＬＲ，Ｈａｍｚａ＂１１－１，ＶａｎＨｏｕｗｅｌｉｎｇｅｎＪＣ，ｅｔａ１．Ｍｕｌｌｉｖａｆｉａｔｅｒａｎｄｏｍｅｆｆｅｃｔｓｒｅｃｔａ－ａｎａｌｙｓｉｓｏｆｄｉａｇｎｏｓｔｉｃｔｅｓｔｓｗｉｔｈｍｕｔｉｐｌｅｔｈｒｅｓｈｏｌｄｓ．ＢＭＣＭｅｄｉｃａｌＲｅｓｅａｒｃｈＭｅｔｈｏｄｏｌｏｇｙ２００９，９：７３．
４．ＴｏｈｒｕＫｏｂａｙａｓｈｉ，ＹｏｓｈｉｎａｒｉＩｎｏｕｅ，ＫａｚｕｏＴａｋｅｕｃｈｉ，ｅｔａ１．ＰｒｅｄｉｃｔｉｏｎｏｆＩｎｔｒａｖｅｎｏｕｓＩｍｍｕｎｏｇｌｏｂｕｌｉｎＵｎｒｅｓｐｏｎｓｉｖｅｎｅｓｓｉｎＰａｔｉｅｎｔｓｗｉｔｈＫａｗａＳａｋｉＤｉｓｅａｓｅ．Ｃｉｒｃｕｌａｔｉｏｎ，２００１６。１１３：２６０６．２６１２．１９１
多变量随机效应模型及其SAS实现
作者：
作者单位：刘文华， 吴家利， 杨杨， 宋婷婷， 章顺悦， 陈远方， 尹平华中科技大学同济医学院公共卫生学院流行病与卫生统计学系 430030
引用本文格式：刘文华.吴家利.杨杨.宋婷婷.章顺悦.陈远方.尹平 多变量随机效应模型及其SAS实现[会议论文] 2011
２０１１年中国卫生统计学年会论文集
多变量随机效应模型及其ＳＡＳ实现
刘文华，吴家利，杨杨，宋婷婷，章顺悦，陈远方，尹平△
华中科技犬学尉济医学院公共卫生学院流行病与卫生统计学系（４３００３０）
摘要目的介绍多变量随机效应模型及其参数估计，促进其在多个阈值诊断试验ｍｅｔａ分析中的应用。方法对双变量随机效应模型进行扩展和延伸，构建多变量随机效应模型，采月｛｝最大似然法，并利用ＳＡＳ中的ＰＲＯＣＮＬＭＩＸＥＤ过程进行统计实现。结果采用多变量随机效应模型可得到每一闽值的合并灵敏度、特异度、诊断比值比，及曲线下面积ＡＵＣ等评价指标。亦可绘制ＳＲＯＣ曲线及其置信带、预测带等。结论多变量随机效应模型能够充分利用各研究的原始数据信息，得到稳定可靠的结果，有着较强的应用价值。
关键词多变量随机效应模型闽值诊断试验ｍｅｔａ分析
Ｍｕｌｔｉｖａｒｉａｔｅｒａｎｄｏｍ
ＬｉｕｅｆｆｅｃｔｓｍｏｄｅｌａｎｄｉｔｓｓＡＳｐｒｏｇｒａｍｓＷｅｎｈｕａ，ＷｕＪｉａｌｉ，ＹａｎｇＹａｎｇ，ＳｏｎｇＴｉｎｇｔｉｎｇ，ＺｈａｎｇＳｈｕｎｙｕｅ，ＣｈｅｎＹｕａｎｇｆａｎｇ，ＹｉｎＰｉｎｇ
ＤｅｐａｒｔｍｅｎｔｏｆＥｐｉｄｅｍｉｏｌｏｇｙａｎｄＨｅａｌｔｈＳｔａｔｉｓｔｉｃｓ，ＳｃｈｏｏｌｏｆＰｕｂｌｉｃ
ｏｆＳｃｉｅｎｃｅａｎｄＨｅａｔｈ，ＴｏｎｇｉｉＭｅｄｉｃａｌＣｏｌｌｅｇｅｏｆＨｕａｚｈｏｎｇＵｎｉｖｅｒｓｉｔｙＴｅｃｈｎｏｌｏｇｙ｛４３００３０），Ｗｕｈａｎ
ｔｏＡｂｓｔｒａｃｔ０Ｍｅｃｔｉｖｅ：Ｔｈｅｓｔｕｄｙｗａｓｃｏｎｄｕｃｔｅｄｉｎｔｒｏｄｕｃｅｔｈｅｍｕｌｔｉｖａｒｉａｔｅｒａｎｄｏｍｅｆｆｅｃｔｓｍｏｄｅｌａｎｄｉｔｓ
ｐａｒａｍｅｔｅｒｅｓｔｉｍａｔｉｏｎ，ａｎｄｔｏｆａｃｉｌｉｔｙｉｔｓａｐｐｌｉｃａｔｉｏｎｓｉｎｍｅｔａ－ａｎａｌｙｓｉｓｏｆｄｉａｇｎｏｓｔｉｃｔｅｓｔｓｗｉｔｈｍｕｌｔｉｐｌｅｔｈｒｅｓｈｏｌｄｓ．Ｍｅｔｈｏｄｓ：Ｔｈｉｓｐａｐｅｒｅｘｐａｎｄｅｄａｎｄｅｘｔｅｎｄｅｄｔｈｅｂｉｖａｒｉａｔｅｒａｎｄｏｍｅｆｆｅｃｔｓｍｏｄｅｌｔｏｄｅｖｅｌｏｐｔｈｅｍｕｌｔｉｖａｒｉａｔｅｒａｎｄｏｍｅｆｆｅｃｔｓｍｏｄｅｌ，ａｎｄｉｍｐｌｅｍｅｎｔｅｄｉｔｉｎｔｈｅＳＡＳＰｒｏｃＮＬＭＩＸＥＤｐｒｏｃｅｄｕｒｅｗｉｔｈｍａｘｉｍｕｍｌｉｋｅｌｉｈｏｏｄ．Ｒｅｓｕｌｔｓ：Ｔｈｅｍｕｌｔｉｖａｒｉａｔｅｒａｎｄｏｍｅｆｆｅｃｔｓｍｏｄｅｌａｌｌｏｗｓｃａｌｃｕｌａｔｉｎｇｔｈｅｐｏｏｌｅｄ
ａｒｅａｓｅｎｓｉｔｉｖｉｔｙ，ｓｐｅｃｉｆｉｃｉｔｙ，ｄｉａｇｎｏｓｔｉｃｏｄｄｓｒａｔｉｏｆｏｒｅａｃｈｔｈｒｅｓｈｏｌｄａｎｄｔｈｅ
ｔｈｅＳＲＯＣｃｕｒｖｅｕｎｄｅｒｔｈｅｃｕｒｖｅ，ｗｈｉｌｅｄｒａｗｉｎｇｗｉｔｈｉｔｓｃｏｎｆｉｄｅｎｃｅｂａｎｄａｎｄｐｒｅｄｉｃｔｉｏｎｂａｎｄ．Ｃｏｎｃｌｕｓｉｏｎｓ：Ｔｈｅｍｕｌｔｉｖａｒｉａｔｅｒａｎｄｏｍ
ｆｕｌｌｙｅｘｐｌｏｉｔｔｈｅｄａｔａ，ｏｂｔａｉｎｓｔａｂｌｅａｎｄｒｅｌｉａｂｌｅｒｅｓｕｌｔｓ，ａｎｄｈａｖｅａｅｆｆｅｃｔｓｍｏｄｅｌｃａｎｇｏｏｄａｐｐｌｉｃａｔｉｏｎ．
Ｋｅｙｗｏｒｄｓ：ｍｕｌｔｉｖａｒｉａｔｅｒａｎｄｏｍｅｆｆｅｃｔｓｍｏｄｅｌ；ｔｈｒｅｓｈｏｌｄ；ｄｉａｇｎｏｓｔｉｃｔｅｓｔｓ；ｍｅｔａ－ａｎａｌｙｓｉｓ
近年来，诊断试验ｍｅｔａ分析在医学研究中应用越来越广泛，其统计学方法主要针对每个研究只有一个诊断阈值，提供一对灵敏度和特异度或一个四格表的资料，一般采用分别合并灵敏度、特异度，或采用固定或随机效应模型联合灵敏度和特异度进行分析【１，２１。但在实际工作中会经常遇到检测结果有多个诊断阈值的情形，针对这类问题，以往ｍｅｔａ分析常用的方法是将其转换为二分类，即每个研究根据一个闽值得到一对灵敏度和特异度，再应用ＳＲＯＣ方法或双变量随机效应模型进行分析。而这种“降维”的方法会损失数据信息，降低数据的利用度。解决这类多阂值问题的有效方法是多变量随机效应模型，其可以充分利用原始数据信息，获得更为可靠的结论。国际上已进行了一些研究
【３Ｊ’但国内的相关研究报道甚少。为此，本文对多变量随机效应模型进行介绍，以期能更好地促进国内该方法在多个阂值诊断试验ｍｅｔａ分析中的应用。△通讯作者：尹平
２０１１年中国卫生统计学年会论文集
一、模型与参数估计
多变量随机效应模型可看作双变量随机效应模型的扩展和延伸．在双变量随机效应模型中，灵敏度和１．特异度经Ｌｏｇｉｔ转换后的真实值服从双变量正态分布，灵敏度和１一特异度的真实值分别服从二项分布。而在多变量随机效应模型中，针对多个阈值，可得到相应的多对灵敏度和１一特异度，经Ｌｏｇｉｔ转换后其真实值服从多变量正态分布，而检测试验取每个结果的频数在疾病人群和非疾病人群中均服从多项分布。根据乘积多项分布理论，可获得其联合概率，从而构造出似然函数，并采用最大似然法进行参数估计。
假设经纳入标准和排除标准筛选，１个研究，每个研究有ｇ－１个阈值，则检测指标可分为ｇ个组段或ｇ个有序变量。对于研究ｆ，按检测结果从高到低，依次列出各组段或等级的观测数，得到一个２ｘｇ列联表，如表１所示。
表１第ｆ研究金标准与检测试验结果的２×ｇ频数分布
列联表中每个格子的观测数目为随机变量，在病例人群中服从参数为ｎｌｉ，刀，∥…，兀岵的多项分布，在非病例人群中服从参数为，＇ｌＯｉ，巧０ｆｌ，…，石ｏｆｇ的多项分布。
以测定值大于阈值判断为检测试验阳性，则对于研究；，阈值，，其真阳性率研碾；，；翌，假。
ｎｌｉ
眦率脚。；擎，将珊玎，肿。进弛ｇ；ｔ转换，参数ｙ俩郴ｕ表示撇盯，脚。转换后的总体真实值，则殉，却，问线性关系为：
ＹⅡ＝口ｆ＋肛目
而ｘｄ，可用如下模型表示：
ｘ盯＝ｘ，＋６ｆ＋Ｐ口（１）（２）
式中，口；表示各研究常数项参数；斜率卢为对称性参数，各研究斜率固定，卢＝１表示ＳＲＯＣ曲线围绕“灵敏度：特异度，，直线对称；；，表示对于诊断阈值．『，‰的平均值；６ｉ表示系统误差，：白表示残差。该模型还作如下假设：Ｏｔｉ。Ⅳ五，仃。２），６，一Ｎ（Ｏ，仃ｄ２），ｅ．｝：－Ｎ（Ｏ，仃。２），仃。＝０，仃＆＝０ａ因此模型也可写成如下矩阵形式：
——————————————————————————————————————————————，——————————————————一一．一！！！！兰！望墨竺篁盐堂兰全笙壅篓一——
…
ａｄｒ口ｄｄｒ甜仃神
ｉ
２２…
ｘｉｌｏ口＋０ｅ仃６２仃６２
：
●
～Ｎｆ３）
Ｘ每仃６‘．Ｄｄ
●
：骥二
Ｚｉ，ｇ一１ｚｇ—ｌ』Ｉ口面…２
仃６口６…口ｄ２＋仃Ｐ２
采用最人似然法进行参数估计，根据乘积多项分布原理，其联合概率密度函数为：厂（，ｚ。玎，，ｚ・玎ｋ。ｒ，九，ｚ，万叫，万ｔ玎）２ｌ＝ｌｎ‘ｉ丽ｎｏｉ！．ｎｌｉ！导ｇｊ万叫疗叫万－玎”ｌ口）ｃ４，
ｅＸＯ，＿『＝１
１４－ｅＸ，Ｊ
其中，石叫＝‰（月孵）＝三一』兰，＿『．２，…，ｇ－１，。（５）。
１４－ｅｂ１＋ｅｘｉ＂卜１
１一』三，＿『；ｇ。。。
１＋ｅＸ，，Ｊ一１
ｇｙｇ
ａｊ；１
１＋ｅ’口
刀埘＝置牙Ｏｌ帮）＝』二一』二，－『：２，…，ｇ一１，。。。～（６）
１＋ｅｙ９１．＋ｅｙｉ．卜１
１一』兰，’『；ｇ。。。
１＋ｅｙｅ，卜１
骞鼽ｎｇ
Ｉ，ＪＪｃｌ－１，善荟～“。Ｊ—ＪＪ。Ｉ
因此，似然函数为：Ｌ＝Ｆ（臼卜ｏｉｊ＇ｎｌＯ）＝，Ｑ嘶，ｎ坷ｎｏｉ，ｎｉｌ∥Ｏｑ＇万坷）。同时，为得到非偏估计。采用迭代法不断修正最大似然估计值。
当然，也可将斜率参数芦设为随机变量，同对，为校正研究间协变量（ｚ）的影响，模型可写成：ｙ可＝绞ｆ＋屈ｘ口＋弦ａ进行参数估计时需增加仃芦２、仃肚２、仃彤２、），等４个参数，但较多的参数在一定程度上会影响模型的收敛性。
利用ＳＡＳ的ＰＲＯＣＮＬＭＩＸＥＤ过程可以完成多变量随机效应模型的拟合，得到每一阂值的合并灵敏度、特异度、诊断比值比，曲线下面积ＡＵＣ等评价指标。还可绘制ＳＲＯＣ曲线及其置信带、预测带，以及每个研究的ＲＯＣ曲线。１８９
二、ＳＡＳ实现
２００６年Ｋｏｂａｙａｓｈｉ等【４】提出了最新的川崎病并发冠脉损伤的高危评分体系，总积分为１１分，Ｏ～３分为低危组，４分以上为高危人群，７分以上为极高危人群。本文模拟Ｋｏｂａｙａｓｈｉ评分系统检测川崎病患者是否并发冠脉损伤的数据，将评分结果分为＝３分，４－６分，＝７分三组，按照川崎病和冠脉损伤诊断标准将川崎病ＩＶＩＧ无反应者分为ＣＡＬ（冠脉损伤）组和ＮＣＡＬ（无冠脉损伤）组，评价该评分体系的诊断准确性。ＳＡＳ过程如下。
ｄａｔａｄａｔａ；
ｌｅｎｇｔｈａｕｔｈｅｒａｒｅａ￥１０：
ｉｎｐｕｔｓｔｕｄｙａｕｔｈｅｒ￥ｙｅａｒａｒｅａ￥ｈｅａ＿ｈｉｈｅａ—ｍｉｄｈｅａ一１０Ｗｄｉｓ—ｈｉｄｉｓ—ｍｉｄｄｉｓ＿ｌｏｗ；
ｈｅａ＿＿ｎ＝ｈｅａ＿ｈｉ＋ｈｅａｍｉｄ＋ｈｅａｌｏｗ；
ｄｉｓ＿ｎ＝ｄｉｓ＿ｈｉ＋ｄｉｓ＿ｍｉｄ＋ｄｉｓ＿ｌｏｗ；
ｄａｔａｌｉｎｅｓ；
１Ｋｏｂａｙａｓｈｉ２００６ＵＳ６３９４８２１００１１６７８
２Ｙａｓｕｎｏｒｉ２００７ＵＳ９５９６０３２９３８８１
３Ｉｎｏｕｅ２００８ＵＳ１０８５７００３０５０９５
４Ｔａｋｅｕｃｈｉ２００８ＵＳ２８１５６２９９１８３０４
５Ａｅｒｔｇｅｅｒｔｓ２００９ＵＳ２１３１０２２０３８５２
６Ｈｕａｎｇ２００９Ｃｈｉｎａ１１８２２８２３２１１９
７０ＧＡＴＡ２００９Ｊａｐａｎ３４１３２２６１２８３９９
８Ｃｈｅｎ２０１０Ｃｈｉｎａ１８１２５３２２９１５
９Ｈｗａｎｇ２０１０Ｋｏｒｅａ２９９５２７６１２２６５
１０Ｋｅｓｔｅｒ２０１１ＵＫ３５１３４３１５２９３８１１
ｒｕｎ；
ｄａｔａｍｕｌｔｉｖａｒ；ｓｅｔｄａｔａ；
ｉｆａｒｅａ＝”ＵＳ”ｏｒａｒｅａ＝”ＵＫ”ｔｈｅｎｚ＝ｌ；ｅｌｓｅｚ＝０；
ｒｅｎａｍｅｈｅａ—ｈｉ＝ｎ０１ｈｅａｍｉｄ＝ｎ０２ｈｅａ—ｌｏｗ＝ｎ０３ｈｅａ＿ｎ＝ｎ０ｄｉｓ—ｈｉ＝ｎｌｌｄｉｓ＿ｍｉｄ＝ｈｉ２ｄｉｓ—ｌｏｗ＝ｎ１３ｄｉｓ＿ｎ＝ｎｌ；
ｒｕｎ；
ｐｒｏｃｎｌｍｉｘｅｄｄａｔａ＝ｍｕｌｔｉｖａｒｔｅｃｈ＝ｎｅｗｒａｐ；
ｐａｒｒｎｓｍａ＝２．８ｍｘｌ＝－３．６７ｒｅｘ２＝．１．２７ｂ＝ｌｃ＝０ｖａ＝０．３ｃｏｖａｄ－一０．２ｖｄ＝０．６ｖｅ＝０．０５；
ｂｏｕｎｄｓｍｘｌ—ｒｅｘ２＜＝０；
ｙｌ＝ａ＋ｂ＋ｘｌ＋ｃ４ｚ；
ｙ２＝ａ＋ｂ害ｘ２＋ｃ８ｚ：
ｐＯｌ＝ｌ／（１＋ｅｘｐ（一（ｘ１）））；
ｐ０２＝ｌ／（１＋ｅｘｐ（一（ｘ２）））－１／（１＋ｅｘｐ（一（ｘ１）））；
ｐ０３＝ｌ－ｌ／（１＋ｅｘｐ（－（Ｘ２）））；
ｐｌｌ＝ｌ／（１＋ｅｘｐ（－（ｙ１）））；
ｐ１２＝ｌ／（１＋ｅｘｐ（一（ｙ２）））－１／（１＋ｅｘｐ（－（ｙ１）））；
ｐ１３－－１－１／（１＋ｅｘｐ（一（ｙ２）”；
ｉｆ（ｐＯｌ“＝Ｏａｎｄｐ０２“＝Ｏａｎｄｐ０３“＝Ｏａｎｄｐ１１“＝０ａｎｄｐ１２“＝Ｏａｎｄｐ１３“＝ｏ）ｔｈｅｎｌｏｇｌ＝ｎ０１＋ｌｏｇ（ｐ０１）＋ｎ０２４ｌｏｇ（ｐ０２）＋ｎ０３＋ｌｏｇ（ｐ０３）＋ｎｌｌ４ｌｏｇ（ｐ１１）＋ｎ１２４ｌｏｇ（ｐ１２）＋ｎ１３４ｌｏｇ（ｐ１３）；
ｅｌｓｅｌｏｇｌ＝１４ｅｘｐ（－３）；
ｍｏｄｅｌｎ０１～ｇｅｎｅｒａｌＯ０９１）；
ｒａｎｄｏｍａｘｌｘ２一ｎｏｒｍａｌ（［ｍａ，ｍｘｌ，ｍｘ２］，
【ｖａ，
ｃｏｖａｄ，ｖｄ＋ｖｅ，
ｃｏｖａｄ，ｖｄ，ｖｄ＋ｖｅ］）
ｓｕｂｊｅｃｔ＝ｓｔｕｄｙｏｕｔ＝ｔｏｒｅｒｏ；
ｅｓｔｉｍａｔｅ”Ｓｅｎｌ”ｌ／（１＋ｅｘｐ（一（ｍａ＋０．９２＋ｍｘｌ）））；
ｅｓｔｉｍａｔｅ”Ｓｐｅｌ”１／（１＋ｅｘｐ（ｍｘｌ））；
ｅｓｔｉｍａｔｅ”ＤＯＲｌ”ｅｘｐ（ｍａ＋０．９２＋ｍｘｌ—ｍｘｌ）；
ｅｓｔｉｍａｔｅ”Ｓｅｎ２”ｌ／（１＋ｅｘｐ（一（ｍａ＋０．９２４ｍｘ２）”；
ｅｓｔｉｍａｔｅ”Ｓｐｅ２”ｌ／（１＋ｅｘｐ（ｍｘ２））；
ｅｓｔｉｍａｔｅ”ＤＯＲ２”ｅｘｐ（ｍａ＋０．９２＋ｍｘ２－ｍｘ２）；
ｒｕｎ；
结果显示，最大似然法估计随机变量ａ均值为２．６６（０．３２），参数卢为ｏ．８７（０．０６），ｙ为・０．６２（０．３１），若以Ｋｏｂａｙａｓｈｉ评分＝７分为阳性判定标准，则其综合灵敏度为０．３８，９５％ＣＩ为（０．２０，０．５７）；特异度为Ｏ．９７，９５％ＣＩ为（０．９６，０．９９）；诊断比值比为２２．８，９５％Ｃ１为（５．Ｏ，４０．５）。若以Ｋｏｂａｙａｓｈｉ评分＝４分为阳性判定标准，则其综合灵敏度为Ｏ．７８，９５％ＣＩ为（Ｏ．６７，０．９０）；特异度为０．８３，０９５％ＣＩ（０．７４，０．９２）；诊断比值比为１７．５，９５％ＣＩ为（３．７，３１．３）。曲线下面积采用梯形面积法计算，每研究的ＲＯＣ曲线，ＳＲＯＣ曲线及其置信带、预测带可根据估计的非线性混合模型转换而成，ＳＡＳ程序较复杂，这里就不详细介绍。
三、讨论
诊断试验中较常见的是一个闽值，获得的是一对灵敏度和特异度。对于每个研究只有一个阈值的ｍｅｔａ分析，采用双变量随机效应模型结合传统的ＳＲＯＣ法，可以得到较准确而可靠的分析结果。但每个研究如果有多个阂值，则可获得多对灵敏度和特异度。此时，若通过降维的方法采用双变量随机效应模型进行评价，则会造成数据信息的损失。为此，国内外的研究者们提出了针对有序检测结果的潜在正态分布模型法、采用分层结合有序回归的ｍｅｔａ分析法及针对癌症分级的多变量随机效应模型等一系列方法［３１。但这些方法并不是双变量随机效应模型的直接扩展，应用存在一定的局限性。本文所介绍的多变量随机效应模型在原理和方法上可看作是双变量随机效应模型的扩展和延伸，可通过ＳＡＳ的ＰＲＯＣＮＬＭＩＸＥＤ过程实现，有着广泛的应用前景。
本模型针对每个研究有相同数目的阈值，对于阈值数目不同的情况，虽然模型允许缺失数据的存在，但是否可以无需校正直接应用该模型，还有待进行进一步探讨和验证。
参考文献
１．ＬｉｔｔｅｎｂｅｒｇＢ，ＭｏｓｅｓＬＥ．Ｅｓｔｉｍａｔｉｎｇｄｉａｇｎｏｓｔｉｃ・ａｃｃｕｒａｃｙｆｒｏｍｍｕｌｔｉｐｌｅｃｏｎｆｌｉｃｔｉｎｇｒｅｐｏｒｔｓ－ａｎｅｗｍｅｔａ－ａｎａｌｙｔｉｃｍｅｔｈｏｄ．ＭｅｄｉｃａｌＤｅｃｉｓｉｏｎＭａｋｉｎｇ．１９９３，１３：３１３－３２１．
２．ＡｒｅｎｄｓＬＲ，ＨａｍｚａＴＨ，ｖａｌｌＨｏｕｗｅｌｉｎｇｅｎＪＣ，ｅｌａ１．Ｂｉｖａｒｉａｔｅｒａｎｄｏｍｅｆｆｅｃｔｓｍｅｔａ－ａｎａｌｙｓｉｓ
Ｄｅｃｉｓｉｏｎｏｆ—ＲＯＣｃ—ｕｒｖｅｓ．ＭｅｄｉｃａｌＭａｋｉｎｇ，２００８，２８（５）：６２１・６３８．
３．ＡｒｅｎｄｓＬＲ，Ｈａｍｚａ＂１１－１，ＶａｎＨｏｕｗｅｌｉｎｇｅｎＪＣ，ｅｔａ１．Ｍｕｌｌｉｖａｆｉａｔｅｒａｎｄｏｍｅｆｆｅｃｔｓｒｅｃｔａ－ａｎａｌｙｓｉｓｏｆｄｉａｇｎｏｓｔｉｃｔｅｓｔｓｗｉｔｈｍｕｔｉｐｌｅｔｈｒｅｓｈｏｌｄｓ．ＢＭＣＭｅｄｉｃａｌＲｅｓｅａｒｃｈＭｅｔｈｏｄｏｌｏｇｙ２００９，９：７３．
４．ＴｏｈｒｕＫｏｂａｙａｓｈｉ，ＹｏｓｈｉｎａｒｉＩｎｏｕｅ，ＫａｚｕｏＴａｋｅｕｃｈｉ，ｅｔａ１．ＰｒｅｄｉｃｔｉｏｎｏｆＩｎｔｒａｖｅｎｏｕｓＩｍｍｕｎｏｇｌｏｂｕｌｉｎＵｎｒｅｓｐｏｎｓｉｖｅｎｅｓｓｉｎＰａｔｉｅｎｔｓｗｉｔｈＫａｗａＳａｋｉＤｉｓｅａｓｅ．Ｃｉｒｃｕｌａｔｉｏｎ，２００１６。１１３：２６０６．２６１２．１９１
多变量随机效应模型及其SAS实现
作者：
作者单位：刘文华， 吴家利， 杨杨， 宋婷婷， 章顺悦， 陈远方， 尹平华中科技大学同济医学院公共卫生学院流行病与卫生统计学系 430030
引用本文格式：刘文华.吴家利.杨杨.宋婷婷.章顺悦.陈远方.尹平 多变量随机效应模型及其SAS实现[会议论文] 2011