580_3948826_混杂效应和随机效应模型.ppt.ppt

随机效应概念与混合效应模型(Concept of Random Effects and Mixed Effects Models),1,2,介绍内容一.统计模型的概念二.随机效应的概念与识别三.混合效应模型四.混合效应模型分析的例子,一.统计模型的概念,统计模型是对资料结构的一种数学表述.数量关系的概念化结构.包含两个元素1.函数表达式:描述结果变量与解释变量之间的关系(固定效应).2.误差表达式:描述结果变量观察值随机变异的概率分布(随机变异).例如:2种药物(A、B)治疗某种疾病的疗效分析。用均衡设计，每种药物治疗的病人数相等，都为n.反应变量：Yij表示生化测定值，i=1,为病例编号，J=1,2为药物编号 自变量：药物种类(A,B)，令Xj=第j种药物，传统的统计分析方法（固定效应模型，效应为常数）为：（1）用单向方差分析模型表示为:Yij=j+eij=+j+eij,eij N(0,e2),Yij N(i,e2)，j=j-,H0:j=0,限制条件：j=0（2）用线形回归模型表示为：Yij=0+iXij+eij,，Yij N(0+iXij,e2),H0:i=0,限制条件：B=0 含随机效应的混合效应模型为:Yij=(0+i)+iXij+eij,，i N(0,2),eij N(0,e2)这时 Yij N(0+iXij,i2+e2),Var(Yij)=Var(i)+Var(eij)=i+e2，,3,4,例1：A.B两种治疗药物在同一病人体内实验,采用区组随机化设计方案(即用药先后顺序是随机化的)，对每种药物处理后的反应变量进行测定.用6例病人.结果如下表.,构造三种模型:完全随机设计模型：不考虑区组（病人）效应:Yij=+j+eij,j 为药物效应随机化区组设计模型：考虑区组（病人）效应：Yij=+j+i+eij 随机效应模型：病人是从病人总体中随机的，也存在随机误差，统计学中用病人间的方差来衡量这种随机误差.。Yij=+j+(i)+eij=(+j+(j+eij),jN(0,2),eij N(0,e2)，Var(Yij)=(2+e2)在此简单情况下，(3)与(2)等价，但解释不同。在有缺失值情况下的结果不同。,5,模型一：完全随机设计模型:J：第J种药物效应,从上表估计模型参数:=20.75,A=22.83-20.75=2.08,B=18.67-20.75=-2.08 差值(difference)=22.83-18.67=4.17(或 A B=2.08-(-2.08),6,完全随机设计模型的PROC ANOVA 计算结果:,Sum ofSource DF Squares Mean Square F Value Pr FModel 1 52.0833333 52.0833333 2.68 0.1325Error 10 194.1666667 19.4166667Corrected Total 11 246.250000,Means with the same letter are not significantly different.SNK Grouping Mean N drug A 22.833 6 A A A 18.667 6 B,PROC ANOVA DATA=example_1;CLASS drug;MODEL y=drug;MEANS drug/SNK ALPHA=0.05;run;,7,7,完全随机设计模型的 PROC GLM 计算结果:,Sum ofSource DF Squares Mean Square F Value Pr FModel 1 52.0833333 52.0833333 2.68 0.1325Error 10 194.1666667 19.4166667 Corrected Total 11 246.2500000,StandardParameter Estimate Error t Value Pr|t|Intercept 18.66666667 B 1.79891943 10.38.0001()drug A 4.16666667 B 2.54405625 1.64 0.1325(A=4.17)drug B 0.00000000 B.(B=0.00),PROC GLM DATA=example _1;/*Model 1:completely randomized design model*/CLASS drug;MODEL y=drug/SOLUTION;RUN;,8,模型二:随机区组设计模型(考虑病人效应i):yij=+j+i+eij eijN(0,2),PROC GLM DATA=example_1;/*model 2:Randomized block design model*/CLASS drug patient;MODEL y=drug patient;RUN;,Sum ofSource DF Squares Mean Square F Value Pr FModel 6 206.8333333 34.4722222 4.37 0.0634Error 5 39.4166667 7.8833333 Corrected Total 11 246.2500000 Source DF Type III SS Mean Square F Value Pr Fdrug 1 52.0833333 52.0833333 6.61 0.0500patient 5 154.7500000 30.9500000 3.93 0.0798,(与模型一比较,残差均方(Mean Square(Error,2)由19.4166667降到7.88),9,模型三:病人为随机效应的模型:,因此,对同一病人的不同观察之间是相关的,具有协方差 2,包含在总方差Var(yij)=e2+2内,2和e2 都称为方差分量.但特别指2。,10,PROC MIXED DATA=example_1;/*Model 3:Random effects model by using PROC MIXED*/CLASS drug patient;MODEL y=drug;RANDOM patient/S;RUN;,用SAS中的 PROC MIXED 计算结果:,Cov Parm Estimatepatient 11.5333(用PROC GLM的RANDOM语句得不到此方差分量)Residual 7.8833(组内相关系数ICC=11.53/(11.53+7.88)=0.59)Type 3 Tests of Fixed Effects Num DenEffect DF DF F Value Pr Fdrug 1 5 6.61 0.0500,在本例中,对drug 的检验,用PROC MIXED的计算结果与用PROC GLM(2)的计算结果同（F=6.61）,即规定病人是固定效应,还是随机效应,对处理效应的检验结果没有影响（这是由于方差的性质决定的，即观察值的方差与中心化值的方差相等）.但如果有缺失值时,其结果不同.,在本例的模型三中,假定病人具有随机效应.病人来自一个具有均值为0,方差为2的正态分布总体.因此它们的期望值为0,但每个病人彼此不同。每个病人都具有相同期望值的假定与直观不符.须根据每例病人的观察值,确定其在正态分布中的一个位点.这一预报值的可信区间较固定效应的可信区间要窄,在统计学上称为收缩”shrunken”估计.,11,在本例的模型三中,假定病人具有随机效应.即规定病人来自一个具有均值为0,方差为2的正态分布总体.因此它们的期望值为0。但每个病人彼此不同。每个病人都具有同一期望值的假定与直观不符.须根据每例病人的观察值,确定其在正态分布中的一个位点.这一预报值的可信区间较固定效应的可信区间要窄,在统计学上称为收缩”shrunken”估计.这一收缩的幅度与病人方差分量和残差方差分量有关。当病人方差分量为0时，所有病人的预报值相等。对每个病人的观察值越少时，收缩的幅度相对越大。,随机效应模型的反应变量估计或预报,从上表可见，随机效应模型的预报值更接近观察值。,随机效应:在一项研究中，如果进入研究的因子的水平数只是其总体中的所有水平数的一个随机代表时,该因子的效应为随机效应。对应于该因子的总体中各水平的效应就构成了一个概率分布总体。样本中的各水平是来自总体中更多水平的一个随机样本.如:一个城市有很多学校,为了解学生体质,抽查了部分学校的学生体质，则所抽查的学校就是一个随机效应因子.一个城市有很多医院,为了解医疗质量,抽查了若干医院的出院病人记录进行医疗质量分析。则所抽查的医院就是一个随机效应因子.一条大河有许多支流的入口，为评价河水的亚硝酸盐浓度，只能抽查一小部分支流入口处的水样作检验，所抽查的支流入口处的水样就是一个随机效应因子.一个县市有很多村镇,为了解村民健康状况,随机抽查了若干村镇的村民，记录了他们的健康状况。则所抽查的这些村镇就是一个随机效应因子.研究随机效应的目的：1.估计随机效应的协方差参数。2.对总体参数作假设检验3.构造总体参数的可信区间。,12,二.随机效应的概念与识别,13,固定效应与随机效应的识别方法:当一个因子(预测变量)对反应变量的效应不易区别是固定效应还是随机效应时,可用可互换性或唯一性规则来作判断.可互换性(exchangeability)判别。一个随机效应因子的水平是随机地或非系统地选自具有更多水平的总体.观察样本中的水平只是总体中包含的更多水平中的一个随机样本.一个随机效应因子的水平是随机的,在不改变实验的基本性质情况下,当重 复实验时,其水平可能发生改变.固定效应因子:特意选择水平的因子.在不改 变实验基本性质的情况下,当重复实验时,其水平不发生改变.(2)从模型理论上区别:如果一个效应水平能够合理地假定为代表一种概率分布的话,则该效应就是 随机效应;如果不代表一种概率分布的话,则该效应就是固定效应.药物疗效试验中的药物品种是特选的,是不能互换的，故为固定效应。而在 药品价格调查中,每类药品选一种作为代表,这时调查的药品名称是可互换的，故为随机效应。世界上没有固定效应因子和随机效应因子之分,而是研究者在设计一个实验 时强加的一种不同结构的模型.以便更好地解释客观存在。,13,14,三、混合效应模型(Mixed effects Models),混合效应模型是一种线性模型，包含有固定效应和随机效应，用于处理非独立观察资料。又称：重复测量模型(Repeated measures models)，多水平模型(Multilevel models)，层次结构模型(Hierarchical models)从统计学归类，混合效应模型包含三种类型的模型:1.随机效应模型(Random effects models)。假定除测量误差导致的变异外，还来自具有某种概率分布的随机效应带来的变异，称随机变异。如分析临床多中心试验中的不同中心之间的变异。2.协方差类型模型(Covariance pattern models)，直接对重复测量之间的相关结构进行分析，分析效应随时间衰减的特点。3.随机系数模型(Random coefficient models)，直接分析反应变量在时间轴上的变化率，但容许协变量的效应具有随机性，协变量对反应变量效应的变化率随观察对象而不同，即具有随机变化的特点。以上三种类型的模型可以联合应用。,15,混合效应模型的方差协方差结构,其中:Y,为反应变量向量,X为固定效应因子的设计矩阵,为固定效应参数向量Z为随机效应因子的设计矩阵,为随机效应参数向量，i N(0,2),e为残差向量.eijN(0,e2),cov(i,e)=0,混合效应模型的参数估计:,16,16,例子：两种药物治疗效果的随机区组设计模型的矩阵表达,其中:Yij,为反应变量,j=1,2代表药物号；，2，。，6代表病例号X为（12行,3列）固定效应因子的设计矩阵,=(,1,2),为(3行1列)固定效应参数向量Z为(12x6)的随机效应因子的设计矩阵,i为(6x1)维随机效应参数向量e为(12x!)维残差向量.,Treatment：为固定效应，离散化为X矩阵Patient：为随机效应，离散化为Z矩阵,/*-Model 3_2:PROC MIXED-*/PROC MIXED DATA=intro;CLASS drug patient;MODEL y=drug;RANDOM patient;RUN;,17,=(，j,2),=(1,2,3,4,5,6),i(0,2),Var(yij)=(2+e2),18,例2:一项治疗高血压的多中心临床药物试验:三种药物(A,B,C),共有29所医疗中心参与,观察病人总数n=288人.(filename of dataset=hypertension)研究目的:在控制治疗前舒张压条件下,分析三种药物的降压效果.,四.混合效应模型分析的例子,19,变量名:记录号:patient,n=288,反应变量:dbp:治疗后舒张期血压,处理因素:treat:三种药物:A=Carvedilol,B=Nifedipine,C=Atenolol;控制因素:1.医疗中心:centre,29 所医院.2.治疗前舒张期血压:dbp1,连续变量.,Obs patient centre treat dbp dbp1 1 1 29 C 86 97 2 2 29 C 72 109 3 3 5 B 109 117 4 4 5 A 87 100 5 5 29 A 85 105 6 7 3 A 100 114 7 8 3 B 80 105 8 9 3 B 90 100 9 10 3 A 100 102 10 11 3 C 94 105,Variable N Mean Minimum Maximum-dbp 288 90.2465278 70.0000000 140.0000000dbp1 288 102.8541667 92.0000000 120.0000000,数据集中的前10例病人的记录,20,模型A:简单药物效应treatk:k=A,B,C.（固定效应）,PROC MIXED DATA=hypertension;/*MODEL A*/CLASS centre treat;MODEL dbp=treat/SOLUTION;RUN;,Covariance Parameter EstimatesCov Parm Estimate Residual 81.5660 Fit Statistics-2 Res Log Likelihood 2076.9AIC(smaller is better)2078.9AICC(smaller is better)2078.9BIC(smaller is better)2082.5,Solution for Fixed Effects StandardEffect treat Estimate Error DF t Value Pr|t|Intercept 88.6211 0.9266 285 95.64 Ftreat 2 285 2.73 0.0670,21,21,模型B:在模型A的基础上加入基础血压dbp1.（固定效应）,PROC MIXED data=hypertension;/*model B/CLASS centre treat;MODEL dbp=treat dbp1/SOLUTION;RUN;,Solution for Fixed Effects StandardEffect treat Estimate Error DF t Value Pr|t|Intercept 58.1490 11.3666 284 5.12 Ftreat 2 284 2.92 0.0558dbp1 1 284 7.23 0.0076,Covariance Parameter EstimatesCov Parm EstimateResidual 79.8201 Fit Statistics-2 Res Log Likelihood 2072.3AIC(smaller is better)2074.3AICC(smaller is better)2074.3BIC(smaller is better)2078.0,Covariance Parameter EstimatesCov Parm Estimate Residual 81.5660(model A)Fit Statistics-2 Res Log Likelihood 2076.9AIC(smaller is better)2078.9AICC(smaller is better)2078.9BIC(smaller is better)2082.5,22,Covariance Parameter EstimatesCov Parm EstimateResidual 71.9213 Fit Statistics-2 Res Log Likelihood 1892.6AIC(smaller is better)1894.6AICC(smaller is better)1894.6BIC(smaller is better)1898.2,Type 3 Tests of Fixed Effects Num DenEffect DF DF F Value Pr Fdbp1 1 256 3.87 0.0501treat 2 256 2.96 0.0535centre 28 256 2.11 0.0013,PROC MIXED DATA=hypertension;/*model C*/CLASS centre treat;MODEL dbp=dbp1 treat centre/SOLUTION;RUN;,模型C:在模型B的基础上加入医疗中心centre.（固定效应）,Covariance Parameter EstimatesCov Parm EstimateResidual 79.8201(model B)Fit Statistics-2 Res Log Likelihood 2072.3AIC(smaller is better)2074.3AICC(smaller is better)2074.3BIC(smaller is better)2078.0,23,Solution for Fixed Effects StandardEffect treat centre Estimate Error DF t Value Pr|t|Intercept 65.5796 12.9401 256 5.07.0001dbp1 0.2230 0.1133 256 1.97 0.0501treat A 2.9907 1.2336 256 2.42 0.0160treat B 1.7937 1.2651 256 1.42 0.1574treat C 0.centre 1 2.0297 4.4980 256 0.45 0.6522centre 2-4.5780 5.0511 256-0.91 0.3656Centre 3 6.1361 5.1958 256 1.18 0.2387Centre 4 0.3923 4.9744 256 0.08 0.9372centre 5 3.5446 4.9810 256 0.71 0.4773 Centre 35-1.6042 6.0736 256-0.26 0.7919Centre 36-1.1272 4.6359 256-0.24 0.8081centre 37-2.5290 5.5541 256-0.46 0.6492centre 40 6.2783 7.3758 256 0.85 0.3955centre 41 0.,24,模型D:在模型C基础上加入医疗中心与治疗的交互作用:centre*treat.（固定效应）,PROC MIXED DATA=hypertension;/*model D*/CLASS centre treat;MODEL dbp=dbp1 treat centre centre*treat/SOLUTION;TITLE MODEL D;RUN;,Covariance Parameter EstimatesCov Parm EstimateResidual 69.2614 Fit Statistics-2 Res Log Likelihood 1558.1AIC(smaller is better)1560.1AICC(smaller is better)1560.1BIC(smaller is better)1563.5,Type 3 Tests of Fixed Effects Num DenEffect DF DF F Value Pr Fdbp1 1 208 0.99 0.3198treat 2 208 1.24 0.2905centre 28 208 1.98 0.0038centre*treat 48 208 1.20 0.1884,交互效应centre*treat作用项不显著,Covariance Parameter EstimatesCov Parm EstimateResidual 71.9213(model C)Fit Statistics-2 Res Log Likelihood 1892.6AIC(smaller is better)1894.6AICC(smaller is better)1894.6BIC(smaller is better)1898.2,25,25,模型E:在模型B的基础上,将医疗中心centre作为随机效应.（与固定效应比较）,Type 3 Tests of Fixed Effects Num DenEffect DF DF F Value Pr Fdbp1 1 256 6.84 0.0095treat 2 256 3.10 0.0466,PROC MIXED DATA=hypertension;/*model E*/CLASS centre treat;MODEL dbp=dbp1 treat/SOLUTION;RANDOM centre/SOLUTION;RUN;,Covariance Parameter EstimatesCov Parm Estimatecentre 7.8248Residual 70.9263Fit Statistics-2 Res Log Likelihood 2056.5,Covariance Parameter EstimatesCov Parm EstimateResidual 79.8201(model B)Fit Statistics-2 Res Log Likelihood 2072.3,Covariance Parameter EstimatesCov Parm Estimate(model D)Residual 69.2614 Fit Statistics-2 Res Log Likelihood 1558.1,五种模型的配合结果比较：,注：X-无统计学意义，?-接近=0.05水平，*-在=0.05水平上有统计学意义,27,ODS GRAPHICS ON;PROC MIXED DATA=hypertension;CLASS centre treat;TITLE Residual Plot of Model E;MODEL dbp=treat dbp1/DDFM=KR RESIDUAL;RANDOM centre/SOLUTION;LSMEANS treat/DIFF PDIFF CL;ID patient centre treat;RUN;ODS GRAPHICS OFF;,模型E的残差分析:,28,例3:三水平混合效应模型的例子。一个制药厂为了解生产的稳定性，进行了抽样研究。两种原料，从原料(1)中抽取4批产品，从原料(2)中抽取4批产品，每批产品中随机抽取3分样品，再从每份样品中随机取出3分检验样品。总检验数为：2x4x3x3=72个检验品。抽样过程为:,原料(source,i=1,2)每种原料抽取4批产品(lot,j=1,2,3,4)每批产品中抽取3分样品(samp,k=1,2,3)再从每份样品中取出3分检验品(lab,m=1,2,3),检验项目：有效物质含量(y)。检验目的：分析产品不稳定当 来源。,30,资料表：,31,proc mixed data=semiconductor method=reml;class source lot sample;model content_y=source/ddfm=kr;random lot(source)sample(source lot);lsmeans source/diff;run;,原料(source,i=1,2)为固定效应每种原料抽取4批次的产品(lot,j=1,2,3,4)为随机效应(_嵌套于原料)每批产品中抽取3分样品(samp,k=1,2,3)为随机效应(_嵌套于原料和批次)再从每份样品中取出3分检验品(lab,m=1,2,3)的化验结果为content_y(反应变量),Covariance Parameter EstimatesCov Parm Estimatelot(source)119.89sample(source*lot)35.8657Residual 12.5694,结论1：批次之间的变异性最大，为样品变异的20倍(=119.89/35.8657),为样品内变异的46倍(=119.86/12.5694)品(lab,m=1,2,3)的化验结果为content_y(反应变量),32,Least Squares Means StandardEffect source Estimate Error DF t Value Pr|t|source 1 1995.11 5.7716 6 345.68|t|source 1 2-10.0833 8.1622 6-1.24 0.2629,结论2：原料之间对产品的质量变异的影响不大Pr|t|=0.2629),33,4:配对病例对照的例子。婴儿卒死综合症（sudden infant death syndrome，SID）的母亲与同期同一所医院内相同年龄组事件前分娩相同年龄组母亲和事件后分娩相同年龄组母亲各一人作1:2配对调查。研究指标为deprivation score（depcat,1-7分，分值越高者越严重）。由于各种原因，配对不是完整的。资料的配对情况为:,The FREQ Procedure Cumulative Cumulativetot Frequency Percent Frequency Percent-A 45 22.39 45 22.39AB 36 17.91 81 40.30ABB 120 59.70 201 100.00,注：A：病例，B:对照,34,Obs id group grp depcat 1 1 A 1 7 2 2 A 1 7 3 3 A 1 5 4 4 A 1 7 5 5 A 1 5 6 6 A 1 2 7 6 B 0 2 8 6 B 0 2 9 7 A 1 4 10 7 B 0 4 11 7 B 0.12 8 A 1 4 13 8 B 0 5 14 8 B 0 6 15 9 A 1 1 16 10 A 1 7 17 11 A 1 6 18 11 B 0 4 19 11 B 0 4 20 12 A 1 7,前20组匹配资料：1-5组只有病例，没有对照。第7匹配组中有一例对照的资料缺失。,35,PROC MIXED DATA=exsample_3 NOCLPRINT;CLASS group id;MODEL depcat=group/DDFM=SATTERTH outp=op outpm=opm;RANDOM ID;*/SOLUTION;ESTIMATE A-B group 1-1;ID id group;title Model 1;run;,DimensionsCovariance Parameters 2Columns in X 3Columns in Z Per Subject 1Subjects 201Max Obs Per Subject 3Observations Used 461Observations Not Used 16Total Observations 477,Covariance Parameter EstimatesCov Parm Subject EstimateIntercept id 0.9576Residual 1.8602,组内相关系数：0.9576/(1.8602+0.9576)=0.34,36,Type 3 Tests of Fixed Effects Num DenEffect DF DF F Value Pr Fgroup 1 315 40.19|t|A-B 0.8446 0.1332 315 6.34.0001,统计学结论：婴儿卒死综合症（sudden infant death syndrome，SID）与母亲deprivation score 有关。病例的母亲平均记分高于对照母亲。,37,谢谢大家,