沈敏學(xué) 胡 明△ 曾 娜,2 孫振球
·方法介紹·
雙重差分模型在醫(yī)學(xué)研究中的應(yīng)用*
沈敏學(xué)1胡 明1△曾 娜1,2孫振球1
雙重差分模型(difference-in-differences,DD)由Ashenfelter和Card于1985年對一項(xiàng)干預(yù)研究進(jìn)行評價(jià)時(shí)提出[1],其后該模型在計(jì)量經(jīng)濟(jì)學(xué)領(lǐng)域得到廣泛應(yīng)用[2-4]。DD模型提出的背景是:一項(xiàng)公共政策的實(shí)施無法保障干預(yù)組與對照組的隨機(jī)化,反應(yīng)變量在基線水平上未必可比,因此如果只通過自身前后對照或兩組在橫斷面上的對比,難以得到政策實(shí)施效果的無偏估計(jì)。此外,大規(guī)模的人群調(diào)查難以對個體進(jìn)行隨訪,因此二次抽樣也限制了統(tǒng)計(jì)學(xué)方法的選擇。DD模型適用于醫(yī)學(xué)研究,尤其是公共衛(wèi)生研究[5-10]。本文在介紹DD模型基本原理和形式推廣的基礎(chǔ)上,將其與重復(fù)測量的方差分析和混合模型進(jìn)行比較,并以實(shí)例說明幾種模型的適用條件與區(qū)別。
DD的核心是模型構(gòu)造雙重差分估計(jì)量(DD estimator)。設(shè)分組變量Gi∈{0,1} ,時(shí)間變量Ti∈{0,1} ,DD的基本模型可表示如下式:
Yi=β0+β1Gi+β2Ti+β3(Gi·Ti)+ei
其中,分組和時(shí)間均為啞變量,Gi·Ti是二者的交互項(xiàng),ei為殘差。
值得注意的是,分組變量并不等同干預(yù),因?yàn)樵诨€水平(即Ti=0時(shí)),兩組均未接受干預(yù)。當(dāng)Gi=1且Ti=1,即啞變量Gi·Ti=1時(shí)才指代干預(yù)。不難得出,干預(yù)組和對照組在干預(yù)實(shí)施前后應(yīng)變量的數(shù)學(xué)期望(均數(shù))分別如下:
E[Y|G=1,T=0]=β0+β1
E[Y|G=1,T=1]=β0+β1+β2+β3
E[Y|G=0,T=0]=β0
E[Y|G=0,T=1]=β0+β2
一般線性模型須滿足Gauss-Markov假設(shè),即殘差的均數(shù)為零且獨(dú)立于解釋變量,因而無殘差項(xiàng)。
雙重差分估計(jì)量是橫向和縱向比較的結(jié)合,即干預(yù)組前后差異與對照組前后差異之差,實(shí)際上是時(shí)間和分組交互項(xiàng)的偏回歸系數(shù)β3。在資料滿足線性回歸條件的基礎(chǔ)上,采用最小二乘法(OLS)即可得到β3的無偏估計(jì)[11]。
應(yīng)用DD模型評估干預(yù)效果時(shí),除應(yīng)符合線性回歸的前提外,還應(yīng)滿足3個假設(shè):干預(yù)措施對對照組不產(chǎn)生影響;干預(yù)之外的因素對干預(yù)組和對照組影響相同;干預(yù)組和對照組中觀察單位的某些特征分布穩(wěn)定,不隨時(shí)間變化[12]。
1.DDD模型
自然實(shí)驗(yàn)中,對照組的選取可能影響效應(yīng)估計(jì)的穩(wěn)健性。此處以Gruber給出的例子來說明[13]:某地政府強(qiáng)制要求企業(yè)雇主為其20~40歲的已婚女性員工購買覆蓋生育費(fèi)用的醫(yī)療保險(xiǎn),分析政策對該年齡段已婚女性的工資是否有影響。Gruber選取了3組對照:本地20~40歲未婚男性及40歲以上男性;非政策實(shí)施地20~40歲的已婚女性;非政策實(shí)施地20~40歲未婚男性及40歲以上男性。該例中,如果僅以非政策實(shí)施地的人群為對照,則可能因各地經(jīng)濟(jì)發(fā)展?fàn)顩r不同而對工資造成不等量的影響,對政策效應(yīng)的估計(jì)是有偏的;如果僅以本地男性為對照,則可能存在其他宏觀因素(如國家層面的政策)對男女工資水平有不等量的影響。以T表示時(shí)間變量,G表示分組變量(政策實(shí)施),S表示性別變量,則回歸模型可表示為下式:
Yi=β0+β1Gi+β2Ti+β3(Gi·Ti)+β4Si+β5(Gi·Si)+β6(Ti·Si)+β7(Gi·Ti·Si)+ei
DDD估計(jì)量為干預(yù)組和對照組DD估計(jì)量之差,在滿足線性回歸條件的基礎(chǔ)上,由OLS估計(jì)的DDD統(tǒng)計(jì)量為各變量交互項(xiàng)的偏回歸系數(shù)β7,如下式:
-(E[Y|G=0,T=1,S=1]+E[Y|G=0,T=0,S=1])}
-{(E[Y|G=1,T=1,S=0]+E[Y|G=1,T=0,S=0]
+(E[Y|G=0,T=1,S=0]-E[Y|G=0,T=0,S=0])}=(β3+β7)-β3=β7
2.一般化模型
當(dāng)分組和時(shí)間變量的水平數(shù)超過2時(shí),基本DD或DDD模型可推廣為一般化模型。此外,在自然試驗(yàn)中,干預(yù)組和對照組的非隨機(jī)化分配將導(dǎo)致觀察單位特質(zhì)分布的不均衡,因此在一般化模型中,通常會考慮分組和時(shí)間之外的解釋變量。下式中,向量λT和αG分別為分組和時(shí)間啞變量的偏回歸系數(shù),向量γ為觀察單位變量(如人口學(xué)特征)ZGT的偏回歸系數(shù),向量β為交互項(xiàng)的偏回歸系數(shù),即DD估計(jì)量。
Yi=λΤ·Τ+αG·G+β·XTG+γ·ΖTG+eι
考慮到個體觀測值可能存在非獨(dú)立性[14],一般化模型還可推廣為多水平模型,如下式。
Yi=λΤ·Τ+αG·G+β·XTG+γ·ΖTG+vgι+uigt
其中,vgt為群組水平隨機(jī)誤差,uigt為個體水平隨機(jī)誤差。
3.廣義線性模型
當(dāng)應(yīng)變量不服從正態(tài)分布時(shí),則應(yīng)對其進(jìn)行數(shù)學(xué)變換,或采用廣義線性模型,如Probit模型、Logit模型等[15]。Ai和Norton指出,有學(xué)者誤將廣義線性模型中交互項(xiàng)的偏回歸系數(shù)等同于干預(yù)效果[16]。實(shí)際上,因連接函數(shù)的存在,DD估計(jì)量已演變?yōu)橄率?Φ為Probit函數(shù)):
Yi=Φ(α·G+β·T+γ·G·T)
含義是:“接受了干預(yù)”的干預(yù)組和(G=1,T=1,G×T=1)“假定沒有接受干預(yù)”的干預(yù)組(G=1,T=1,G×T=0)在第2個時(shí)間點(diǎn)上數(shù)學(xué)期望的差值。
將雙重差分模型(此處指線性模型)與醫(yī)學(xué)研究中常用的重復(fù)測量方差分析及混合模型進(jìn)行比較,見表1。
此處列舉兩個例子,第一例為重復(fù)測量數(shù)據(jù),第二例為混合橫斷面數(shù)據(jù)(二次抽樣),分別用雙重差分模型、重復(fù)測量資料的方差分析和混合模型進(jìn)行估計(jì)和結(jié)果的比較。
例1 將20名高血壓病患者隨機(jī)分為兩組,對處理組予以某種治療,對照組不予處理,在治療前后分別測定血壓。以舒張壓為結(jié)局變量,試評價(jià)治療的效果(摘自研究生規(guī)劃教材《醫(yī)學(xué)統(tǒng)計(jì)學(xué)》第四版第12章[11])。
本例是簡單的重復(fù)測量設(shè)計(jì)資料,現(xiàn)用雙重差分模型和重復(fù)測量方差分析分別對其治療效果進(jìn)行估計(jì),結(jié)果見表2和表3。經(jīng)Mauchly檢驗(yàn),該資料滿足“球?qū)ΨQ”假設(shè),即重復(fù)測量誤差的協(xié)方差矩陣為球?qū)ΨQ結(jié)構(gòu),因此兩個模型的估計(jì)都是有效的。重復(fù)測量資料的方差分析顯示,分組和時(shí)間存在交互作用,處理組的降壓效果優(yōu)于對照組,兩組血壓的前后差值分別為-16.0和-4.2,二者之差也剛好等于DD估計(jì)量-11.8,因此兩個模型的結(jié)論一致。二者的區(qū)別在于:重復(fù)測量方差分析中,時(shí)間和分組變量的統(tǒng)計(jì)學(xué)檢驗(yàn)是對其主效應(yīng)的檢驗(yàn)(盡管當(dāng)交互作用存在時(shí),分析主效應(yīng)意義不大);而在DD模型中,“分組”并不等同于“處理”,因?yàn)樵诨€上兩組都未接受治療,其偏回歸系數(shù)是基線時(shí)兩組舒張壓均數(shù)之差。
需注意的是,在實(shí)際研究中,干預(yù)組和對照組并不總是均衡可比的;重復(fù)測量值之間往往具有相關(guān)性,擬合線性模型可能虛增檢驗(yàn)效能,增加I類錯誤的概率;隨機(jī)誤差可能與協(xié)變量有相關(guān)關(guān)系,使用最小二乘法將導(dǎo)致模型的錯誤估計(jì)。因此需要充分考察數(shù)據(jù)的特征和模型的使用前提,并適時(shí)采用加權(quán)最小二乘法、廣義最小二乘法等進(jìn)行參數(shù)估計(jì)。
例2 研究者對西部地區(qū)12所小學(xué)4~6年級學(xué)生進(jìn)行營養(yǎng)知識的健康教育,以學(xué)校為單位進(jìn)行隨機(jī)化分組。干預(yù)實(shí)施前,從干預(yù)組和對照組中共抽取了378名學(xué)生進(jìn)行營養(yǎng)知識測試;干預(yù)一年后,再次隨機(jī)抽取了478名學(xué)生進(jìn)行測試[18]。現(xiàn)以測試總分為評價(jià)指標(biāo),試估計(jì)健康教育的效應(yīng)。
該例沒有對樣本進(jìn)行隨訪,進(jìn)行了二次抽樣,因此無法使用重復(fù)測量資料的方差分析?,F(xiàn)以雙重差分模型和混合模型分別對干預(yù)效應(yīng)進(jìn)行估計(jì)[19-20]。由表4可知,DD模型的雙重差分估計(jì)量為2.90,混合模型估計(jì)的交互項(xiàng)固定效應(yīng)為4.04,高于DD模型的估計(jì)值。實(shí)際上,引入交互項(xiàng)后,混合模型是多水平的DD模型。由隨機(jī)效應(yīng)可知,該數(shù)據(jù)具有層次結(jié)構(gòu),測試分?jǐn)?shù)在學(xué)校水平上存在一定聚集性。盡管兩種方法估計(jì)所得的干預(yù)效果較接近,但由于學(xué)校水平殘差不為零,因此個體水平殘差的分布不再滿足OLS估計(jì)的前提,單水平DD模型可能導(dǎo)致錯誤的推論。
在醫(yī)學(xué)研究尤其是公共衛(wèi)生研究中,未必總能通過隨機(jī)化分配或匹配來保證干預(yù)組和對照組的可比性。雙重差分模型通過構(gòu)造雙重差分估計(jì)量來控制和消除其他協(xié)變量對干預(yù)效果的影響,且不要求個體水平的重復(fù)測量,簡單易行,是評價(jià)自然實(shí)驗(yàn)的良好方法。實(shí)際應(yīng)用時(shí),應(yīng)充分考慮模型協(xié)變量的選取,消除不匹配的因素,以保證效應(yīng)的正確估計(jì);應(yīng)考察模型估計(jì)方法的前提條件和適用性,選取適當(dāng)?shù)膮?shù)、半?yún)?shù)或非參數(shù)估計(jì);當(dāng)資料存在層次結(jié)構(gòu)時(shí),應(yīng)考慮誤差在不同水平上的分布以及變量之間的協(xié)方差和相關(guān)關(guān)系,并可將雙重差分模型與多水平模型結(jié)合應(yīng)用;當(dāng)應(yīng)變量不滿足正態(tài)分布時(shí),還應(yīng)進(jìn)行數(shù)學(xué)轉(zhuǎn)換,或采用廣義線性模型來擬合資料。
[1]Ashenfelter O,Card D.Using the longitudinal structure of earnings to estimate the effect of training programs.The Review of Economics and Statistics,1985,67(4):648-660.
[2]Meyer BD.Natural and quasi-experiments in economics.Journal of Business and Economic Statistics,1995,13(2):151-161.
[3]Bertrand M,Duflo E,Mullainathan S.How much should we trust differences-in-differences estimates.The Quarterly Journal of Economics,2004,119(1):249-275.
[4]Athey S,Imbens GW.Identification and inference in nonlinear difference-in-differences models.Econometrica,2006,74(2):431-497.
[5]Stuart EA,Huskamp HA,Duckworth K,et al.Using propensity scores in difference-in-differences models to estimate the effects of a policy change.Health Serv Outcomes Res Methodol,2014,114(4):166-182.
[6]Zhang X,Wu Q,Liu G,et al.The effect of the National Essential Medicines Policy on health expenditures and service delivery in Chinese township health centres:evidence from a longitudinal study.BMJ Open,2014,4(12):e006471.
[7]Morris S,Hunter RM,Ramsay AI,et al.Impact of centralising acute stroke services in English metropolitan areas on mortality and length of hospital stay:difference-in-differences analysis.BMJ,2014,349:g4757.
[8]Fillmore H,DuBard CA,Ritter GA,et al.Health care savings with the patient-centered medical home:Community Care of North Carolina′s experience.Popul Health Manag,2014,17(3):141-148.
[9]Nolan A.An extension in eligibility for free primary care and avoidable hospitalisations:a natural experiment.Soc Sci Med,2011,73(7):978-985.
[10]Alderman H,Ndiaye B,Linnemayr S,et al.Effectiveness of a community-based intervention to improve nutrition in young children in Senegal:a difference in difference analysis.Public Health Nutr,2009,12(5):667-673.
[11]孫振球,徐勇勇.醫(yī)學(xué)統(tǒng)計(jì)學(xué).第4版.北京:人民衛(wèi)生出版社,2014.
[12]葉芳,王燕.雙重差分模型介紹及其應(yīng)用.中國衛(wèi)生統(tǒng)計(jì),2013,30(1):131-134.
[13]Gruber J.The incidence of mandated maternity benefits.Am Econ Rev,1994,84(3):622-641.
[14]Hansen LP,Sargent TJ.Recursive robust estimation and control without commitment.Journal of Economic Theory,2007,136(1):1-27.
[15]Karaca-Mandic P,Norton EC,Dowd B.Interaction terms in nonlinear models.Health Serv Res,2012,47(1 Pt 1):255-274.
[16]Ai C,Norton EC.Interaction terms in logit and probit models.Economics Letters,2003,80(1):123-129.
[17]Puhani PA.The treatment effect,the cross difference,and the interaction term in nonlinear “difference-in-differences” models.Economics Letters,2012,115(1):85-87.
[18]沈敏學(xué),彭真,胡婧璇,等.中國西部地區(qū)小學(xué)生營養(yǎng)與食品安全教育的整群隨機(jī)試驗(yàn)效果評價(jià).中華流行病學(xué)雜志,2013,34(9):879-883.
[19]張巖波,何大衛(wèi),劉桂芬,等.重復(fù)測量數(shù)據(jù)的混合模型及其MIXED過程實(shí)現(xiàn)——混合線性模型及其SAS軟件實(shí)現(xiàn)(二).中國衛(wèi)生統(tǒng)計(jì),2001,18(5):272-275.
[20]高萌,張強(qiáng),鄧紅,等.方差分析和混合線性模型在重復(fù)測量數(shù)據(jù)中的應(yīng)用探討.現(xiàn)代預(yù)防醫(yī)學(xué),2008,35(7):1221-1223.
(責(zé)任編輯:郭海強(qiáng))
附錄:例題SAS程序
例1
/*例1:重復(fù)測量方差分析*/
data example1a;input t1 t2 group @@;
cards;/*每例患者一條記錄*/
130 114 1
……
134 128 2;
proc glm;class group;model t1 t2=group;repeated time 2
contrast(1)/summary;
run;
/*例1:DD模型*/
data example1b;input BP time group @@;
cards;/*每次測定一條記錄,同一患者可有多條記錄*/
130 1 1
114 2 1
……
134 1 2
128 2 2;
data example1b;set example1b;interact=time * group;
run;
proc reg;model BP=time group interact/std;
run;
例2
/*例2:混合模型*/
data example2;input school time group score @@;
cards;/*每次測定一條記錄*/
1 1 1 14
1 1 1 13
……
12 2 2 18
12 2 2 21;
data example2;set example2;interact=time * group;
run;
proc mixed;class school;model score=time group interact/s;random intercept/sub=school type=vc;
run;
/*例2:DD模型,數(shù)據(jù)格式同上*/
proc reg;model score=time group interact/std;
run;
國家自然科學(xué)基金(81402770)
1.中南大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計(jì)學(xué)系(410078)
2.中南大學(xué)湘雅三醫(yī)院臨床營養(yǎng)科(410013)
△通信作者:胡明,E-mail:huming0129@126.com