混杂效应和随机效应模型

随机效应概念与混合效应模型

(ConceptofRandomEffectsandMixedEffectsModels)112

介绍内容一.统计模型的概念二.随机效应的概念与识别三.混合效应模型四.混合效应模型分析的例子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

，334Patient(i)yijDifference(yi1–yi2))

ӯi(PatientMean)A(j=1)B(j=2)12012816.022624225.031617-116.542921825.052221121.562417720.5Mean22.8318.674.1720.75例1：A.B两种治疗药物在同一病人体内实验,采用区组随机化设计方案(即用药先后顺序是随机化的)，对每种药物处理后的反应变量进行测定.用6例病人.结果如下表.构造三种模型:完全随机设计模型：不考虑区组（病人）效应:Yij=μ+βj

+eij

,βj为药物效应随机化区组设计模型：考虑区组（病人）效应：Yij=μ+βj

+αi+eij

随机效应模型：病人是从病人总体中随机的，也存在随机误差，统计学中用病人间的方差来衡量这种随机误差.。

Yij=μ+βj

+(γi)+eij

==(μ

+βj+

(γj+eij

),γj~N(0,τγ2),eij

~N(0,σe2)，Var(Yij)=

(τγ2+σe2)

在此简单情况下，(3)与(2)等价，但解释不同。在有缺失值情况下的结果不同。45模型一：完全随机设计模型:

β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))PatientTreatmentDifference(A–B)PatientMeanAB12012816.022624225.031617-116.542921825.052221121.562417720.5Mean22.8318.674.1720.7556完全随机设计模型的PROCANOVA计算结果:

SumofSourceDFSquaresMeanSquareFValuePr>FModel152.083333352.08333332.680.1325Error10194.166666719.4166667CorrectedTotal11246.250000Meanswiththesameletterarenotsignificantlydifferent.SNKGroupingMeanNdrugA22.8336AAA18.6676BPROCANOVADATA=example_1;CLASSdrug;MODELy=drug;MEANSdrug/SNKALPHA=0.05;run;677完全随机设计模型的PROCGLM计算结果:

SumofSourceDFSquaresMeanSquareFValuePr>FModel152.083333352.08333332.680.1325Error10194.166666719.4166667CorrectedTotal11246.2500000StandardParameterEstimateErrortValuePr>|t|Intercept18.66666667B1.7989194310.38<.0001(μ)drug

A4.16666667

B2.544056251.640.1325

(βA=4.17)drugB0.00000000B..(βB=0.00)PROCGLMDATA=example_1;/*Model1:completelyrandomizeddesignmodel*/

CLASSdrug;

MODELy=drug/SOLUTION;RUN;78模型二:随机区组设计模型(考虑病人效应αi):yij=μ+βj+αi+eijeij~N(0,σ2),PROCGLMDATA=example_1;/*model2:Randomizedblockdesignmodel*/

CLASSdrugpatient;

MODELy=drugpatient;RUN;SumofSourceDFSquaresMeanSquareFValuePr>FModel6206.833333334.47222224.370.0634Error539.41666677.8833333CorrectedTotal11246.2500000SourceDFTypeIIISSMeanSquareFValuePr>Fdrug

1

52.083333352.0833333

