二分类资料的多水平模型.ppt

二分类资料的多水平分析,张菊英 教授,问 题 背 景,数据的层次结构(hierarchical structure)现象是广泛存在的，这种结构可以是自然的，亦可以是人为形成的。,试验研究,致畸试验 常用孕鼠作试验 ，观察每个孕鼠所产子代中发生畸形的情况 。,层次结构：,孕鼠1,孕鼠2,孕鼠p,窝别效应,某省调查其农村居民的卫生服务 随机抽取30个乡镇，每个乡镇分别抽取2个行政村，每个村再随机抽取33户（家庭），对每个家庭前半年内的常住人口进行问卷调查。,调查研究,资料的特点,具有明显的层次结构(乡镇 行政村 户 个体) ； 在经济水平、生活方式、生活习惯上都具有某种程度上的相似性或聚集性 ；个体的数据是非独立的 。,通常处理方式,结果变量为连续性或定量 如研究某药物对仔鼠体重的影响，可用方差分析分离出区组间效应（如窝别效应），然后比较试验组与对照组仔鼠体重的差异,通常处理方式,结果变量为二项分布 两组的畸形发生率，用两个率的差别的Z检验或 检验进行比较 如需要控制混杂因素，采用多因素的logistic回归等方法 。,通常处理方式,不足之处 忽略了数据层次结构的特征,二分类反应变量单水平模型,二分类反应变量的单水平模型,例5.1 某公共卫生学院营养与食品卫生学教研室用孕鼠进行口服花粉的致畸实验，将26 只孕鼠随机分为甲、乙两组，甲组在孕早期每天给予辐射花粉2000mg/kg , 乙组每天给予500mg/kg 。待孕鼠分娩后，观察其子鼠骨骼畸形的发生情况。,不同处理的孕鼠所产仔鼠骨骼畸形发生结果,注：分母表示该窝仔鼠总数，分子表示仔鼠中骨骼畸形发生数,忽略窝别效应，直接进行两组仔鼠畸形发生率的比较，可采用单水平的logistic回归,模型,为应变量：1 表示发生畸形 0 表示未发生畸形 为处理因素 ：1 表示高剂量花粉(甲组) 0 表示低剂量花粉(乙组),各软件对例5.1拟合单水平logistic回归模型结果,结果解释,甲组相对于乙组，仔鼠发生畸形的比数比的对数值(ln OR)为0.773，P0.043。 由结果可见，忽略窝别效应，处理因素(甲、乙两种喂养方式)对结果的影响有统计学意义，高剂量花粉比低剂量花粉更容易发生仔鼠畸形。,讨论,对于例5.1，单水平logistic回归分析处理效应的检验结果P值接近检验水准0.05，这一结果可靠程度值得怀疑,讨论,试验中仔鼠畸形在窝别间存在聚集性，仔鼠发生畸形的概率不是完全独立的，数据不再服从二项分布 若仍用单水平logistic回归处理这类资料，则会低估两个率差别的标准误，从而增大犯型错误的概率，即将本来无差别的两个率判为有差别,广义线性模型,实际上logistic回归模型只是广义线性模型(generalized linear model，GLM)的一种 广义线性模型是一般线性模型的直接推广 很多模型属于非线性模型，如指数模型、logistic回归模型等，它们通过一定的变量变换，可以转化成线性模型，并满足或近似满足线性模型分析的要求,广义线性模型的一般形式,是已知的设计矩阵 是待估的未知参数 是连接函数,广义线性模型表示为以下三部分,随机部分： 结果变量Y，服从指数族概率分布，,系统部分,通过协变量 产生线性预测值 ，,连接函数,随机部分和系统部分由下式联系在一起 其中 被称为连接函数。,连接函数有多种，应根据应变量的特点加以选择。 应变量与不同的连接函数就构成了不同的回归模型,应变量为二项分布时常用的连接函数,probit模型,probit模型在教育领域应用较多 probit模型中，与 相应的是正态分布的下侧累计概率函数。因此系数 的解释是当其它自变量不变时，自变量改变一个单位时，所导致的“概率单位”的改变。,probit模型,概率单位在数学上的意义是很明确的，但在生物学上没有合适的解释，因此probit回归在流行病学应用上尚未普及的原因 对多数资料来说，logistic回归模型与probit模型的估计概率相近,各软件对例5.1拟合单水平probit回归结果,二分类反应变量两水平模型,优势,处理具有层次结构特征的数据资料，可将传统模型中的随机误差项分解到与数据层次结构相应的水平上，使得个体的随机误差更纯 。,二项分布,当反应变量为比数(率)，一般采用二项分布。 标准二项分布假定的比数(率)：,i指示水平1单位，j指示水平2单位， 为比数(率)的分母，在未分组数据的情况下， 。,两水平logit模型,两水平logit模型,为处理因素的效应参数，又称固定效应(fixed effect)参数 为水平2单位的logit均值 与总均值 之差，又称为随机效应(random effect)或高水平的残差。,两水平logit模型,的方差 又称为随机参数(random coefficient)，反映了高水平单位间的比数(率)的差别。 越大说明数据在高水平单位内的聚集性越强。 为0时，该模型演变为一般的logistic回归模型。,两水平logit模型,为尺度参数。当反应变量确定服从二项分布，则尺度参数应该为1或接近1。即当模型的固定效应参数和随机效应参数的估计值确定后，反应变量的方差估计值为,两水平logit模型,拟合模型时，若假设二项方差成立，则设置尺度参数 为1 可允许 为待估参数，对水平1方差是否“超二项变异” 进行检验，即考察水平1方差是否满足二项分布的假定，可根据估计的尺度参数值和1.0的差值与的估计标准误之比作正态性Z检验而得。,判断是否存在高水平效应,密切结合专业知识和具体情况进行判断 对随机参数 的估计值做检验 用VPC(variance partition coefficient)来进行度量,VPC,当应变量为连续型变量时，VPC等价于组内相关系数(intra-class correlation) 在两水平的方差成分模型中， VPC表示了水平2的方差占总方差的比例,VPC,当应变量为离散型变量时，VPC与ICC不等价 以二项分布的资料为例，水平1的方差依赖于模型中解释变量的值，因此没有一个简单的方法来计算VPC,假设 m个水平2方差 的样本 在已知的固定效应 估计下，m个水平2单位的方差均值的经验估计为：,水平1误差的经验值,操作过程,1. 产生m个正态分布的随机数，记为r，例如m5 000，通常大样本为佳,2. 确定变量的值，根据已拟合的模型计算模型中的线性部分(L指示的部分)，由此获得m个所估计模型的样本 方差成分模型，L为 随机斜率模型，L为,3. 计算m个模型样本的估计反应率,4. 对m个计算经验方差获得,5.计算m个的水平1的二项分布误差，然后取其平均值得到水平1误差的经验值，,6.计算VPC,两水平模型的应用,对例5.1，在进一步考虑了“窝别效应”后，该资料可看作具有两个水平的层次结构，即：孕鼠水平(水平2 单位)与仔鼠水平(水平1 单位) ；甲、乙两种喂养方式是解释变量；仔鼠是否发生畸形是反应变量。,各软件对例5.1拟合方差成分模型结果,以MLwiN的拟合结果为例，水平2 单位(不同孕鼠之间)的随机效应残差(即窝别效应)方差为3.423，标准误为1.554， Wald检验结果为4.852，P0.028。可以认为不同窝别的仔鼠骨骼畸形的发生概率是不同的，“窝别效应”存在,在考虑了“窝别效应”后，处理因素的系数为1.078，标准误为0.993, OR 2.939，由于P0.05，因此，尚不能认为给孕鼠喂养大剂量的辐射花粉会引起仔鼠骨骼畸形发生增加,进一步分析窝别效应的强度,当处理因素取值为0时，VPC0.293，因此在乙组中，有29.3的残差方差归因于窝别效应 当处理因素取值为1时，VPC0.344，说明在甲组中，有34.4的残差方差归因于窝别效应,由此例分析可见，当数据存在随机效应时，忽略此效应拟合的单水平模型增大了处理效应的I型误差，使得本来无统计学意义的因素误认为有统计学意义。,通过解释变量(本例为甲、乙两个处理组)的各种组合，可以得到VPC的取值范围。本例VPC为0.2930.344之间。因此大概有30的残差变异归因于窝别效应。,三水平数据结构应用,例5.2,某省进行了农村贫困居民的家庭卫生服务调查。先随机抽取乡镇，每个乡镇分别抽取2个行政村，每个村再随机抽取一定数量的家庭进行调查。共调查了30个乡镇，832户贫困家庭户，共计2369名15岁及以上的居民。 现拟探讨该省农村贫困居民卫生服务需要的影响因素。,以两周是否患病作为应变量。结合资料的层次结构特点，采用二分类多水平logit回归模型探讨农村贫困居民两周是否患病的影响因素。,变量赋值表（1）,变量赋值表（2）,变量赋值表（3）,家庭人均居住面积和乡镇人均可耕地面积以连续性变量形式中心化之后纳入模型； 无序多分类变量婚姻状况和职业以哑元形式纳入； 有序多分类变量年龄、文化程度和自身健康状况评价，是以分组线性变量或哑元的形式纳入，依据似然比检验结果加以判断。,有序多分类变量纳入模型方式检验结果,年龄和文化程度以哑元形式纳入；自我健康状况评价以分组线性变量形式纳入模型。,哑变量名定义,age1：45，age2：65 ； marriage1：已婚，marriage2：离婚， marriage3：丧偶； education1：小学，education2：初中及以上； occupation1：农业劳动者，occupation2：学生， occupation3：离退休 occupation4：无业、失业、半失业,特别注意,MLwiN拟合模型之前，应对数据从高水平到低水平逐步排序，即按照乡镇、户、个人的顺序对全部数据进行排序。,拟合二水平不含任何解释变量的零模型,以户作为高水平即水平2，个体作为低水平即水平1，采用2阶 PQL算法（后面具体讲解）拟合二水平不含任何解释变量的零模型，固定尺度参数为1。,两周是否患病两水平零模型,零模型的水平2方差具有统计学意义，结合专业知识，数据在高水平上具有聚集性，其层次结构不能忽略。故进一步引入解释变量拟合二水平方差成分模型。,两周患病两水平方差成分模型（1）,结果解释,在控制其他因素不变的情况下，年龄在4564岁组和65岁及以上组，两周患病率均高于1544岁人群。 女性两周患病率高于男性。 半年内患有慢性病的居民两周患病率高于未患慢性病者。,结果解释,经常饮酒的居民两周患病率高于不饮酒或偶尔饮酒者 已婚、离异的居民两周患病率高于未婚者，丧偶者的患病率与未婚者相比未见统计学差异 随着家庭人均住房面积增大，两周患病率降低,结果解释,小学文化的居民两周患病率高于文盲半文盲者，初中及以上文化程度居民的患病率与文盲比较差异无统计学意义,结合例5.2资料尚具有更高层次结构(即乡镇水平)的特点，可构造三水平零模型(3-level model),三水平零模型,三水平零模型,两周患病3水平方差成分模型,虽然高水平方差具有统计学意义，但是否在模型中保留乡镇水平应根据研究者实际关注的重点来加以取舍。,以上模型估计的协变量的系数均为固定的，若在考虑这些协变量系数的固定效应的同时考虑其随机效应，即协变量对反应变量的效应在不同的高水平单位间是不同的，此时的模型即为随机系数模型(random coefficient model),两水平随机系数模型,随机系数模型与方差成分模型的区别在于，方差成分模型中 的系数固定为 。随机成分模型中假定 的效应在水平2单位间随机变化，且服从均数为 ，方差为 的边际正态分布。 的随机效应与截距的随机效应存在协变异时，用 来衡量。,注意,模型随机部分的解释变量可以是固定部分的一个子集，也可以不是，即可以在模型的固定部分或随机部分纳入任何水平上测量的解释变量。,参数估计方法,参数估计方法,多水平模型中的参数估计采用“迭代广义最小二乘法”(Iterative Generalized Least Squares，IGLS)或者“限制性迭代广义最小二乘法”(Restricted Iterative Generalized Least Squares，RIGLS，Goldstein)。,参数估计方法,当模型的随机变量在每个水平上均服从多变量正态分布，则IGLS等价于最大似然估计，RIGLS等价于限制性极大似然估计(Restricted Maximum Likelihood，REML)。,其它的参数估计方法,Longford 1987年提出的基于“费歇尔得分”的算法 ； Raudenbush 1994年证明它等价于IGLS 。 Liang和Zeger 1984年提出的广义估计方程(Generalized Estimating Equations，GEE) ； “马尔科夫链蒙特卡罗”(Markov Chain Monte Carlo，MCMC)方法，尤其是吉布斯抽样(Gibbs Sampling)的发展，使完全贝叶斯技术在计算上变得更为可行，进而完善了小样本的参数估计方法。,参数估计方法,由于不知道 的真实值，于是利用每次迭代中参数的当前估计值 来预测 ，再计算的 值 ，由于迭代中仅用到二项分布的均值和方差进行估计，所以这种估计方法称为“拟似然法”(quasi-likelihood)。,参数估计方法,建立线性化模型，形成泰勒级数展开时 一种算法是将当时的水平2残差估计值加到非线性函数的新型成份中，即“预测性拟似然法”(predictive quasi-likelihood，PQL)； 一种算法则不纳入这些估计值，称为“边际拟似然法”(marginal quasi-likelihood，MQL)。,参数估计方法,MQL方法计算速度较快，但只考虑了固定效应，当高水平单位具有较大方差且其低水平单位数较小的时候，趋向于低估固定和随机参数的值； PQL方法同时利用了水平二的残差，计算具有较小的偏性，但PQL算法不稳定，在某些情况下不容易收敛，且如果水平2内包含的水平1单位数偏少时，PQL估计的固定效应虽然无误，但水平2变异的估计值却可能是有偏的。,三种不同的随机效应回归模型,它们都用来描述随机效应的大小 Beta二项回归模型。Crowder 1978提出此线性模型。该模型假定同一处理组不同区组发生某事件的概率服从beta分布，则第j组中第i个区组事件发生数的边缘分别服从beta二项分布。 logistic正态回归模型。Pierce和Sands 1975年提出此线性模型。该模型假定随机效应与固定效应一样，与logit p 呈线性关系，且服从与标准正态分布。,三种不同的随机效应回归模型,logistic二项回归模型。Mauritsen1984年提出此线性模型。该模型假定随机效应与logit p呈线性关系，且服从二项分布。,与多水平模型的区别,这些模型只关注于用一个参数反映随机效应的综合效应，以便得到较好的固定效应估计值及其方差。它们不能有效地分析方差成分，也不能对各水平上的随机效应或参数估计做出有意义的解释。,双反应变量的多水平模型,对例5.2的资料拟合了多水平的logistic回归模型，了解农村贫困居民两周患病与否的影响因素，发现女性、年龄在45岁以上、经常饮酒者的患病率高于其对照组居民，而且半年内患有慢性病的居民其两周患病率显著高于未患慢性病者。,半年内患慢性病为医生诊断的患慢性病现状，两周患病为自己报告的近期急性病或身体不适的反映，两者均为反应变量。我们可进一步探讨两者的关系及其影响因素。例如，与两周患病率有关的因素，即性别、年龄和饮酒等是否也与半年内患有慢性病有关？其相关的程度是否相似？,自报告两周患病与半年内患慢性病诊断的相关性,两周患病率为24.7(590/2 388)，半年内慢性病患病率 为18.8(449/2 388)。 用Kappa系数和Pearsons列联系数来测量两反应变量 的相关性分别为0.468和0.475，表明两者间有相关性。,两个二分类变量的比例关系,则可知两变量均为阳性时的概率 两周患病为阳性的概率 半年内患慢性病为阳性的概率,双变量的二项分布可总结为有两个参数 和 及一个协变异参数的联合分布，其协变异参数表达为：,其值的大小反映两反应变量间相关性的强弱，理论上与Pearsons相关系数等同,对这两个有相关性的反应变量拟合多水平模型的基本思想是它们的数据具有基本的两水平结构：重复测量值(水平1)嵌套于个体内(水平2),最简单的两个二分类反应变量的模型可由以下两水平模型表达：,（式5.7）,上述模型可以直接推广到以下几种数据的模型分析,（1）研究同一协变量对不同反应变量相应的影响效应及其效应间的差别,对于例5.2的资料，前面的分析发现女性报告的两周患病率高于男性。 问题： 1.是否女性的半年内慢性病的患病率也高于男性呢？ 2.由于女性可能普遍较男性对身体发生不适 更敏感，这是否会导致自报的两周患病率对女性的影响更大于对半年内慢性病的影响？即女性两周患病的OR值大于女性半年内患慢性病的OR值呢？,构造扩展模型,性别两周患病的OR由 估计 半年内患慢性病的OR由 估计 性别效应在两个反应变量间的差别则可用广义的Wald检验：,如果已知性别效应在两反应变量间无差别，尚可在上述扩展的模型中限定(constraint)令 。 多个协变量可直接加到模型中,用MLwiN拟合的双反应变量logistic回归模型估计值,由模型(5.7)的估计值可见，个体内两反应变量的协方差估计值即为Pearsons相关系数值，对截距估计值分别进行反logit变换则获得原始的人群患病率24.7(两周患病)和18.8(半年内患慢性病),模型(5.8)的估计值提示女性和经常饮酒居民的两周患病的危险性大于对照组(男性，很少饮酒，年龄小于45岁)，而这种危险性并不存在于半年内患慢性病的情况。年龄效应对两种患病测量则相似。,对性别与饮酒效应在两种患病指标间的差别做假设检验获得Wald值分别为7.42(P0.006 5)和17.96(P0.000 1)。可以认为女性和经常饮酒者对近期身体不适较之于慢性病对身体的影响更为敏感,(2) 数据具有更高水平的层次结构,本章前面的单反应变量模型已显示个体两周患病在家庭内有聚集性，进而在乡镇内有聚集性。这里的家庭和乡镇是在个体水平上的高层水平单位，使得两反应变量在本例的数据有四个层次的结构：个体内重复测量单位 个体