如何用SAS软件正确分析生物医学科研资料XXII. 结果变量为二值变量的高维列联表资料的统计分析与SAS软件实现（二）

鲍晓蕾，胡良平



如何用SAS软件正确分析生物医学科研资料XXII. 结果变量为二值变量的高维列联表资料的统计分析与SAS软件实现（二）

鲍晓蕾，胡良平

100850 北京，军事医学科学院生物医学统计学咨询中心

编者按

生物统计学是生物学领域科学研究和实际工作中必不可少的工具，在分子生物学迅速发展的今天，生物统计学更显示出了它的重要性。实验设计与数据统计分析是现代生物学的基石，是生物学研究者检验假说、寻找模式、建立生物学理论的有利工具，也是生物学研究者探索微观和宏观生物世界的必备基础知识。对于每天甚至是每时每刻涌现的大量的、以天文数字计量的分子遗传数据，必须借助统计学知识加以分析处理，才能从中获得有意义的信息。“生物多样性数据分析”是开展生物多样性研究的一个重要方面，数据分析能力的高低极大地影响着我们对各种生态学现象认识的深度和广度。现在，电子计算机的普及使得生物统计分析过程大大简化，生物统计分析软件包的普及将生物统计学从统计学家的书本里解放了出来，简化了生物统计分析过程，使之成为生物学研究者的常用工具。本刊特邀军事医学科学院生物医学统计学咨询中心主任胡良平教授，以“如何用SAS软件正确分析生物医学科研资料”为题，撰写系列统计学讲座，希望该系列讲座能对生物医学科研工作者有所帮助。

当观测结果是定性资料时，人们习惯将资料整理成列联表形式。比如“2´2 列联表资料”、“R´C 列联表资料”和“高维列联表资料”等。所谓高维列联表，也就是表中涉及到的定性变量的个数≥ 3。对于高维列联表资料，根据结果变量的性质可将高维列联表分为以下三类：一是结果变量为二值变量的高维列联表；二是结果变量为多值有序变量的高维列联表；三是结果变量为多值名义变量的高维列联表。在上一期中，我们已经介绍了用CMH2检验和加权2检验处理结果变量为二值变量的高维列联表资料，本文将继续介绍结果变量为二值变量的高维列联表资料的其他分析方法及 SAS 实现。

处理不同的高维列联表资料所用的统计分析方法不尽相同，对于结果变量为二值变量的高维列联表资料可使用 CMH2检验、多重 logistic 回归分析、对数线性模型等统计分析方法。本文重点介绍对数线性模型和多重 logistic 回归分析如何处理结果变量为二值变量的高维列联表资料。

1 对数线性模型

在上一期中，我们介绍了用 CMH2检验和加权2检验来处理结果变量为二值变量的高维列联表资料，这两种方法的本质都是控制其中一个原因变量对结果的影响，从而分析另一个原因变量的作用。这样的做法往往不够准确，特别是当被合并的变量是一个混杂因素，即与留下的变量之间不独立时，极易得出错误的结论。那么，当同时分析各变量的主效应及变量间的交互效应时，对数线性模型便显示出其巨大的潜力。

对数线性模型是分析高维列联表资料行之有效的方法，最先由 Yule、Bartlett 利用 Yule（1900 年）定义的交叉乘积比分析三维交互作用，然后在 Kullback（1968 年）引入方差分析的思想发展而来[1]。

对数线性模型的基本思想：把各分组变量（包括自变量和因变量）水平组合下期望（理论频数）的自然对数表示为各分组变量及其交互作用项和定性结果变量的线性函数，通过迭代计算求得模型中参数的估计值，进而运用方差分析的思想检验各主效应和交互作用的效应大小[1]。

主效应与交互效应：对数线性模型主要考察各分类变量间的交互作用（关联性），对主效应的分析相当于单变量分析。交互效应按其因素多少可分为两变量间交互效应和多变量间交互效应，它们依次又被称为一阶交互效应、二阶交互效应……，以此类推。

