上QQ阅读APP看书,第一时间看更新
第二章 数据描述、变异建模与统计实践
杜鸿雁 1 谭铭 2 陈豪 3 金华 3
1美国NorthShore医疗健康系统临床信息研究中心
2美国乔治城大学生物统计、生物信息与生物数学系,Lombardi综合癌症中心
3华南师范大学数学科学学院概率统计系
第一节 概述
在分子生物学和遗传学快速发展的年代,医学研究在各级水平上(从基础学科到转化性研究乃至临床研究)以惊人的速度产生了大量数据,倘若这些数据不转换成信息与知识,其本身没有什么用处。作为科学工具的统计学的重要性就在于能够通过分析样本数据对未知的总体作出推断。变异的概念始终贯穿于统计学理论与实践中,它是我们的日常生活、我们的数据以及由数据所得到的统计估计中所固有的。由于人和人之间是不同的,即便是非常好的药物与治疗方法也可能仅对一些病人有效,而对另一些人无效。我们的血压值总是在变,无论你担心与否,或你是否处于良好健康状态,以及如何控制其他的未知因素(随机变异),它总是受时间、测量方法以及由谁测量等影响。随机变异是无法解释的变异。实际上,如何控制由于不同可能的因素而造成的变异正是统计实验设计的任务。由于随机变异永远都存在,实际工作中,我们不可能控制所有的影响因素。统计方法使我们能判断需要多少样本(受试者数)才能从趋势中分离出噪音。变异的统计估计可以使我们定量地确定生物医学研究结果中的不确定性。这些知识在实际工作中可以用来决定病人的治疗策略。这种努力的目标是基于遗传学信息来达到个体化用药。
很明显,数据的描述,变异的理解是同实验设计和产生数据的生物医学过程以及研究设计(为科学探索而产生数据的筹划过程)是紧密联系的。理解生物和医学的过程对于理解数据的意义非常重要。以简明的方式对变异进行描述和建模的方法与工具,使我们很容易地把数据转换为信息,其关键在于理解这些数据的变异。常用的分析方法在Moore和McCabe [1]1989年编写的书中第一讲以及Altman [2]1991年编写的书中已经介绍和总结了,这一章的目的在于更好地描述和理解变异,介绍几种比较先进的方法。因此本章节将包含基本的统计方法,以及一些最新的进展,比如离散混合效应模型。另外,我们为从研究设计到针对预测和过分拟合等统计分析实践,提供一些有用的一般性指导。并且,通过回顾分析一个以生活质量为结局变量的大型随机化试验,其中删失值是很大的挑战,我们给读者展示统计建模的过程,同时强调统计实践中需要特别小心。
第二节 描述数据的方法
对数的理解并对其作出任何推断的第一步是要了解所要处理的数据属何种类型,因为不同的数据类型需要用不同的统计方法来分析。实际工作中非统计学家有时会很容易忽略它。在本节中我们将首先回顾一些在生物医学研究中常见的数据类型并特别强调不常见于教科书但却在医学研究中越来越常见的某些类型,我们还将指出应当用来分析它们的方法。
一、数据类型与测量
尽管生物医学研究中的数据常常是复杂的,但它们确实可以被分为几种常见的类型,对它们的理解将会指导我们选择正确的方法来总结和分析它们。数据的类型决定了用什么样的方法来分析数据作出推断。下面除了回顾一些基本的数据类型外(Altman [1]),我们还将叙述其他的一些在现代生物医学研究中越来越多见的数据类型。
1.分类数据
当一个病人的病情或他的某些状况能够分成几种类别时,对他们的观测可以产生分类数据(有时也称为二分类或定性数据)。最简单的例子是yes或no两个分类的观测,如一个病人对肿瘤治疗是否有反应、吸烟还是不吸烟,或这个病人是否患有结肠癌等。这种类型的数据有时被称为二分类或0-1变量数据。三分类或更多分类的数据有如血型(A、B、AB、O);复合分类如男性或女性白血病和非白血病病人。由于这些分类型没有明显的顺序(血型、性别、疾病分类),这种类型的数据也称为名义数据。它们可以用列联表或广义线性模型的方法来分析。另一种分类数据如吸烟的分类(完全不吸、有时吸、严重吸烟),乳腺癌的分期(Ⅰ、Ⅱ、Ⅲ、Ⅳ),治疗后改善的程度(无改善、中度改善、明显改善、完全改善),以及作为主观评价疼痛的等级(轻微、中等、严重、无法忍受)。像这样在所有类别之间有很明显的等级时成为有序数据。然而就像名义数据那样,尽管某些有序数据也会以数字形式出现如乳腺癌的分期(Ⅰ、Ⅱ、Ⅲ、Ⅳ),但有序数字的算术运算没有意义,如我们很难说无法忍受的疼痛是严重疼痛的两倍。
有一种不同的情况是用分数来表示一定结果的计分数据,这种数据中从一个点到另一个较高点的增量确实是等距的。
2.连续性数据
连续性数据通常是指通过测量得到的数据,如体重、气温、血压和大多数血液化学检验的数据(如胆红素、血红素、胆固醇)。这些观测或它的转换值(如对数转换值、平方根转换值)通常认为是服从正态分布的,目前连续性数据的统计方法和分析模型是发展最全面的,然而这些测量的准确性和可靠性对于作出有效的推断非常重要(Gleser [3])。尤其应当注意的是当在分析中这些观测被用作自变量时,有时需采用变量误差(或测量误差)模型(Carroll等 [4])。一般教科书通常没有强调测量误差的影响及其重要性。最近的研究将测量误差模型进一步推广到广义线性模型和生存分析方面。连续性数据也包括区间标度数据(如一列温度10、20、30)和比率标度数据(如年龄10、20、30)。诚如其名,将年龄上的30岁解释为10岁的3倍之老(30/10=3),这样是说得过去的;但将温度上的30度说成是10度的3倍之热,这样的说法是没有意义的,因为温度没有绝对零度,只能说30度比10度热20度(30-10=20)。
3.比
当我们取两个变量的比值时就产生了比数据。例如排除率(一个重要的心脏功能指数),是指心脏舒张末期容量的差值同心脏收缩末期容量之比。还有与一定基值相比的肾脏功能的百分比变化(如肾小球滤过率)。最近,微点阵基因表达比成为许多边缘医学研究的焦点。微点阵技术可以允许大范围(上千个基因)快速的基因表达分析。在这些实验中,每个点即每个基因的表达情况是以疾病组织样本的基因表达参照样本的基因表达之比来表示的。对这些基因表达比的分析必须考虑到比的来源需要利用适当的分布(实际上多为Gamma分布)(Chen等 [5],Newton等 [6])。
4.连续性比例数据
它实际上是当比为0到1的百分数时的一种比值数据的子类型,例如,与基线相比不同时间肾功能降低的百分比,关于某些生理变量及分子或基因靶体在治疗前后变化的百分比。用Barndorff-Nielsen和Jorgesen的单纯形分布(Jorgesen [7])直接针对比例反应的均值建模的统计方法才刚刚出现(Song和Tan [8],Tan [9])。单纯形分布考虑到这类反应变量是0与1之间的百分比,而且有很大的离散型。
5.重复测量
在一些疾病的自然研究或疗效观察中,受试者常在不同时间被重复地随访或观察,或对一些实验单位(如双眼或两侧肢体)进行多项测量或观测。这类研究中所获得的数据称为重复测量数据(repeated measures)。当这些数据是对同一个体不同时间的观察值,也称为纵向数据(longitudinal data)。这种实验设计对于评价疾病长时间的变化情况常常是必要的。例如,我们可能会对一些生理变量(如胆红素滤过率)或遗传变量(如端粒的长度)的长时间的变化情况进行观察,或观测长时间内某些事件是否发生(如耳朵的感染)。
在这类事件设计和数据分析中,关键是要考虑病人自身和器官群间的相关性。例如,某个儿童在一侧耳朵感染时,另一个耳朵也容易感染。因此,10个病人(每人测量2只耳朵)的功效不如20个病人各测量1只耳朵的功效。
根据反应变量的类型,我们可能有连续性数据或分类数据或等级数据。不同类型的重复测量数据需用不同的统计模型来分析,尽管目前都可以统一用广义线性模型(Diggle,Liang和Zeger [10])来分析。但对于重复等级资料,还可使用广义线性模型以外的模型,例如,比例优势模型(Qu和Tan [11],Tan等 [12])。
6.删失数据
当我们不能精确地测量一个观测,而仅仅知道这个观测是超过一定的阈值时,就称这个观测被删失。医学研究中最常见的删失数据是生存数据,就是到一定事件发生时的时间长度数据,这事件可以是病毒感染,或病人的死亡。这或许也是医学研究中最常见的数据类型。因为在医学实践中我们经常想知道某种新药或某种外科手术或某个医学手段是否比常规方法挽救更多生命,或更多的延长生存时间,由于以下几个原因,生存数据的分析常需要特殊的方法。首先,生存数据的分布不是对称的,因此也不服从正态分布,在建模时用其他的分布会更加满意;第二,在分析时间数据时常有某些病人的生存终点还没有被观测到,甚至某些病人的生存状态由于失访而无法知道。
二、变异
变异是所有统计理论和方法中最重要的一个基本概念。世界充满了不确定性,有幸有统计方法来科学地研究不确定性的存在。由于实验中,尤其是人体试验中变异的存在,通常一个具有生物活性的药物只有5%的机会进入临床。变异使医学研究中统计学和统计学家都必不可少。所以与此有关的另一个基本概念就是用来描述和分析数据的概率分布。我们经常假定观测值来自只有一些参数未知的某种分布,最重要的分布是正态分布,它在统计学中具有非常重要的地位,因为中心极限定理表明许多常见的统计量都近似服从正态分布。不假定参数形式的方法称为非参数方法。参数模型的优点是简单有效。有时也采用一种折中(半参数)的方法,只是将感兴趣的主要特征假定为一个参数模型。目前,参数和半参数方法是医学研究中最常用的方法。
三、基本方法
最常用的数据描述指标是均值与标准差,它常同数据的参数描述联系在一起。正态分布完全由它的均值与标准差确定。均值用于测量中心位置而标准差描述变异性。由于正态分布在统计推断理论中的重要性,这两个指标具有特殊的意义。但如果研究变量的分布不是正态的,它们则不一定能给出好的推断值。有时变异会超出所假定分布能够描述的范围(即所谓的过度离散)。
另一些常用于描述数据的统计量是5个数值综合统计量,即最小值、最大值、75%、50%(中位数)和25%(分位数)。同均值和标准差一起,这5个数值能很好概括数据的分布。例如,如果分布是对称的,则均值与中位数应该相等。若均值比中位数大,分布呈负偏峰,若均值比中位数小,分布呈正偏峰。
例2-1 一个关于topotecan在实体癌中药代动力学研究的Ⅰ期临床试验。Topotecan是一种新的分子靶向性抗肿瘤药物。这项研究旨在确定Ⅰ期临床试验中是否可通过药代动力学指导下的剂量调整措施而降低Topotecan内酯系统暴露的变异,暴露情况是用血浆浓度-时间曲线下的面积来衡量的。试验是给15个有复发实体癌的儿童在连续两周内、每周5天,静脉注射topotecan,通过调整药物的剂量以使每天topotecan内酯血浆浓度-时间曲线下面积(AUC)尽量保持在某个特定的目标范围内,队列1的目标为(150±30)ng/ml*hr(9个病人)。其中队列1的两个病人由于毒性过大而移至队列2(10天后),即把AUC的目标范围由(150±30)ng/ml*hr降至(150±20)ng/ml*hr。试验期间,分别在第1、3、6、8、10天测量AUC。
四、作图法
事实上有时一幅图胜过数千个文字。在统计学和医学研究中经常使用作图的方法来描述和说明。
例2-2 前述的5个指标就常常被集中反映在箱须图(box-and-whisker)中(见Altman [2]),中线代表中位数,箱子的上下两个边缘代表25和75百分位数(或者更低到更高的百分位数),两个须分别指的是最小和最大值。为了观察15个病人AUC数据的总体分布和变异,图2-1给出了在第1天(第一次给药)和第3天(第二次给药)15个病人AUC的箱须图。如图可见第1天AUC的分布偏向右侧,很不对称;而第3天的分布则比较对称,这部分显示了基于药代动力学的药物剂量调整的效果。
图2-1 第1天与第3天AUC的比较
为了描述数据的分布和变异,常常使用直方图和一些平滑技术。样条平滑估计是对分布密度的一种非参数估计,它可以较好地描述分布的概率密度。用现代的统计软件是很容易产生这种估计并将其重叠在直方图上。图2-2给出了例2中队列1的8个病人的AUC直方图及拟合的密度曲线。其中“固定剂量组(fixed)”指的是8个病人的36个AUC值,其AUC值是以一个固定剂量(4mg/m 2)除以每个病人topotecan内酯清除指标来计算的;“剂量调整组(targeted)”指的是同样这8个病人的27个实际AUC值,而第一次给药的8个AUC值及第二次给药的1个AUC值被排除了。
图2-2 固定剂量组(Fixed)与剂量调整组(Targeted)AUC的比较
第三节 通过模型调整因素来描述数据
就像前面提到的,数据的有效描述依赖于设计。有时对复杂设计下产生数据的总结不能够直接得到,比如,当观测之间存在依赖性或存在缺失数据时,一个直接由原始数据得到的均值和标准差可能会误导。如在例2-1中病人经历了多次化疗和药代动力学研究,而某些病人缺失了若干次,这样就产生了非平衡重复测量的数据结构。我们可以使用混合效应模型来估计达到目标(AUC值达到控制目标范围内)和没有达到目标的病人的药代动力学参数。对于这类数据,综合统计量(如均值与标准差)的计算需要考虑数据间的相关。为了便于比较,表2-1列出了考虑相关(调整后)及未考虑相关(未调整)时取得估计值。由表中可以看出,基于模型,利用所有数据所得的综合统计量的数值与简单计算得到的统计量数值不同,如作统计推断应使用考虑相关的统计量。
表2-1 基于混合效应模型及简单计算的均值与标准差
为了避免偏性,需要作基于模型的精心估计。从Meyers,Nelson和Tan [13]以及Nelson [14]还可看到其他的一些例子。例如利用混合效应样条模型来估计糖尿病在生活不同阶段的平均肾小球滤过率及标准误。我们将在4.4节与第5节中分别对连续性比例数据和不完整纵向数据探讨基于模型的方法。
第四节 过度离散问题
过度离散指的是观测到的变异(方差)大于某个假定模型下名义上变异的现象。从统计学角度上说,过度离散表明所假设分布的均值-方差关系是不正确的。尽管通常认为过度离散常发生在二项分布或Poisson分布假定下的离散数据模型中,近年来在Song和Tan [8]的研究中发现它也可以发生在连续性比例数据中。实际上在统计学中,很早以前就注意到了过度离散的存在。1951年Fisher就注意到在实践中许多数据是过度离散的。这就自然对我们提出了问题,分析中忽视过度离散的存在将会有什么后果?有哪些适当的方法可以用来对过度离散的数据进行检测和建模。在这一节中,我们将在几种分布中讨论这些问题,包括已熟悉的二项分布和Poisson分布数据,以及关于比例数据的新进展。
一、过度离散的二项分布数据
两分类结果,如治疗的成功/失败,对肿瘤药物有无反应等,是医学研究中最常见的。通常我们定义成功的概率为 p, n分类序列(例如, n个细胞, n个小鼠或 n个病人)中的每一个二分类(0-1)结果为 Y i,则二项分布的结果可以表示为 。当经验方差大于二项分布的方差 np(1- p)时,数据就存在过度离散。方差是均值 p的一个函数,这时分布是完全由参数 p所决定的。二项分布的方差为
因此,当二分类数列彼此间不独立时,也就是说一些协方差cov( Y i, Y j)不为零时,就会产生过度离散现象。其后果取决于过度离散程度的轻重。一般说来,过度离散是不能忽略的。
过度离散可以用广义线性模型来检验。随着近年来发展起来的广义线性混合效应模型和Bayesian等级模型,可以直接通过建模来考虑过度离散。
二、过度离散的Poisson分布数据
类似于二项分布数据,Poisson分布也是由它的均值参数来决定的。由于Poisson模型属于广义线性模型,与二项分布中类似的检验统计量和建模方法可以用来对Poisson分布数据中的过度离散进行检验。
三、过度离散的连续性比例数据
连续性比例数据和方向性数据在以前的文章中很少被谈到。当反应变量是位于0到1的百分数时就会产生连续性比例数据,例如同基线相比不同时间的肾脏功能降低百分数,或者同基线相比血压降低的百分数。实际工作中常把它们当成正态分布。然而,就如Song和Tan [8]指出的,反应百分比的变异超出了正态分布能够描述的范围。尽管当离散参数很小时,离散模型近似于正态分布(Jorgensen [7]),然而真实世界的数据却经常存在很大的离散。这时用正态分布来描述和分析是不合适的,因为当两个变量服从正态分布时(这个假定常被认为是合理的),而两个变量的比值一般却并不服从正态分布。
例2-3 本例是一个研究视网膜修复中眼内气体作用的前瞻性眼科学研究(Meyers等 [13])。结果变量是眼内气体残留百分比。31个病人在手术前向眼内注射了气体。在随后的三个月期间随访了3到8次(平均5次)。每一次随访眼内气体容量都以其占初始气体容量的百分比来记录。研究问题是要估计气体消失的动力学(如气体衰减率)。反应变量定义在0到1之间。尽管logit变换后的反应变量可以使用线性回归模型,但是非线性转换后的反应变量常常很难解释。尤其,非线性转换后反应变量的系列相关结构常难以再转回到原来的反应变量。研究目的是为气体衰减均值对若干协变量的依赖直接建模。通常假定反应变量是正态分布的而忽视反应变量是定义在0到1间的百分比,然而经过以下的分析我们会发现百分比反应变量的变异超出了正态分布能够描述的范围。
离散参数 σ 2的矩估计可以通过 d( Y; μ)的期望值等于 σ 2来估计
若 ij是一致性估计,当 m趋于无穷大时,它就是 σ 2的一致估计。
在例3中,离散参数 σ 2的估计等于14.2。基于自由度为2的 χ 2分布 p值为0.0008,提示离散参数明显大于0,也就是,明显大于正态分布的离散程度。提示气体容量百分比不服从正态分布。实际上,在图中也可以看出离散参数很大时单纯形密度函数大多集中在0.8到1之间,这一点同数据的分布特性相符,也就是说多于40%的观测值位于这个范围内。另外,Song,Qiu和Tan [15]考察了离散程度随时间和其他协变量是如何变化的,并将边际模型推广到允许离散程度随时间而改变的情形。因此,在这种数据的分析中应当考虑过度离散。
四、连续性比例数据的离散性建模
对于上面的例子,已经知道气体容量并不服从正态分布,存在过度离散的现象,因此下个问题就是如何分析这种过度离散的纵向数据。因为一个重要的目标是对特定个体进行推断,所以一个自然的模型是混合效应模型,其中结局变量服从单纯形分布。下面简要介绍Qiu,Song和Tan [16]提出的这种单纯形混合效应模型(SMM)。
令 y ij是第 i个个体在时间 j的反应变量,而 x ij与 z ij分别为对应固定效应和随机效应的解释变量。假设对于第 i个个体,给定其 q维随机效应向量 b i, y ij是条件独立的,均值为 μ ij= E( y ij| b i)。因为结局变量服从单纯形分布(见附录),故有 y ij| b i~ S -( μ ij, σ 2),而条件均值 μ ij与预测变量(固定效应和随机效应)相关关系为:
其中 g为连结函数, β为 p维固定效应向量,随机效应 b 1, b 2,…, ,且 D= D( θ)依赖于未知参数 θ。对于二分类和连续性比例数据最常用的连接函数为logit连接,它可表示为 g( μ)= log{ μ/(1- μ)}。记 ,且 对应第 i个个体。联合密度函数为:
其中 p( y ij| b i)为单纯形分布的密度, p( b i)为正态分布的密度。那么参数 β, σ 2和 θ的对数似然函数正比于
因此,积分掉随机效应就得到 β与 θ的边际对数似然函数:
其中
从增广似然函数积分掉随机效应得到的边际对数似然函数往往没有明确的解析表达式。为了对模型进行估计,Qiu,Song和Tan [16]针对广义线性混合效应模型的近似拟似然函数方法(PQL/REML)做了一些推广。这种方法在概念上是简单的,数值上是稳定的,并且适用于任何维数的随机效应。尽管这种近似推断参数估计可能有偏差,但其算法速度常常很快。Qiu,Song和Tan [16]证明了这种估计偏差在使用四阶拉普拉斯逼近时能够减小到满意的水平。
例2-4 为分析气体冷却数据,我们利用下面单纯形混合效应模型(SMM):
logit( μ ij)= β 0+ b 0 i+ β 1log( t ij)+ β 2log 2( t ij)+ β 3 x ij,
其中随机截距 b 0 i~ N(0, θ 0)。另外,我们也考虑包括随机斜率 b 1 i~ N(0, θ 1)的如下模型:
logit( μ ij)= β 0+ b 0 i+ β 1log( t ij)+ β 2log 2( t ij)+( β 3+ b 1 i) x ij
其中 x ij为标准化气体浓度协变量(取值为1、0或-1), t ij为时间协变量(手术后的天数)。我们利用校正的PQL方法得到 θ 0=0.26(0.19), θ 1=0.09(0.25)。结果表明,两个方差参数无统计学意义,这意味着刻画序列相关的参数可能过多,因此我们考虑简单的随机截距模型进行统计推断,细节可参见Qiu,Song和Tan [16]。PQLs分析发现:气体浓度水平 x ij和时间对数的平方log 2 t ij是有统计学意义的( p﹤0.05)。然而,如果只是对气体容量做logit变换,再使用常规的线性混合效应模型分析,那么气体的浓度并无统计学意义( p=0.14)。
第五节 统计实践
这一节我们将探讨统计实践,特别是统计应用于生物医学研究中的问题,并讨论为什么需要一个适当的研究设计、严谨的分析和合理的推断。统计应用于生物医学研究常常遵循三个步骤:首先确定研究设计的性质,明确研究的问题;然后确定数据的类型(变量的测量尺度)和要比较的组数,从而形成一个完整的统计分析方案。
一、研究设计
一项研究最基本要做到设计优良,因为没有什么统计方法能使设计不当的研究起死回生。一个研究可以是实验性的,也可以是观察性的。在实验研究中,研究者给研究对象(如病人,动物或临床观测点)分配治疗或干预,以确定该分配或治疗是否有效。随机化常被使用,其目的是得到对病人所有特征(如年龄或性别,是否接受治疗这点除外)具有可比性的分组。如果感兴趣的结局变量组间存在差异,那么研究者可以有信心认为这个治疗或干预是有效的。另外,具有统计学意义并非意味有临床意义(实际意义)。比如,对于糖尿病患者而言,改善22mg血糖或许没什么意义。因此,在研究设计时确定一个有实际意义的差异是很重要的。举个例子,临床试验设计必须要有足够的功效(如80%)能够检测预先确定的具有实际意义的差异,同时样本容量不能超过必要样本容量。另一方面,在观察研究中治疗或干预的分配不受研究者控制。无论是实验研究还是观察研究都有局限性,这可由严格分层设计(Concato [17])解释。一项研究也可设计成准实验设计,此时实验者可能很少或者不能控制治疗的分配。我们除了要确定研究是实验性的还是观察性的之外,还要弄清数据收集是回顾性的还是前瞻性的,以及数据是在一个时点(截面的)还是多个时点(纵向的)获得的,因为需要特定的统计方法刻画相关和协方差结构,用于分析对同一个体重复测量而产生的相关观测结果。
二、研究的问题
研究的问题要有重点、明确并且具体。在许多情况下,它很容易转化成两个统计假设——原假设与备择假设。原假设是一个陈述,假定它是对的,直到出现充分的相反证据(即治疗无效),而备择假设代表研究者的兴趣(即治疗有效)。重要的是,在进行探索性数据分析与统计假设检验之前,都需要确定一个检验水平(α)。通常 p﹤0.05或 p﹤0.01被认为具有统计学意义, p﹤0.1被认为是有沾边的统计学意义。一般来说,研究报告的置信度是与统计检验的水平连在一起的。研究者的兴趣或要研究的问题是明确的,比如,刻画变量间的相关性,评估不同实验方法的测量值的一致性,检测某一时点以及随时间变化时不同组之间平均值的差异。在进行统计检验前明确研究假设是很重要的,这样按照设计好的方案,研究人员就可以避免危及研究结论有效性的陷阱。同时,在某些情况下,没有统计学意义与研究者特意寻找的统计学意义同样重要。为了行之有效,统计学家需要与临床医生和医学工作者一道紧密合作,以保证大家对分析结果达成共识。因此,除了严格的统计训练,在生物医学、合作研究经验和沟通技巧的跨学科培训也很重要。
随着近来遗传学与人口学数据的爆炸性增长,对研究问题缺乏清晰的理解对统计推断造成了大量的混乱状况。牢记研究的目的并运用恰当的统计模型是非常重要的。目前一个值得关注的问题是,研究的目的是探索性的还是为了预测?这在统计技术工具箱日益丰富的今天常被忽略。下面章节将详细讨论这个问题。
三、研究目标是探索性的还是预测性的
就本身的术语而言,探索性建模常用来描述一个变量和我们要寻找的与之相关的一些解释变量之间的关系;只有因果因子发生在结局变量之前,才能建立一个因果关系的解释。另一方面,预测性的建模是基于已有数据,得到统计模型,以预测新病人的未来结果(响应)。因此,探索性统计模型的功效不同于预测性统计模型的功效。
一般而言,统计建模是要寻找最简单的模型,变量最少,数值上要稳定,还要容易推广。临床工作者在流行病学研究中趋向于包括尽可能多的混杂因素,即使这些因素没有统计学意义,这会使得模型更依赖于给定的数据。这样建模的优点是模型常常有好的拟合性,然而这会导致“过度拟合”,参数增多,标准误增大,因为模型中变量越多,变异就越大。所以,如此构建的模型可能有很好的探索功效但是缺少预测能力。在医学研究中,特别是在评估用于早期癌症检测的生物标志的诊断医学中,重要的是构建一个预测法则,而不是描述或解释数据间的关联。
常用在流行病学研究中的(线性和非线性)回归分析,经常将 R 2(决定系数)作为模型可解释的方差百分比,它可度量模型的解释能力。另一方面,相同的模型可用于其他目的。比如,考虑是否有某种疾病(如结肠癌)的logistic回归概率模型,常用于诊断医学(带有预测目的)中的疾病分类。接受者特征曲线下的面积(AUC),也被称为(Wilcoxon-Mann-Whitney检验中的)c-统计量或一致性概率,是模型预测性判别功效的一个好的概括性度量。为得到评判模型预测性能的无偏估计常要用到重抽样技术。
探索性建模与预测性建模的本质区别很微妙但是很重要。最近,Shmueli [18]在这方面做过深入的探讨,提出在构建模型与检验模型上两者都很必要,只是分别扮演不同的角色而已,正如与解释相关的不确定性本质上不同于与预测有关的不确定性(Helmer and Rescher [19])。例如,流行病学研究的趋势是通过调整所有潜在的危险因素去评估关联性,而另外一些人则致力于寻找生物标志,希望能用于预测病人未来的结局。统计上区别探索性建模与预测性建模的是偏差-方差的折中。正如Hastie,Tibshirani和Friedman [20](223页)所示,预测误差估计(EPE)可分解为噪音、偏差和估计方差,即
EPE=噪音+偏差+估计方差
其中噪音是随机误差(无论模型有多正确、估计有多精确),偏差是系统误差,刻画模型偏离给定数据的程度,而估计方差是用样本估计模型引起的,即是真实模型拟合相似数据集的抽样误差。用一个简单的比喻来解释这个概念,假设有个并不标准的温度计,那么它测量人的体温常常(系统地)会高1度,这个系统误差就是噪音。再假设无论某人发高热或体温正常,都有一个温度计(真实模型)能准确测量体温:一个温度计(模型A)对于体温较低的人表现良好(偏差较小),另一个温度计(模型B)对于体温较高的人表现较好(偏差较小)。尽管模型A与B的偏差都小,但是其中只有一个用于预测,因为事先并不知道人的体温是多少。模型A或B与真实模型之间的差异就是估计方差。上面的等式揭示了探索能力与预测能力的不同之处,其中探索能力用偏差表示,可由给定的数据得到,而预测能力用估计方差表示,需要未来的数据评估。如果预测误差估计(EPE)固定且噪音得到控制,那么减少偏差就会增加估计方差,反之亦然。因此,从探索目的得到的模型可能预测效果不好,为了预测,避免过度拟合很关键。在分析大量基因数据时,观测的数据量常比预测变量要少很多,因此过度拟合容易发生。受此问题驱动,出现了很多模型与方法。例如,典则回归或统计学习模型(Hastie等 [20])已经成为分析这类问题的有效工具。最近,在高维数据情形,也有一些模型直接最大化接受者工作特征曲线(ROC)的效用函数(如ROC曲线下的面积), F-度量值,以及灵敏度与特异度的线性组合(见Liu等 [21,22])。
四、数据分布,正态性假设与稳健性
一旦明确了研究设计与要研究的问题,并且将它准确地转化为统计假设,我们就要处理数据了。如前所述,需要评估数据的分布以确保要进行的统计检验的假设条件完全满足。如果检验需要的假设条件不满足,可以选择别的检验方法。例如,参数统计方法,如 t检验,方差分析(ANOVA)以及Pearson相关系数,都需要满足下面三个假设:正态性、独立性与方差齐性;然而非参数方法,如Wilcoxon秩和检验,Krushal-Wallis检验以及Spearman相关系数,都与分布无关,并且不需要严格的假设条件。一般而言,参数方法更加敏感,功效更高,但需要较大的样本量,而非参数方法较不敏感,功效较低。
如果参数方法涉及到正态反应变量,那么常常需要评估模型的正态性假设。一旦不满足正态性,则需要进行一些适当的变量变换,如对数变换,平方根变换或逆变换。也可以对反应变量 Y进行Box-Cox变换(Box and Cox [23]): T( Y)=( Y λ-1)/ λ,其中 λ为最优功效系数。自然对数变换是它的特例: λ=0。
例2-5 本例是一个随机化纵向临床试验,结局变量是与健康有关的生存质量,如最新诊断出的患有慢性粒细胞白血病患者在12个月内的、由癌症治疗生物反应调节剂功能评估(Bacik J等 [24])度量的试验结果指数(TOI)。病人在试验开始被随机分配到两个治疗组中。TOI得分介于0(最差)到108(最好)之间,每个病人测量多达九次:包括试验开始以及1、2、3、4、5、6、9和12月后。研究的基本目的是确定两组TOI得分均值是否随着时间的推移有差异。主要协变量包括治疗组(Ⅰ和Ⅱ),年龄(均值为中心)和性别。在开始接受随机分配的1049位病人中,共有979位确认有初始TOI得分,他们包含在Du等 [25]的分析报告中。TOI的分布呈正偏,需要用平方根转换进行校正,因为这种情况对数转换不够充分,而逆变换则有点矫枉过正。另外,为要假设仅在开始观测到的总体分布与那些后面观测的(6月或9月后)相同,需要在6月与9月后评估病人特征(性别和年龄)的相似性,而这些特征都进行了矫正,无论模型有没有包含结局变量TOI。鉴于删失数据超过10%,需要评估删失对于健康相关的生存质量的估计的影响。已有联合模型来分析具有正态假设随机效应的纵向TOI(混合效应)和中途退出时间。
尽管混合效应模型被广泛用于分析纵向和相关数据,其中随机效应常假设为高斯分布,检验这种联合模型的基本假设及其稳健性是很重要的。一些作者已经考虑了线性混合效应模型中的有关问题(Zhang and Davidian [26])。最近,McCulloch和Neuhaus [27]基于实例、理论计算以及随机模拟考虑了用于聚类或纵向数据分析的广义线性混合效应模型广泛应用中的稳健性,包括预测、协变量效应、随机效应的预测以及随机效应方差的估计。他们认为:对随机效应分布错误指定的担忧通常是瞎操心,因为①灵敏度仅限于估计方面,而这往往不是感兴趣的;②考虑的情况过于极端;③发表的结果实际上并不支持错误指定的敏感性。
Du等 [25]用一个更加复杂的联合模型来考虑事件发生用的时间和生存质量纵向结局数据,并研究了随机效应方面的稳健性。我们利用关于生存和纵向结局的联合模型理论与方法、通过纵向测量次数的分布来评估纵向测量有多密集(Hsieh等 [28];Huang,Stefanski和Davidian [29])。给定数据TOIs以及固定效应参数值,随机效应的后验分布在每个个体纵向测量次数适当大时近似正态。由于随机效应结构和纵向测量(TOIs)都只通过此后验密度与生存参数相关,当纵向测量比较密集时,最大似然估计(MLE)就会稳健,即使关于生存时间(或中途退出时间)和纵向数据的联合模型中的随机效应偏离正态性假设;另一方面,如果纵向测量比较稀疏(如,每个个体只有3个观察值),那么最大似然估计可能会偏差很大(Hsie等 [28])。表2-2显示如何通过纵向测量次数的分布来考察纵向测量的密集程度。
表2-2 纵向测量次数的分布
由表2-2可知,超过90%的病人至少观测4次(或以上);中位数大于8次。分组考虑,大约第一组56.5%和第二组39.2%的人员所有九次测量都有TOI得分。因此,就随机效应的正态性假设而言,纵向TOI测量次数确实不少,最大似然估计是相当稳健的。
对于二分类结局进行logistic回归,连续性协变量需要满足logit线性假设(Hosmer和Lemshow [30])。Box-Tidwell变换可用来检验logit中的线性性,即只要在模型中添加x*log(x)这样一个变量,然后考察其参数估计的统计学意义:这项参数估计有意义就表明非线性。
五、选择正确的统计检验
在统计实践中,针对涉及的问题选择正确的统计检验显然至关重要。在这一节,我们进行一些总结并给出实例去说明该如何根据研究设计与问题选择正确的统计检验;同时特别关注不同教科书中零乱、有时还被忽略的一些方法。
1.评估关联
两个变量的相关性:
当去评估对同一个体进行观测得到的两个变量(如成像参数和生物标志物)的关系时,常用的方法是将数据呈现在散点图上以确定这两变量间的相关性,其中解释变量放在x轴,而反应变量放在y轴上。散点图是确定变量间整体关系的有用工具,一旦线性关系被确定,研究者就可以进一步计算相关系数,度量相关强度及方向,并作出拟合数据最优的最小二乘回归直线。
两次观测的一致性:
Bland和Altman [31]指出,高度相关并非意味着好的一致性。Bland-Altman给出了对同一个体的两次观测(或两个方法、两个评分者)一致性的评估方法:作一个散点图,横坐标是两次观测的平均值,纵坐标是两次观测的差值,然后画出一致性的范围,其中界限由平均差值±1.96倍的差值的标准差。对连续型变量,可用组内相关系数(ICC)定量刻画一致性,而对分类变量,常用Kappa统计量去度量。
疾病与暴露的关系:
在控制混杂因素去评估暴露与疾病的关系时,常用Mantel-Hantzel检验,当然也有一些更高级的方法(Kosinski和Flanders [32])。
2.检测组间差异
横截面观测:
对于两组的比较,在正态性、独立性、方差齐性的假设条件下可用 t检验,否则要用Wilcoxon秩和检验;对于三组或者更多组的比较,在上述提到的假设下可用方差分析(ANOVA),否则用Kruskall-Wallis检验。
纵向或重复观测:
在纵向研究中,每个个体都是重复测量的,因此独立性假设(观测之间相互独立,单个个体的观测值不受其他个体的影响)不满足,需要利用特定的方差-协方差结构去解释个体内的相关性。另外,组间差异不仅表现为每个时点上的差异,而且也表现为重复观测下(或纵向观测下随时间的)结局变量的整个变化。对于匹配数据(即每个个体测量两次),连续性变量在满足正态性假设的前提下可使用配对 t检验,二分类变量(即是与否)可利用McNemar检验。当多于两个观测时,可使用重复测量的方差分析(ANOVA),或基于针对特定个体推断的似然估计的广义混合效应模型,或基于在总体水平推断的广义估计方程(GEE)的总体平均模型(Diggle等 [10])。特别地,对于涉及双胞胎(同卵与异卵)的研究,除GEE方法外还可使用经典的结构方程模型,例如关于压力性尿失禁与非遗传因素关系的流行病学研究(Nguyen et al [33])。在经典的双胞胎模型中,同卵与异卵双胞胎在表现型方差的差异可以归结于三个基本因素:可加的遗传性影响,公共环境影响以及特定个体或独特的环境因素,假设为:同卵双胞胎有相同的可加遗传变异,而异卵双胞胎仅仅共享一半的可加遗传变异,并且两组双胞胎都有相同的公共环境变异(Neale等 [34,35])。
六、建立模型
建模的第一步需要进行全面、仔细的单因素分析,以识别任何感兴趣的模式、关系和关联,然后建立多元主效应模型:刚开始包含所有单因素分析中 p值小于0.25的协变量,然后利用常用的逐步法、向后法或向前法,以0.05为阈值对候选变量进行筛选和剔除。多元模型中每个协变量回归系数的大小需要与单因素分析的结果对比,来决定重要的变量是否不应该剔除。一旦多元主效应模型确定,就要检验有临床意义的交互项和评估潜在的混杂因素。注意混杂因素仅仅能被评估而不能被检验,且应该与结局变量和危险因素相关联;一个协变量常常可作为混杂因素,如果校正它使得主要危险因素的影响幅度改变10%甚至更多,例如在校正年龄之后,吸烟对出生体重轻的回归系数从0.5变为0.56(变化﹥10%),因此,年龄被认为是一个混杂因素,应该包含在最终的多元模型中,即使某些情况下年龄可能不具有统计学意义(即 p﹥0.05)。最后,在做统计推断之前,需要评估模型的充分性和拟合程度。
高维遗传数据的观测次数通常远远大于变量/特征数,为分析这些高维遗传数据,已有大量文献使用统计学习方法,如LASSO特征选择(Hastie等 [20])。这些方法同样适用于经典的数据分析模型,因此可以使用基于LASSO的变量选择方法代替向前或向后筛选法。
此外,模型的建立通常是连贯的。在考查了数据及其分布之后,就需要更加复杂的模型进行适当的数据分析,同时重要的是用一种容易理解的方式(如图示)将结果向更广大的读者有效地展现出来。为阐明这一点,我们在例6中继续分析Du等 [25]的模型建立和数据分析过程。
例2-6 续例5。建立的第一个模型是通常的随机缺失(MAR)机制下的混合效应模型: ,其中第 i个个体在时间 j处测量, i=1,…, n且 j=1,…, 表示固定效应, W 1 i( s ij)表示个体的随机效应,且 。此外,建立一个混合模型(pattern-mixture model),其中个体根据缺失的模式分为若干组,而缺失模式作为一个个体间的变量被包含在纵向模型中。此模型与第一个模型一样,只是包括了中途退出(dropout)状态这个协变量。因为无12月的TOI的模式是稀疏的,我们把所有这些模式都当作是中途退出的,中间断断续续缺失的观测都作为随机缺失,并将那些可在最后12月得到TOI观测值的定义为完成者(Hedeker和Gibbons [36])。因为中途退出时间(time-to-dropout)可能很重要,所以我们构建一个参数分享联合模型,其中纵向TOI与中途退出时间过程同时建模。对中途退出时间建立指数回归模型,其中 t(月)时刻的危险率为: ( t))。这里 同样表示固定效应, W 2 i( s ij)表示个体随机效应。这两个模型的连结函数是 W 1 i( s)= U 0 i+ U 1 i s和 W 2 i( t)= γ 0 U 0 i+ γ 1 U 1 i,其中 γ 0与 γ 1分别为刻画这两个模型关联度的随机截距与斜率。为简单起见,不考虑截距与斜率之间的相关性,因为其大小可以忽略( r﹤0.013)。从研究开始日期到跟踪完成者12月之间的时间定义为缺失时间。对于中途退出者,从研究开始日期到最后一次HRQL收集时间与计划最后一次访问时间的中点之间的时间定义为中途退出时间。协变量包括年龄、性别、时间的四次多项式,而随机截距与斜率的方差-协方差并非结构化。在联合模型中,中途退出时间假设服从指数分布。联合模型揭示:除了研究开始外的每次访问,组间差异都有统计学意义( p﹤0.001)。参数估计(详见Du等 [25])类似于单独的纵向与生存子模型:两者间的关联参数有统计学意义( p=0.039),意味着TOI的斜率与中途退出的危险率之间呈负相关,因此不能忽略中途退出。研究表明,当纵向数据中存在不可忽略的缺失数据时,联合模型是量化中途退出与结局之间关系的有效手段。而且,对有不可忽略缺失的复杂纵向数据而言,这也是检测模型对基本假设敏感程度的好方法。
致谢
作者感谢他的同事Victor Santana博士和Clinton Stewart博士为本章的例题提供了宝贵的数据,Peter Song博士在连续型比例数据方面的合作以及Kevin Liu和Catherine Billups所做的数据分析。我们还要感谢美国国家癌症研究中心的资助(编号为P30 CA21765)以及美国,黎巴嫩,叙利亚相关慈善机构(ALSAC)的帮助。
附录
均值(位置参数)为 μ∈(0,1)和离散参数为 σ 2﹥0的单纯形分布的密度(Jorgensen [7])为:
p( y; μ, σ 2)=[2 πσ 2{ y(1- y)} 3] -1 / 2exp{- d( y; μ)/(2 σ 2)}, y∈(0,1)
其中
应用该分布的好处是:单纯形分布是一个分散模型(Jorgensen [7]),其中反应变量的密度函数式为:
a( y; σ 2)exp{- d( y; μ)/(2 σ 2)}, y∈(0,1)
该分散模型的密度看起来像正态分布(详见Jorgensen [7]),并且它还包括了一大类定义在(0,1)范围内的分布,从高度偏峰的分布到非常平坦的分布。这种分散模型要比基于指数族分布的广义线性模型更为一般。
参考文献
1.Moore DS,McCabe GP. Introduction to the Practice of Statistics. New York,1989.
2.Altman DG. Practical Statistics for Medical Research. Chapman and Hall,London,1991.
3.Gleser LJ. The importance of assessing measurement reliability in multivariate regression. Journal of the American Statistical Association,1994,87:696-707.
4.Carroll R,Ruppert D,Stefanski LA. Measurement Error in Nonlinear Models. Chapman and Hall,London,1995.
5.Chen Y,Dougherty ER,Bittner ML. Ratio-based decisions and the quantitative analysis of cDNA microarrays. Nature Genetics Supplement,1997,21:33-37.
6.Newton MA,Kendziorski CM,Richmond CS,et al. On differential variability of expression ratios:Improving statistical inference about gene expression changes from microarray data. Journal of Computational Biology,2001,8(1):37-52.
7.Jorgenson B. Dispersion Models Chapman and Hall/CRC,London,1997.
8.Song P,Tan M. Marginal model for continuous proportional data. Biometrics,2000,56:496-502.
9.Tan M. Using dispersion models in molecular pharmacology and genetics. Invited Presentation at Joint Statistical Meetings,Atlanta,GA,August 7,2001.
10.Diggle P,Heagerty P,Liang KY,et al. Analysis of Longitudinal Data. Oxford University Press,UK,2002.
11.Qu Y,Tan M. Analysis of clustered ordinal data with subclusters via a Bayesian hierarchical model. Communications in Statistics A: Theory & Method,1998,27:1461-1475.
12.Tan M,Qu Y,Mascha E,et al. A Bayesian hierarchical model for multi-level repeated ordinal Data:Analysis of oral practice examinations in a large anesthesiology training program. Statistics in Medicine,1999,18:1983-1992.
13.Myers BD,Nelson RG,Tan M,et al. Progression of overt nephropathy in non-insulin-dependent diabetes. Kidney International,1995,47:1781-1789.
14.Nelson RG,Bennett PH,Beck GJ,et al. Development and progression of renal disease in Pima Indians with non-insulin-dependent diabetes mellitus. New England Journal of Medicine,1996,335:1636-1642.
15.Song P,Qiu Z,Tan M. Modeling heterogeneous dispersion in marginal simplex models for continuous longitudinal proportional data. Biometrical Journal,2004,46:540-553.
16.Qiu Z,Song PXK,Tan M. Simplex Mixed-Effects Models for Longitudinal Proportional Data. Scandinavia Journal of Statistics,2008,35:577-596.
17.Concato J. Observational Versus Experimental Studies:What’s the Evidence for a Hierarchy? Neuro RX: The Journal of the American Society for Experimental Neuro Therapeutics,2004,1:341-347.
18.Shmueli G. To Explain or to Predict? Statistical Science,2010,25(3):289-310.
19.Helmer O,Rescher N. On the epistemology of the inexact sciences. Manag. Sci,1959,5:25-52.
20.Hastie T,Tibshirani R,Friedman JH. The Elements of Statistical Learning: Data Mining, Inference, and Prediction. 2nd ed. Springer,New York,2009.
21.Liu Z,Tan M. ROC based utility function maximization for feature selection and classification with application to high dimensional protease data. Biometrics,2008,64:1155-1161.
22.Liu Z,Chen X,Gartenhaus RB,et al. Survival Prediction and Gene Identification with Penalized Global AUC Maximization. Journal of Computational Biology,2009,16(12):1661-1670.
23.Box GEP,Cox DR. An analysis of transformations. Journal of the Royal Statistical Society Series B,1964,26:211-252.
24.Bacik J,Mazumdar M,Murphy BA,et al. The functional assessment of cancer therapy-BRM(FACT-BRM):a new tool for the assessment of quality of life in patients treated with biologic response modifiers. Qual Life Res,2004,13(1):137-154.
25.Du H,Hahn EA,Cella D. The impact of missing data on estimation of health-related quality of life outcomes:an analysis of a randomized longitudinal clinical trial. Health Services and Outcomes Research Methodology,2011,11:134-144.
26.Zhang D,Davidian M. Linear mixed models with flexible distributions of random effects for longitudinal data. Biometrics,2001,57:795-802.
27.McCulloch CE,Neuhaus JM. Misspecifying the Shape of a Random Effects Distribution:Why Getting It Wrong May Not Matter. Statistics Science,2011,26(3):358-402.
28.Hsieh F,Tseng YK,Wang JL. Joint modeling of survival and longitudinal data:likelihood approach revisited. Biometrics,2006,62:1037-1043.
29.Huang X,Stefanski LA,Davidian M. Latent-model robustness in joint models for a primary endpoint and a longitudinal process. Biometrics,2009,65:719-727.
30.Hosmer D,Lemeshow S. Applied Logistic Regression. 2 nd ed. John Wiley & Sons,INC. New York,2000.
31.Bland JM,Altman DG. Statistical methods for assessing agreement between two methods of clinical measurement. The Lancet,1986,1(8476):307-310.
32.Kosinski AS,Flanders WD. Evaluating the exposure and disease relationship with adjustment for different types of exposure misclassification:a regression approach. Stat Med,1999,18(20):2795-2808.
33.Nguyen A,Sarit A,Sand P,et al. Nongenetic factors associated with stress urinary incontinence. American Journal of Obstetrics & Gynecology,2011,117(2 Pt 1):251-255.
34.Neale MC,Boker SM,Xie G,et al. Mx:statistical modeling. 6th ed. Richmond(VA):Department of Psychiatry,Medical College of Virginia,2002.
35.Neale MC,Maes HM. Methodology for genetics studies of twins and families. Dordrecht(The Netherlands):Kluwer Academic,2004.
36.Hedeker D,Gibbons RD. Application of random-effects pattern-mixture models for missing data in longitudinal studies. Psychol Methods,1997,2(1):64-78.
主要作者简介
杜鸿雁,现为美国NorthShore医疗健康系统临床信息研究中心统计专家、美国统计协会认证专家(PStat®)。1993年毕业于华西医科大学,2003年获伊利诺伊大学芝加哥分校生物统计学硕士学位。在与健康结果相关的研究设计与数据分析方面具有丰富的合作研究经验,曾与癌症、糖尿病、慢性肾病、艾滋病以及其他复杂疾病的研究者进行过广泛的合作,在Cancer Research,Stroke,Clinical Trials,Statistics in Medicine and Health Services and Outcomes Research Methodology等发表论文70余篇。研究兴趣包括带不可忽略缺失数据的纵向模型、生物标记物评价(ROC等)以及预后建模。
谭铭博士,现为乔治城大学生物信息学、生物数学和生物统计系主任、教授。1990年获普度大学统计学博士学位。曾任Cleveland诊所生物统计和流行病学系助理教授、副教授,St. Jude儿童研究医院生物统计系准会员/教授及此医院的实体恶性肿瘤治疗计划的生物统计学主任,马里兰大学医学院生物统计学教授及生物信息学、生物统计学主任,马里兰大学Greenebaum癌症中心生物统计学主任。他是Biometrics和Statistics in Medicine副主编,美国FDA(食品和药物管理局)顾问委员会会员及多个美国国立卫生研究院专家组成员,美国统计协会资深会员。目前研究兴趣包括多种药物组合的设计与分析、生物标记物及诊断方法的评估、自适应临床试验设计、纵向数据的随机效应与贝叶斯模型。