混雜效應和隨機效應模型

混雜效應和隨機效應模型當前1頁，總共35頁。一.統(tǒng)計模型的概念統(tǒng)計模型是對資料結構的一種數(shù)學表述.數(shù)量關系的概念化結構.包含兩個元素1.函數(shù)表達式:描述結果變量與解釋變量之間的關系(固定效應).2.誤差表達式:描述結果變量觀察值隨機變異的概率分布(隨機變異).例如:2種藥物(A、B)治療某種疾病的療效分析。用均衡設計，每種藥物治療的病人數(shù)相等，都為n.反應變量：Yij表示生化測定值，i=1,…,為病例編號，J=1,2為藥物編號自變量：藥物種類(A,B)，令Xj=第j種藥物，傳統(tǒng)的統(tǒng)計分析方法（固定效應模型，效應為常數(shù)）為：（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)=γi2+σe2

，2當前2頁，總共35頁。3Patient(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兩種治療藥物在同一病人體內(nèi)實驗,采用區(qū)組隨機化設計方案(即用藥先后順序是隨機化的)，對每種藥物處理后的反應變量進行測定.用6例病人.結果如下表.構造三種模型:完全隨機設計模型：不考慮區(qū)組（病人）效應:Yij=μ+βj

+eij

,βj為藥物效應隨機化區(qū)組設計模型：考慮區(qū)組（病人）效應：Yij=μ+βj

+αi+eij

隨機效應模型：病人是從病人總體中隨機的，也存在隨機誤差，統(tǒng)計學中用病人間的方差來衡量這種隨機誤差.。

Yij=μ+βj

+(γi)+eij

==(μ

+βj+

(γj+eij

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

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

(τγ2+σe2)

在此簡單情況下，(3)與(2)等價，但解釋不同。在有缺失值情況下的結果不同。當前3頁，總共35頁。4模型一：完全隨機設計模型:

βJ：第J種藥物效應從上表估計模型參數(shù):

μ=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.75當前4頁，總共35頁。5完全隨機設計模型的PROCANOVA計算結果:

SumofSourceDFSquaresMeanSquareFValuePr>FModel152.083333352.08333332.680.1325Error10194.166666719.4166667CorrectedTotal11246.250000Meanswiththesameletterarenotsignificantlydifferent.SNKGroupingMeanNdrugA22.8336AAA18.6676BPROCANOVADATA=example_1;CLASSdrug;MODELy=drug;MEANSdrug/SNKALPHA=0.05;run;當前5頁，總共35頁。66完全隨機設計模型的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;當前6頁，總共35頁。7模型二:隨機區(qū)組設計模型(考慮病人效應α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)當前7頁，總共35頁。8模型三:病人為隨機效應的模型:因此,對同一病人的不同觀察之間是相關的,具有協(xié)方差σγ2,包含在總方差Var(yij)=σe2+σγ2內(nèi),σγ2和σe2

都稱為方差分量.但特別指σγ2。當前8頁，總共35頁。9PROCMIXEDDATA=example_1;/*Model3:RandomeffectsmodelbyusingPROCMIXED*/

CLASSdrugpatient;

MODELy=drug;

RANDOMpatient/S;RUN;用SAS中的PROCMIXED計算結果:CovParmEstimatepatient11.5333(用PROCGLM的RANDOM語句得不到此方差分量)Residual7.8833(組內(nèi)相關系數(shù)ICC=11.53/(11.53+7.88)=0.59)

Type3TestsofFixedEffectsNumDenEffectDFDFFValuePr>Fdrug1

5

6.61

0.0500在本例中,對drug的檢驗,用PROCMIXED的計算結果與用PROCGLM(2)的計算結果同（F=6.61）,即規(guī)定病人是固定效應,還是隨機效應,對處理效應的檢驗結果沒有影響（這是由于方差的性質決定的，即觀察值的方差與中心化值的方差相等）.但如果有缺失值時,其結果不同.在本例的模型三中,假定病人具有隨機效應.病人來自一個具有均值為0,方差為σα2的正態(tài)分布總體.因此它們的期望值為0,但每個病人彼此不同。每個病人都具有相同期望值的假定與直觀不符.須根據(jù)每例病人的觀察值,確定其在正態(tài)分布中的一個位點.這一預報值的可信區(qū)間較固定效應的可信區(qū)間要窄,在統(tǒng)計學上稱為收縮”shrunken”估計.當前9頁，總共35頁。10在本例的模型三中,假定病人具有隨機效應.即規(guī)定病人來自一個具有均值為0,方差為σα2的正態(tài)分布總體.因此它們的期望值為0。但每個病人彼此不同。每個病人都具有同一期望值的假定與直觀不符.須根據(jù)每例病人的觀察值,確定其在正態(tài)分布中的一個位點.這一預報值的可信區(qū)間較固定效應的可信區(qū)間要窄,在統(tǒng)計學上稱為收縮”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從上表可見，隨機效應模型的預報值更接近觀察值。當前10頁，總共35頁。隨機效應:在一項研究中，如果進入研究的因子的水平數(shù)只是其總體中的所有水平數(shù)的一個隨機代表時,該因子的效應為隨機效應。對應于該因子的總體中各水平的效應就構成了一個概率分布總體。樣本中的各水平是來自總體中更多水平的一個隨機樣本.如:一個城市有很多學校,為了解學生體質,抽查了部分學校的學生體質，則所抽查的學校就是一個隨機效應因子.一個城市有很多醫(yī)院,為了解醫(yī)療質量,抽查了若干醫(yī)院的出院病人記錄進行醫(yī)療質量分析。則所抽查的醫(yī)院就是一個隨機效應因子.一條大河有許多支流的入口，為評價河水的亞硝酸鹽濃度，只能抽查一小部分支流入口處的水樣作檢驗，所抽查的支流入口處的水樣就是一個隨機效應因子.一個縣市有很多村鎮(zhèn),為了解村民健康狀況,隨機抽查了若干村鎮(zhèn)的村民，記錄了他們的健康狀況。則所抽查的這些村鎮(zhèn)就是一個隨機效應因子.研究隨機效應的目的：1.估計隨機效應的協(xié)方差參數(shù)。2.對總體參數(shù)作假設檢驗3.構造總體參數(shù)的可信區(qū)間。11二.隨機效應的概念與識別當前11頁，總共35頁。12固定效應與隨機效應的識別方法:當一個因子(預測變量)對反應變量的效應不易區(qū)別是固定效應還是隨機效應時,可用可互換性或唯一性規(guī)則來作判斷.可互換性(exchangeability)判別。

一個隨機效應因子的水平是隨機地或非系統(tǒng)地選自具有更多水平的總體.觀察樣本中的水平只是總體中包含的更多水平中的一個隨機樣本.一個隨機效應因子的水平是隨機的,在不改變實驗的基本性質情況下,當重

復實驗時,其水平可能發(fā)生改變.固定效應因子:特意選擇水平的因子.在不改

變實驗基本性質的情況下,當重復實驗時,其水平不發(fā)生改變.(2)從模型理論上區(qū)別:

如果一個效應水平能夠合理地假定為代表一種概率分布的話,則該效應就是

隨機效應;如果不代表一種概率分布的話,則該效應就是固定效應.

藥物療效試驗中的藥物品種是特選的,是不能互換的，故為固定效應。而在

藥品價格調查中,每類藥品選一種作為代表,這時調查的藥品名稱是可互換的，

故為隨機效應。

世界上沒有固定效應因子和隨機效應因子之分,而是研究者在設計一個實驗

時強加的一種不同結構的模型.以便更好地解釋客觀存在。12當前12頁，總共35頁。13三、

混合效應模型(MixedeffectsModels)混合效應模型是一種線性模型，包含有固定效應和隨機效應，用于處理非獨立觀察資料。又稱：

重復測量模型(Repeatedmeasuresmodels)，

多水平模型(Multilevelmodels)，

層次結構模型(Hierarchicalmodels)從統(tǒng)計學歸類，混合效應模型包含三種類型的模型:1.隨機效應模型(Randomeffectsmodels)。假定除測量誤差導致的變異外，還來自具有某種概率分布的隨機效應帶來的變異，稱隨機變異。如分析臨床多中心試驗中的不同中心之間的變異。2.協(xié)方差類型模型(Covariancepatternmodels)，直接對重復測量之間的相關結構進行分析，分析效應隨時間衰減的特點。3.隨機系數(shù)模型(Randomcoefficientmodels)，直接分析反應變量在時間軸上的變化率，但容許協(xié)變量的效應具有隨機性，協(xié)變量對反應變量效應的變化率隨觀察對象而不同，即具有隨機變化的特點。以上三種類型的模型可以聯(lián)合應用。當前13頁，總共35頁。14

混合效應模型的方差協(xié)方差結構其中:Y,為反應變量向量,X為固定效應因子的設計矩陣,β為固定效應參數(shù)向量Z為隨機效應因子的設計矩陣,γ為隨機效應參數(shù)向量，γi~N(0,σγ2),e為殘差向量.eij~N(0,σe2),cov(γi,e)=0

混合效應模型的參數(shù)估計:當前14頁，總共35頁。1515

例子：兩種藥物治療效果的隨機區(qū)組設計模型的矩陣表達其中:Yij,為反應變量,j=1,2代表藥物號；，2，。。。，6代表病例號X為（12行,3列）固定效應因子的設計矩陣,β=(μ,β1,β2),為(3行1列)固定效應參數(shù)向量Z為(12x6)的隨機效應因子的設計矩陣,γi為(6x1)維隨機效應參數(shù)向量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;當前15頁，總共35頁。16Patientdrug反應變量μβ1β2γ1γ2γ3γ4γ5γ6ijYX0x1x2z1z2z3z4z5z6參數(shù)向量1122012111001110000000000μ2122624111001001100000000β1β23121617111001000011000000γ1γ24122921111001000000110000γ3γ4

5122221111001000000001100γ5γ66122417111001000000000011β

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

),

γi~

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

(σγ2+σe2)

當前16頁，總共35頁。17例2:一項治療高血壓的多中心臨床藥物試驗:三種藥物(A,B,C),共有29所醫(yī)療中心參與,觀察病人總數(shù)n=288人.(filenameofdataset=hypertension)研究目的:在控制治療前舒張壓條件下,分析三種藥物的降壓效果.29所醫(yī)療中心的每種藥物治療病人數(shù)及總病人數(shù)MedicalCenterTreatAllMedicalCenterTreatAllABCABC113141239231.2323431024..1133328253227444412263431054521127.11262125291.237666183012258222631121212369..1132211411444123521141243411369682313112437312614888244011.21544311412114182226All1009395288四.混合效應模型分析的例子當前17頁，總共35頁。18變量名:

記錄號:patient,n=288,

反應變量:dbp:治療后舒張期血壓,

處理因素:treat:三種藥物:A=Carvedilol,B=Nifedipine,C=Atenolol;

控制因素:1.醫(yī)療中心:centre,29所醫(yī)院.2.治療前舒張期血壓:dbp1,連續(xù)變量.Obspatientcentretreatdbpdbp11129C86972229C721093

3

5

B109117445A871005529A85105673A1001147

8

3

B801058

9

3

B901009103A10010210113C94105VariableNMeanMinimumMaximum---------------------------------------------------------------dbp28890.246527870.0000000140.0000000dbp1288102.854166792.0000000120.0000000數(shù)據(jù)集中的前10例病人的記錄當前18頁，總共35頁。19模型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.0670當前19頁，總共35頁。2020模型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.5當前20頁，總共35頁。21CovarianceParameterEstimatesCovParmEstimateResidual71.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的基礎上加入醫(yī)療中心centre.（固定效應）CovarianceParameterEstimatesCovParmEstimateResidual79.8201(modelB)

FitStatistics-2ResLogLikelihood2072.3AIC(smallerisbetter)2074.3AICC(smallerisbetter)2074.3BIC(smallerisbetter)2078.0當前21頁，總共35頁。22SolutionforFixedEffectsStandardEffecttreatcentreEstimateErrorDFtValuePr>|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....當前22頁，總共35頁。23模型D:

在模型C基礎上加入醫(yī)療中心與治療的交互作用: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.2當前23頁，總共35頁。2424模型E:

在模型B的基礎上,將醫(yī)療中心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.1當前24頁，總共35頁。ModelFixedeffectsRandomeffectsLogLikelihoodAIC(smallerisbetter)CovparmResparmATreat*2076.92078.981.56BTreat*

dbp1*2072.32074.379.82Ctreat?dbp1?centre*1892.61894.671.92DtreatX

dbp1X

centre*centre*treatX1558.11560.169.26ETreat*dbp1*centre2056.52060.57.824870.93五種模型的配合結果比較：注：X-無統(tǒng)計學意義，

?-接近α=0.05水平，*-在α=0.05水平上有統(tǒng)計學意義

當前25頁，總共35頁。26ODSGRAPHICSON;PROCMIXEDDATA=hypertension;

CLASScentretreat;

TITLE‘ResidualPlot

ofModelE';

MODELdbp=treatdbp1/DDFM=KR

RESIDUAL;

RANDOMcentre/SOLUTION;

LSMEANStreat/DIFFPDIFFCL;IDpatientcentretreat;RUN;ODSGRAPHICSOFF;模型E的殘差分析:當前26頁，總共35頁。27當前27頁，總共35頁。例3:三水平混合效應模型的例子。一個制藥廠為了解生產(chǎn)的穩(wěn)定性，進行了抽樣研究。兩種原料，從原料(1)中抽取4批產(chǎn)品，從原料(2)中抽取4批產(chǎn)品，每批產(chǎn)品中隨機抽取3分樣品，再從每份樣品中隨機取出3分檢驗樣品。總檢驗數(shù)為：2x4x3x3=72個檢驗品。抽樣過程為:原料(source,i=1,2)↓每種原料抽取4批產(chǎn)品(lot,j=1,2,3,4)↓每批產(chǎn)品中抽取3分樣品(samp,k=1,2,3)↓再從每份樣品中取出3分檢驗品(lab,m=1,2,3)檢驗項目：有效物質含量(y)。檢驗目的：分析產(chǎn)品不穩(wěn)定當來源。當前28頁，總共35頁。29nosrslotsamplabynosrslotsamplabYnosrslotsamplaby111112006211121999311132007411211980511221988611231982711312000811321998911332007101211199111121219901212131988131221198714122219891512231988161231198517123219831812331989191311200020131220042113132004221321200123132219962413232004251331199926133220002713332002281411199729141219943014131996311421199632142220003314232002341431198735143219903614331995372511201338251220043925132009402521202341252220184225232010432531202044253220234525332015462611203247261220364826132030492621201850262220225126232026522631200953263220105426332011552711198456271219935727131993582721199259272219926027231990612731199662273219936327331987642811199665281219896628131996672821199768282219936928231996702831199071283219897228331992資料表：當前29頁，總共35頁。30procmixeddata=semiconductormethod=reml;classsourcelotsample;modelcontent_y=source/ddfm=kr;randomlot(source)sample(sourcelot);lsmeanssource/diff;run;原料(source,i=1,2)為固定效應每種原料抽取4批次的產(chǎn)品(lot,j=1,2,3,4)為隨機效應(_嵌套于原料)每批產(chǎn)品中抽取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),為樣品內(nèi)變異的46倍(=119.86/12.5694)品(lab,m=1,2,3)的化驗結果為content_y(反應變量)當前30頁，總共35頁。31

LeastSquaresMeansStandardEffectsourceEstimateErrorDFtValuePr>|t|source11995.115.77166345.68<.0001source22005.195.77166347.43<.0001DifferencesofLeastSquaresMeansStandardEffectsource_sourceEstimateErrorDFtValuePr>|t|source12-10.08338.16226-1.240.2629結論2：原料之間對產(chǎn)品的質量變異的影響不大Pr>|t|=0.2629)當前31頁，總共35頁。324:配對病例對照的例子。嬰兒卒死綜合癥（suddeninfantdeathsyndrome，SID）的母親與同期同一所醫(yī)院內(nèi)相同年齡組事件前分娩相同年齡組母親和事件后分娩相同年齡組母親各一人作1:2配對調查。研究指標為deprivationscore（depcat,1-7分，分值越高者越嚴重）。由于各種原因，配對不是完整的