对数线性模型参照线性模型和方差分析的思想，将列联表的总效应（期望频数）取自然对数后分解成主效应和各因素间的交互效应。对数线性模型虽与线性模型和方差分析类似，但又有其自身特点。简单起见，以三维列联表为例介绍对数线性模型的构建原理，对于四维乃至更高维度的列联表，不失一般性。

设 3 个分类变量分别为 A、B、C，则包含所有交互效应的饱和模型为：

只包含 A 与 B、A 与 C 交互效应的谱系模型为：

表 1 对数线性模型参数意义与个数

2 多重 logistic 回归分析

Logistic 回归分析是 1970 年 Cox 提出的[1]。设（y = 1│X）（简记为）表示暴露因素为 X 时个体发病的概率。称发病的概率与未发病的概率 1-之比为“优势（odds）”。对作 logit 变换，logit() 定义为优势之对数（log odds），即：

因此，多重 logistic 回归方程的线性表达式定义为：

logit() =+βx+ βx+ … +βx

logit()与各因素间呈线性关系，β为 logistic 回归的回归系数，表示在其他变量都固定的情况下，对即 logit() 影响的大小。取β的反对数可得优势比= exp(β)，表示其他变量不变的情况下，每变化一个单位时所引起的优势比的自然对数改变量。当某种疾病的发病率或死亡率很低时，≈（即相对危险度）[1]。

同多重线性回归分析一样，当比较暴露因素对因变量相对贡献大小时，由于各自变量取值单位不同，也不能用回归系数的大小直接进行比较，而需用标准化回归系数来比较。建立 logistic 回归方程的过程也就是求常数项及各回归系数β的过程。

Logistic 回归分析按照因变量（即反应变量）的类型可分为：因变量为二值变量的 logistic 回归分析、因变量为多值有序变量的 logistic 回归分析、因变量为多值无序（或名义）变量的 logistic 回归分析；logistic 回归按研究设计类型可分为非条件（成组设计）logistic 回归分析和条件（匹配设计）logistic 回归分析。

对于 logistic 回归方程的拟合一般利用 SAS 中logistic 过程完成，一些较复杂的logistic 回归分析则需要通过调用 SAS 中的 catmod 过程或 phreg 过程完成。本文仅介绍调用 SAS 中的 logistic 过程完成一般的 logistic 回归分析。

下面仍沿用上一期的实例，向读者介绍如何通过 SAS 软件使用对数线性模型和多重logistic 回归模型处理结果变量为二值变量的高维列联表资料。

【例 1】 为调查眼科门诊患者的用药依从性及其影响因素，并提出相应对策，提高眼科门诊药房的药学服务水平，在 5 家医院的眼科门诊患者中，随机发放 610 份调查问卷，进行汇总，数据见表 2。

表 2 不同告知者、不同告知程度间患者的依从情况

⑴对数线性模型

SAS 程序如下，程序名为li3.sas。

data a1;do a=1 to 2;do b=1 to 3;do c=1 to 2;input f @@;output;end; end; end; cards;240 24 202 60 31 16 188 19 238 60 47 21;run;ods html;proc catmod;weight f;model a*b*c=_response_/noglsnoparm noresponse ml;loglin a|b|c;run;ods html close;

程序说明：a 表示告知者，a = 1 表示医生，a = 2 表示药师；b 表示告知程度，b = 1 表示详细告知，b = 2 表示告知，b = 3 表示没告知；c 表示是否依从，c = 1 表示依从，c = 2 表示不依从。过程步调用 catmod 过程，若需输出对数线性模型中的各参数的估计值和误差，只需去掉 model 语句中的 noparm 选项。ml 选项表示采用最大似然法进行参数估计，nogls 表示不采用加权最小二乘法计算。Noresponse 表示不输出反应矩阵。Loglin 定义线性模型的效应项。

SAS 程序运行结果：

Maximum likelihood analysis of variance

SourceDFChi-SquarePr > ChiSq a1 0.450.5014 b2173.20< .0001 a*b2 5.420.0665 c1241.62< .0001 a*c1 0.290.5929 b*c2 42.45< .0001 a*b*c2 0.210.8988 Likelihood Ratio0..