6.610.0500patient5154.750000030.95000003.930.0798(与模型一比较,残差均方(MeanSquare(Error,σ2)由19.4166667降到7.88)89模型三:病人为随机效应的模型:因此,对同一病人的不同观察之间是相关的,具有协方差σγ2,包含在总方差Var(yij)=σe2+σγ2内,σγ2和σe2

都称为方差分量.但特别指σγ2。910PROCMIXEDDATA=example_1;/*Model3:RandomeffectsmodelbyusingPROCMIXED*/

CLASSdrugpatient;

MODELy=drug;

RANDOMpatient/S;RUN;用SAS中的PROCMIXED计算结果:CovParmEstimatepatient11.5333(用PROCGLM的RANDOM语句得不到此方差分量)Residual7.8833(组内相关系数ICC=11.53/(11.53+7.88)=0.59)

Type3TestsofFixedEffectsNumDenEffectDFDFFValuePr>Fdrug1

5

6.61

0.0500在本例中,对drug的检验,用PROCMIXED的计算结果与用PROCGLM(2)的计算结果同（F=6.61）,即规定病人是固定效应,还是随机效应,对处理效应的检验结果没有影响（这是由于方差的性质决定的，即观察值的方差与中心化值的方差相等）.但如果有缺失值时,其结果不同.在本例的模型三中,假定病人具有随机效应.病人来自一个具有均值为0,方差为σα2的正态分布总体.因此它们的期望值为0,但每个病人彼此不同。每个病人都具有相同期望值的假定与直观不符.须根据每例病人的观察值,确定其在正态分布中的一个位点.这一预报值的可信区间较固定效应的可信区间要窄,在统计学上称为收缩”shrunken”估计.1011在本例的模型三中,假定病人具有随机效应.即规定病人来自一个具有均值为0,方差为σα2的正态分布总体.因此它们的期望值为0。但每个病人彼此不同。每个病人都具有同一期望值的假定与直观不符.须根据每例病人的观察值,确定其在正态分布中的一个位点.这一预报值的可信区间较固定效应的可信区间要窄,在统计学上称为收缩”shrunken”估计.这一收缩的幅度与病人方差分量和残差方差分量有关。当病人方差分量为0时，所有病人的预报值相等。对每个病人的观察值越少时，收缩的幅度相对越大。随机效应模型的反应变量估计或预报观察值与完全随机设计固定效应模型预报值及随机效应模型预报值的比较病人号123456drugAB均值AB均值AB均值AB均值AB均值AB均值观察值201216.0262425.0161716.5292125.0222121.5241720.5固定效应预报22.818.720.822.818.720.822.818.720.822.818.720.822.818.720.822.818.620.8随机效应预报19.315.117.226.021.823.919.715.517.626.021.823.923.419.221.322.718.520.6从上表可见，随机效应模型的预报值更接近观察值。11随机效应:在一项研究中，如果进入研究的因子的水平数只是其总体中的所有水平数的一个随机代表时,该因子的效应为随机效应。对应于该因子的总体中各水平的效应就构成了一个概率分布总体。样本中的各水平是来自总体中更多水平的一个随机样本.如:一个城市有很多学校,为了解学生体质,抽查了部分学校的学生体质，则所抽查的学校就是一个随机效应因子.一个城市有很多医院,为了解医疗质量,抽查了若干医院的出院病人记录进行医疗质量分析。则所抽查的医院就是一个随机效应因子.一条大河有许多支流的入口，为评价河水的亚硝酸盐浓度，只能抽查一小部分支流入口处的水样作检验，所抽查的支流入口处的水样就是一个随机效应因子.一个县市有很多村镇,为了解村民健康状况,随机抽查了若干村镇的村民，记录了他们的健康状况。则所抽查的这些村镇就是一个随机效应因子.研究随机效应的目的：1.估计随机效应的协方差参数。2.对总体参数作假设检验3.构造总体参数的可信区间。12二.随机效应的概念与识别1213固定效应与随机效应的识别方法:当一个因子(预测变量)对反应变量的效应不易区别是固定效应还是随机效应时,可用可互换性或唯一性规则来作判断.可互换性(exchangeability)判别。

一个随机效应因子的水平是随机地或非系统地选自具有更多水平的总体.观察样本中的水平只是总体中包含的更多水平中的一个随机样本.一个随机效应因子的水平是随机的,在不改变实验的基本性质情况下,当重

复实验时,其水平可能发生改变.固定效应因子:特意选择水平的因子.在不改

变实验基本性质的情况下,当重复实验时,其水平不发生改变.(2)从模型理论上区别:

如果一个效应水平能够合理地假定为代表一种概率分布的话,则该效应就是

随机效应;如果不代表一种概率分布的话,则该效应就是固定效应.

药物疗效试验中的药物品种是特选的,是不能互换的，故为固定效应。而在

药品价格调查中,每类药品选一种作为代表,这时调查的药品名称是可互换的，

故为随机效应。

世界上没有固定效应因子和随机效应因子之分,而是研究者在设计一个实验

时强加的一种不同结构的模型.以便更好地解释客观存在。131314三、

混合效应模型(MixedeffectsModels)混合效应模型是一种线性模型，包含有固定效应和随机效应，用于处理非独立观察资料。又称：

重复测量模型(Repeatedmeasuresmodels)，

多水平模型(Multilevelmodels)，

层次结构模型(Hierarchicalmodels)从统计学归类，混合效应模型包含三种类型的模型:1.随机效应模型(Randomeffectsmodels)。假定除测量误差导致的变异外，还来自具有某种概率分布的随机效应带来的变异，称随机变异。如分析临床多中心试验中的不同中心之间的变异。2.协方差类型模型(Covariancepatternmodels)，直接对重复测量之间的相关结构进行分析，分析效应随时间衰减的特点。3.随机系数模型(Randomcoefficientmodels)，直接分析反应变量在时间轴上的变化率，但容许协变量的效应具有随机性，协变量对反应变量效应的变化率随观察对象而不同，即具有随机变化的特点。以上三种类型的模型可以联合应用。1415

混合效应模型的方差协方差结构其中:Y,为反应变量向量,X为固定效应因子的设计矩阵,β为固定效应参数向量Z为随机效应因子的设计矩阵,γ为随机效应参数向量，γi~N(0,σγ2),e为残差向量.eij~N(0,σe2),cov(γi,e)=0

混合效应模型的参数估计:151616

例子：两种药物治疗效果的随机区组设计模型的矩阵表达其中:Yij,为反应变量,j=1,2代表药物号；，2，。。。，6代表病例号X为（12行,3列）固定效应因子的设计矩阵,β=(μ,β1,β2),为(3行1列)固定效应参数向量Z为(12x6)的随机效应因子的设计矩阵,γi为(6x1)维随机效应参数向量e为(12x!)维残差向量.Patient（i）YijA(j=1)B(j=2)120122262431617429215222162417Mean22.8318.67Treatment：为固定效应，离散化为X矩阵Patient：为随机效应

，离散化为Z矩阵/*---

Model3_2:PROCMIXED---*/PROCMIXEDDATA=intro;

CLASSdrugpatient;

MODELy=drug;

RANDOMpatient;RUN;1617β

=(μ，βj,β2),γ=(γ1,γ2,γ3,γ4,γ5,γ6

),

γi~

(0,σγ2),Var(yij)=

(σγ2+σe2)

1718例2:一项治疗高血压的多中心临床药物试验:三种药物(A,B,C),共有29所医疗中心参与,观察病人总数n=288人.(ofdataset=hypertension)研究目的:在控制治疗前舒张压条件下,分析三种药物的降压效果.四.混合效应模型分析的例子1819变量名:

记录号:patient,n=288,

反应变量:dbp:治疗后舒张期血压,

处理因素:treat:三种药物:A=Carvedilol,B=Nifedipine,C=Atenolol;

控制因素:1.医疗中心:centre,29所医院.2.治疗前舒张期血压:dbp1,连续变量.Obspatientcentretreatdbpdbp11129C86972229C721093

3

5

B109117445A871005529A85105673A1001147

8

3

B801058

9

3

B901009103A10010210113C94105VariableNMeanMinimumMaximum---------------------------------------------------------------dbp28890.246527870.0000000140.0000000dbp1288102.854166792.0000000120.0000000数据集中的前10例病人的记录1920模型A:简单药物效应treatk:k=A,B,C.（固定效应）PROCMIXEDDATA=hypertension;/*MODELA*/CLASScentretreat;

MODELdbp=treat/SOLUTION;RUN;CovarianceParameterEstimatesCovParmEstimateResidual81.5660

FitStatistics-2ResLogLikelihood2076.9AIC(smallerisbetter)2078.9AICC(smallerisbetter)2078.9BIC(smallerisbetter)2082.5

SolutionforFixedEffectsStandardEffecttreatEstimateErrorDFtValuePr>|t|Intercept88.62110.926628595.64<.0001treatA3.00891.29392852.330.0208treatB1.79831.31742851.360.1733treatC0....

Type3TestsofFixedEffectsNumDenEffectDFDFFValuePr>Ftreat22852.730.0670202121模型B:

在模型A的基础上加入基础血压dbp1.（固定效应）PROCMIXEDdata=hypertension;/*modelB/

CLASScentretreat;

MODELdbp=treatdbp1/SOLUTION;RUN;

SolutionforFixedEffectsStandardEffecttreatEstimateErrorDFtValuePr>|t|Intercept58.149011.36662845.12<.0001treatA3.04381.28012842.380.0181treatB2.03871.30632841.560.1197treatC0....dbp10.29540.10982842.690.0076

Type3TestsofFixedEffectsNumDenEffectDFDFFValuePr>Ftreat22842.920.0558dbp112847.230.0076CovarianceParameterEstimatesCovParmEstimateResidual79.8201

FitStatistics-2ResLogLikelihood2072.3AIC(smallerisbetter)2074.3AICC(smallerisbetter)2074.3BIC(smallerisbetter)2078.0CovarianceParameterEstimatesCovParmEstimateResidual81.5660(modelA)

FitStatistics-2ResLogLikelihood2076.9AIC(smallerisbetter)2078.9AICC(smallerisbetter)2078.9BIC(smallerisbetter)2082.52122CovarianceParameterEstimatesCovParmEstimateResidual71.9213

FitStatistics-2ResLogLikelihood1892.6AIC(smallerisbetter)1894.6AICC(smallerisbetter)1894.6BIC(smallerisbetter)1898.2Type3TestsofFixedEffectsNumDenEffectDFDFFValuePr>Fdbp112563.870.0501treat22562.960.0535centre282562.110.0013PROCMIXEDDATA=hypertension;/*modelC*/

CLASScentretreat;

MODELdbp=dbp1treatcentre/SOLUTION;RUN;模型C:

在模型B的基础上加入医疗中心centre.（固定效应）CovarianceParameterEstimatesCovParmEstimateResidual79.8201(modelB)

FitStatistics-2ResLogLikelihood2072.3AIC(smallerisbetter)2074.3AICC(smallerisbetter)2074.3BIC(smallerisbetter)2078.02223SolutionforFixedEffectsStandardEffecttreatcentreEstimateErrorDFtValuePr>|t|Intercept65.579612.94012565.07<.0001dbp10.22300.11332561.970.0501treatA2.99071.23362562.420.0160treatB1.79371.26512561.420.1574treatC0....centre12.02974.49802560.450.6522centre2-4.57805.0511256-0.910.3656Centre

36.13615.19582561.180.2387Centre

40.39234.97442560.080.9372centre

53.54464.98102560.710.4773…⁞…⁞…⁞…⁞…⁞…⁞…Centre

35-1.60426.0736256-0.260.7919Centre

36-1.12724.6359256-0.240.8081centre

37-2.52905.5541256-0.460.6492centre406.27837.37582560.850.3955centre410....2324模型D:

在模型C基础上加入医疗中心与治疗的交互作用:centre*treat.（固定效应）PROCMIXEDDATA=hypertension;/*modelD*/

CLASScentretreat;

MODELdbp=dbp1treatcentrecentre*treat/SOLUTION;

TITLE'MODELD';RUN;CovarianceParameterEstimatesCovParmEstimateResidual69.2614

FitStatistics-2ResLogLikelihood1558.1AIC(smallerisbetter)1560.1AICC(smallerisbetter)1560.1BIC(smallerisbetter)1563.5

Type3TestsofFixedEffectsNumDenEffectDFDFFValuePr>Fdbp112080.990.3198treat22081.240.2905centre282081.980.0038centre*treat482081.200.1884交互效应centre*treat作用项不显著

CovarianceParameterEstimatesCovParmEstimateResidual71.9213(modelC)

FitStatistics-2ResLogLikelihood1892.6AIC(smallerisbetter)1894.6AICC(smallerisbetter)1894.6BIC(smallerisbetter)1898.2242525模型E:

在模型B的基础上,将医疗中心centre作为随机效应.（与固定效应比较）

Type3TestsofFixedEffectsNumDenEffectDFDFFValuePr>Fdbp112566.840.0095treat22563.100.0466PROCMIXEDDATA=hypertension;/*modelE*/CLASScentretreat;

MODELdbp=dbp1treat/SOLUTION;

RANDOMcentre/SOLUTION;RUN;CovarianceParameterEstimatesCovParmEstimatecentre7.8248Residual70.9263FitStatistics-2ResLogLikelihood2056.5CovarianceParameterEstimatesCovParmEstimateResidual79.8201(modelB)

FitStatistics-2ResLogLikelihood2072.3CovarianceParameterEstimatesCovParmEstimate(modelD)Residual69.2614

FitStatistics-2ResLogLikelihood1558.125五种模型的配合结果比较：注：X-无统计学意义，

?-接近α=0.05水平，*-在α=0.05水平上有统计学意义

2627ODSGRAPHICSON;PROCMIXEDDATA=hypertension;

CLASScentretreat;

TITLE‘ResidualPlot

ofModelE';

MODELdbp=treatdbp1/DDFM=KR

RESIDUAL;

RANDOMcentre/SOLUTION;

LSMEANStreat/DIFFPDIFFCL;IDpatientcentretreat;RUN;ODSGRAPHICSOFF;模型E的残差分析:272828例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)。检验目的：分析产品不稳定当来源。2930nosrslotsamplabynosrslotsamplabYnosrslotsamplaby111112006211121999311132007411211980511221988611231982711312000811321998911332007101211199111121219901212131988131221198714122219891512231988161231198517123219831812331989191311200020131220042113132004221321200123132219962413232004251331199926133220002713332002281411199729141219943014131996311421199632142220003314232002341431198735143219903614331995372511201338251220043925132009402521202341252220184225232010432531202044253220234525332015462611203247261220364826132030492621201850262220225126232026522631200953263220105426332011552711198456271219935727131993582721199259272219926027231990612731199662273219936327331987642811199665281219896628131996672821199768282219936928231996702831199071283219897228331992资料表：3031procmixeddata=semiconductormethod=reml;classsourcelotsample;modelcontent_y=source/ddfm=kr;randomlot(source)sample(sourcelot);lsmeanssource/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(反应变量)CovarianceParameterEstimatesCovParmEstimatelot(source)119.89sample(source*lot)35.8657Residual12.5694结论1：批次之间的变异性最大，为样品变异的20倍(=119.89/35.8657),为样品内变异的46倍(=119.86/12.5694)品(lab,m=1,2,3)的化验结果为content_y(反应变量)3132

LeastSquaresMeansStandardEffectsourceEstimateErrorDFtValuePr>|t|source11995.115.77166345.68<.0001source22005.195.77166347.43<.0001DifferencesofLeastSquaresMeansStandardEffectsource_sourceEstimateErrorDFtValuePr>|t|source12-10.08338.16226-1.240.2629结论2：原料之间对产品的质量变异的影响不大Pr>|t|=0.2629)32334:配对病例对照的例子。婴儿卒死综合症（suddeninfantdeathsyndrome，SID）的母亲与同期同一所医院内相同年龄组事件前分娩相同年龄组母亲和事件后分娩相同年龄组母亲各一人作1:2配对调查。研究指标为deprivationscore（depcat,1-7分，分值越高者越严重）。由于各种原因，配对不是完整的。资料的配对情况为:TheFREQProcedureCumulativeCumulativetot