结果显示，b 与 c 的交互作用存在统计学意义（< 0.0001），说明患者被告知的程度不同，其依从性不同。而因素 a（告知者）对结果的影响没有统计学意义，此时把不同告知者的数据进行简单合并后进行分析是可以接受的。若 a*b*c 这个高阶交互作用有统计学意义，则说明不同告知者、不同告知程度下患者的依从性有差别，此时不适合用压缩法消除“告知者”的影响而研究“告知程度”对患者依从性的作用。

专业结论：患者被告知的程度不同，依从性不同（实际资料显示：被详细告知者依从性高）；而不同的告知者对患者依从性的影响没有统计学意义。

⑵多重 logistic 回归分析

程序如下，设程序名为 li4.sas。

data a2;do a=1 to 2;do b=1 to 3;do c=1 to 2;input f @@;output;end; end; end; cards;240 24 202 60 31 16 188 19 238 60 47 21;run;ods html;proc logistic;class b/param=ref;freq f; model c=a b a*b/selection=forward;oddsratio b;run;proc logistic;class b/param=ref;freq f;model c=a b a*b/selection=backward;oddsratio b;run;proc logistic;class b/param=ref;freq f;model c=a b a*b/selection=stepwise;oddsratio b;run;ods html close;

程序说明：数据步与程序 li3.sas 一样，不再重复。调用 logistic 过程进行结果变量为二值变量的多重 logistic 回归分析。Class 语句指定对变量 b 赋哑变量，param = ref表示以其中的一个水平为基准建立哑变量，默认以输入的最后一个水平为基准，本例则以变量b 的最后一个水平（即“未告知”）为基准建立哑变量。语句freq f 指定f 变量为频数变量。通过model 语句进行建模，等号前面表示因变量，等号后面表示自变量，其中a*b 表示变量a 和b 的交互作用，斜杠后面的selection 选项表示进行变量筛选。程序中的三个logistic 过程分别表示用forward（前进法）、backward（后退法）和stepwise（逐步筛选法）进行变量筛选，运用多种变量筛选方法的目的是为了使筛选结果更可信。第一个logistic 过程中selection = forward 表示用前进法进行变量筛选，其中forward 法默认的入选变量进入方程的检验水准为0.05，若想自定义该水准，可在该语句后加入选项sle = xx 来指定；第二个logistic 过程中selection = backward 表示用后退法进行变量筛选，其默认的从方程中剔除变量的检验水准为0.05，若想自定义该水准，可在该语句后加入选项sls = xx 来指定；第三个logistic 过程中selection = stepwise 表示用逐步筛选法进行变量筛选，其默认的入选水准和排除水准分别为0.05，若想自定义该值，可在该语句后加入sle = xx sls = xx 来指定，前者表示入选水准，后者表示剔除水准，里面的“xx”实际使用时都应给予具体的数值，如0.3、0.1、0.03、0.01 等。oddsratio b 表示输出变量b 任意两个水平间的优势比，若不写该语句，程序将自动输出变量b 各水平相对于基准水平的优势比。

SAS 程序运行结果：

Class level information

ClassValueDesign variables b110 201 300

以上结果为通过class 语句对变量b 赋哑变量的结果，以第三水平（即“未告知”）为基准，其赋值方式与写法b1 = 0; b2 = 0; if b = 1 then b1 = 1; else if b = 2 then b2 = 1;等价。

Summary of forward selection

StepEffect enteredDFNumber inScore Chi-SquarePr > ChiSq 1b2146.0932< .0001

以上为前进法（forward）筛选变量的结果，只有变量b入选方程（< 0.0001）。

Summary of backward elimination

StepEffect removedDFNumber inWald Chi-SquarePr > ChiSq 1a*b220.21330.8988 2a110.54780.4592

以上为后退法（backward）筛选变量的结果，因素a 以及a 和b 的交互作用，a*b 都被排除在方程外（值分别为0.4592 和0.8988），说明只有变量b 保留在方程中，该结果与前进法的筛选结果一致。

Summary of stepwise selection

StepEffectDFNumber inScore Chi-SquareWald Chi-SquarePr > ChiSq EnteredRemoved 1b 2146.0932 < .0001

以上为逐步筛选法（stepwise）筛选变量的结果，只有因素b 最终保留在方程中（< 0.0001），该结果与前进法、后退法的筛选结果一致。综合以上三种变量筛选法的结果，可以认为变量b 对结果的影响有统计学意义。

Type 3 analysis of effects

EffectDFWald Chi-SquarePr > ChiSq b242.6689< .0001

以上为对变量b 进行整体假设检验的结果，wald2= 42.6689，< 0.0001，说明变量b 对结果的影响有统计学意义。

Analysis of maximum likelihood estimates

Parameter DFEstimateStandarderrorWaldChi-SquarePr > ChiSq Intercept 10.74580.199613.95810.0002 b111.55210.255836.8130< .0001 b210.55350.22466.07220.0137

以上是关于截距项以及变量b 产生的两个哑变量的参数估计及假设检验的结果。截距项的参数估计值为0.7458，= 0.0002 < 0.05；哑变量b1的参数估计值为1.5521，< 0.0001；哑变量b2的参数估计值为0.5535，= 0.0137 < 0.05。据此可以写出该模型的表达式为：

Odds Ratio Estimates

EffectPoint estimate95% Wald confidence limits b 1 vs 34.7212.8607.795 b 2 vs 31.7391.1202.701

以上为程序默认的优势比的输出结果，给出了告知程度的1 水平（详细告知）和2 水平（告知）相对于3 水平（未告知）的优势比。1 水平（详细告知）相对于3 水平（未告知）的优势比OR= 4.721，其95% 的置信区间为（2.860，7.795），说明被详细告知患者的依从性是未被告知患者依从性的4.721 倍。2 水平（告知）相对于3 水平（未告知）的优势比OR= 1.739，其95% 置信区间为（1.120，2.701），说明被告知患者的依从性是未被告知患者依从性的1.739 倍。

Wald confidence interval for odds ratios

LabelEstimate95% Confidence limits b 1 vs 22.7141.8703.941 b 1 vs 34.7212.8607.795 b 2 vs 31.7391.1202.701

以上为语句oddsratio b 对应的输出结果。可以看到，除了给出变量b（告知程度）1 水平（详细告知）和2 水平（告知）相对于3 水平（未告知）的优势比外，还给出了1 水平（详细告知）相对于2 水平（告知）的优势比OR，比值为2.714，其95% 置信区间为（1.870，3.941），说明被详细告知患者的依从性是被告知患者依从性的2.714 倍。

专业结论：根据logistic 回归分析的结果，告知者对患者依从性没有影响，告知程度对患者的依从性有影响（< 0.0001）。OR= 2.714，其95% 置信区间为（0.352，0.567），说明被详细告知患者的依从性是被告知患者的2.714 倍；OR= 4.721，其95% 置信区间为（2.860，7.795），说明被详细告知患者的依从性是未被告知患者依从性的4.721 倍；OR= 1.739，其95% 置信区间为（1.120，2.701），说明被告知患者的依从性是未被告知患者依从性的1.739 倍。

[1] Hu LP. Medical statistics-analysis of quantitative and qualitative data applying the triple-type theory. Beijing: People’s Military Medical Press, 2009:325-351. (in Chinese)

胡良平. 医学统计学-运用三型理论分析定量与定性资料. 北京: 人民军医出版社, 2009:325-351.

[2] Hu LP. Medical statistics-analysis of regression models using the triple-type theory. Beijing: People’s Military Medical Press, 2010:203- 220. (in Chinese)

胡良平. 医学统计学-运用三型理论进行现代回归分析. 北京: 人民军医出版社, 2010:203-220.

10.3969/cmba.j.issn.1673-713X.2013.02.016

胡良平，Email：lphu812@sina.com