回归分析面板在社会科学中的应用

第七讲史分析（Event 第一节史分析概计分析方法被称为史分析（EventHistoryysis）或生存分析（survivalysis即观察该从某一观测时点起直到的时间。若，我们记成发生(=1)，若生存，则记成未发生(=0)。当然，生存分析并不仅仅局限于医学数据。只要“”的定义广义化，如定义一个公司的为“”，则公司的生存时间就是该公司从成立到所需要的时间，即记为发生(=1)；到的时间，即公司记为发生(=1)，这样处理就会产生生存数据。典型的连续型数据。若我们关心哪些因素(如治疗方案、、等)能够影响则，似乎OLS就可以，但远没有那么简单!(甚至超过研究者的观察年限)，我们不可能密切每个并确切知道他8年，这就出现了删截数据，即与数据向上删截(right现在引入具体的案例来说明：Krall,UthoffHarley(1975)收集并分析了一个某17表变量变量名 解释意 从确诊到的生存时间(单位：月 生存状态(0=生存；1= 确诊时血蛋白含①在社会统计学派中，一般将社会 在Stata中可以简单地读入数据如下：40140107014014010701070108010801139010900090109012101001090121010010111121122在以上数据中，根据患者的生存状态(Vstatus)，可以将其分为两大类一类是，其生存状态为“”(即Vstatus=1)，如第20个，他的生存时间确实只有2个月(因为Time=2)；另一类是，其生存状态为“生存”(即Vstatus=0)，如第1个。这说明什么指在该研究结束之前，还没观测到第1个的确切时间。因此，我们从数据中仅仅知道他的生存时4个月(Time=4)，但是到底大多少我对于这种数据或在史分析中，我们感的不仅仅是某一是否发生(事究方法称之为生存分析(因此类方法最常用于研究)或社会史分析。所谓史分析，是指分析一系列研究的发生及其发生时间的统计方法，包有坏事，但按照习惯称之为“失败(failure)或”。于是，失败或死models学中被称为“失效时间分析”（failuretimeysis史分析对研究社会科学和自然科学中的许多(events)都很有帮助。如疾数据（longitudinaldata）中的发生。被定义为在时间序列中发生了性质上的变化。称为“同期群研究/面板数据研究”(panelstudies。data，件。而且，也有可能在提供一些时变变量（time-dependentcovariates）的信息方此外，受访者样本也可能是一个有偏误的样本（biasedsample。它可能只包括covariate什么时候发生。可以想象，一些自变量并不随着时间变化（比如，，另一些可能会在这段时间中发生变化（比如受雇状况和状况。Logit模型（被逮捕或者未被逮捕，但是这样就些在52都没有被抓的人，在被逮捕的可能性上是不一样的。然后对这个因变量估计一个线性回回归模型。但是，那些在52周之内都没有被models）以及Cox模型（Coxproportionatemodels。第二节史分析的数据结构及建史分析中的每个个案数据都至少要包括两个特点（1）标志是否发生的变量（2）发生的长度(duration)或是否被删截在给某一赋值(assignnumber)的时候，我们实际上已经选择了相应的尺（origin点上的不同会造成估计结果上的差别。时间的起点可以基于：、日历时间time险集（Riskset。用-stset-命令设置为生存时间(survival-time)。Stset只需运行一次，且数据存盘即针对同一对象的单一记录.stsettimevar,failure***变量timevar(如week或year)即提供了或是到某一特定(failure)发生前或到观测结针对同一观测对象的多条记录.stsettimevar,failure(failvar)id(patient)enter(time***patient是个识别码。在此数据中，同一观测对象可以有多于一条的记录，但相应记录案例分析：数据（recidivism数据说明：432个从放出的罪犯在放出后被了一年。Week=thenumberofweekssincearrest=firstarrestsincereleaseage=agewexp=hadworkexperience(1=yes,0=otherwise)mar=marriedornot(1=yes,0=otherwise) ***source:P.H.RossiandR.A.BerkandK.J.Lenihan(1980).Money,WorkandCrime:SomeExperimentalResults,AcademicPress:NewYork.***RecidivismData:Thisdatacomprises432maleinmatesreleasedfromMarylandstateprisonsintheearly1970's.Datesofarrest,whichconstitutesanevent,wererecordedforeachinmateupto1yearfollowingrelease.Amajorgoaloftheysiswastodetermineifthosewhoreceivedfinancialaidhadlowerarrestratesthanthosewhodidnot.Moreover,theeffectofhavingfinancialaidislikelytobemoreimportantimmediayfollowingarrestandlessimportantastimepasses.Adataframecontaining:weeksurvivalarrestcensoring financialaid,0=no1=yes ageinyears race(1=black, fulltimeworkexperience(0=no,1=yes) married(0=no,1=yes) releasedonparole(0=no,1=yes) numberofpriorconvictions.listweekarrestagewexpmarin②(n=1445)参见：use |weekarrestagewexpmar | 0| 0| 0| 1| 0 | 0| 1| 0| 0| 0 | 0| 0| 0| 0| 0 | 0| 0| 0| 0| 0 | 0022.011|23.110|24.010|25.010|.stsetweek,failure(arrest==1)failureevent: arrest==1obs.timeinterval:(0,week]orbefore:total0obs.remaining, failuresinsinglerecord/singlefailure ysistimeatrisk,atrisk t=earliestobserved t=lastobserved t=该案 第三节生存函数和风险函数FunctionS(t)=Pr{T>t}=1-Pr{T≤t}=1-前“”的可能性（累积分布函数。0≤S(t)≤1,S(t)不会随着时间而变大（一）Kan-Meier估KM估计Sˆtt的个案在样本中的比例。比如，如果25%的个案的生存时间超过了10，那么③timeatrisk:Foraparticularsubject,timeatriskisthetimeinyearsbetweenentryintoanysisandexitfromanysis.Thisisthetimeduringwhichasubjectisatriskofhavinganevent(whethersecondcancer,death,etc.)asdefinedbytheysis.Total timeatriskforagivensubcategoryiscalculatedbyadding timeatriskcountedinthatsubcategoryforeachsubjectintheysis.医学研究上定义 -time),是为了计算在研究中所观察的样本,其未罹患疾病,或未发生之时间的总和,因此人-年( -years)是计算每一位参与研究的,对研究的时间作出了不同的贡献,研究者以“年”(years)来做为测量,例如人-年可以下列式子计算：10-yearsatrisk=1peoplefor10=2peoplefor5=10peoplefor1=100peoplefor1/10在单个“右删截”（rightcensoring）的情况下（所有的个案删截于C），对于所的t≤c，Sˆt仍然等于样本中生存时间超过时点t的个案所占的比例。对于t>c，生存函数无法定义（如数据中一些个案的情况。KM估计一般被定义为：Sˆ(t)niti .stsfailure_d:arrest==1ysistime_t:week Time Fail [95%Conf.1102103104105106107108509201020201030202030302050201010403030202020102020204030401020402040202040102050304 (weeksince我们还可以分组来比较生存函数（结婚与否、工作与否.stsgraph,ylabel(0(.2)1)xlabel(0(5)52)l1title(ProportionUnarrest)b2title(weeksincerelease)by(mar)stsgraph,ylabel(0(.2)1)xlabel(0(5)52)l1title(ProportionUnarrest)b2title(weeksincerelease)by(wexp).ststestfailure_d:arrestysistime_t:weekLog-ranktestforequalityofsurvivor mar|observed +0| 1 +Total| chi2(1)=Pr>chi2=.ststesfailure_d:arrestysistime_t:weekLog-ranktestforequalityofsurvivor|Events wexp|observedexpected+0| 1| +Total| chi2(1)=Pr>chi2=h(t)ff(tlimPr(tTtt 换句话说，风险率（hazard/risk）等于某段时间内“”的非条件概率比上从该Function时段（t,t+Δt）的瞬间风险率（instantaneousrisk）可以被定义为：h(t)limPr(tTtt|TTt 第四节史分析的参数模型（Parametric时还假设，对于第i个患者，我们还知道与其生存时间ti相关的解释性变量 ,xip，那么我们应该怎么量化生存时间tixi首先注意到，生存时间ti是一个非负的随量，而一般的回归模型要求logtx' 的情况下的对数生存时间，记为logi0。因此，以上的线性模型可以变形为logtx'logtttexp(x' 这说明，若没有协变量的影响，该的实际生存时间应该为ti0(BaselineSurvival'i对于加 模型，我们还有一个重要的问题没有回答，即我们对残差项(eilogti0)作elog N(0,2 设为Weibull分布，即： e1exp

ei 其中，,让人欣喜的是，Weibull分布非常灵活，若我们定义1Weibull分布则变化为另一个重要的生存分布，即指数分布(exponentialdistribution)： expei/一、Exponential模型logh(t)=β0+model着时间t变化。二、Gompetz模型：h(t)=exp(α+βX)exp(γtlogh(t)=β0+β0x1+…+βkxkexponential三、Weibull模型：h(t)=ptp-1-logh(t)=β0+β0x1+…+βkxk成了exponential模型。其中，四、其他log-normallog-logistic模modelsgoodness-of-fit。【各参数模型STATA实现】例子:数据（RecidivismDataexponential模.stregagewexpmar,dist(exp)nolognohrfailure_d:arrestysistime_t:Exponentialregression--logrelative-hazardNo.ofsubjects=432 Numberofobs=432No.offailures=114Timeatrisk=LRchi2(3)=Loglikelihood=- Prob>chi2=_t Std. P>|z|[95%Conf.+

Cox模型相比较。age|wexp|mar|_cons

-2.550.011- ---- -1.46---- -1.040.299-1.136397 -7.740.000-4.592454-.stcurve,.stregagewexpmar,dist(exp)nologfailure_d:arrestysistime_t:Exponentialregression--logrelative-hazardNo.ofsubjects=432 Numberofobs=432No.offailures=114Timeatrisk=LRchi2(3)= _t|Haz.RatioStd.Err. P>|z|[95%Conf.+age| . -2.55 wexp| mar|.- 表中的风险率比（Hazardratio）是回归系数的指数形式(如0.95=exp-0.055)，.estatGompertz模.stregagewexpmar,dist(gompertz)nologfailure_d:arrestysistime_t:Gompertzregression--logrelative-hazardNo.ofsubjects=432 Numberofobs=432No.offailures=114Timeatrisk=LRchi2(3)= _t|Haz.RatioStd. P>|z|[95%Conf.+age..-0.009 wexp..-0.131 mar-0.297 + | 3.03 .stcurve, GompertzWeibull模.stregagewexpmar,dist(weibull)nologfailure_d:arrestysistime_t:Weibullregression--logrelative-hazardNo.ofsubjects= Numberofobs=No.offailures=114Timeatrisk=19809LRchi2(3)=Loglikelihood=- Prob>chi2=_t|Haz. Std.Err. [95%Conf.+age|.-0.010 wexp|.-0.129 mar|+.-0.295 /ln_p|+. .147469p|6 .1588971/p| .stregagewexpmar,dist(weibull)nolog.stcurve,注：pWeibullLognorm&Log-Log模.stregagewexpmar,dist(lognorm)nologfailure_d:arrest==1ysistime_t:Log-normalregression--acceleratedfailure-timeNo.ofsubjects= Numberofobs=No.offailures=114Timeatrisk=19809LRchi2(3)= _t| Std. P>|z|[95%Conf.+age|. 0.075-. wexp|. 0.027. mar| 0.127- _cons|3.868942 10.710.0003.160807+/ln_sig| .130644+sigma|1.324466 1.139562.stregagewexpmar,dist(loglog)nologfailure_d:arrest==1ysistime_t:Log-logisticregression--acceleratedfailure-timeNo.ofsubjects=432 Numberofobs=432No.offailures=114Timeatrisk=LRchi2(3)=Loglikelihood=- Prob>chi2=|+Std.z[95%|.006168|.-.|.-_cons.10.26 2.883079....../ln_gam|- -4.67 - -+ | 多个模型中所估计的系数很相近如果估计结果报告的是风险率比，对系数的理解就很简单。比如，在exponential模型中，每增加一岁，在其他变量不变的情况下，被逮捕的风险率5.3%1-0.947)；同理，控制其他变量之后，结了婚的人被逮捕的风险率是没结婚的人的67.4%。很难说哪个模型最好，其选择基于我们对发生情况的假设Gompertz模型在社会科学中很少用到（虽然它最早在率研究中出现Weibullexponential引入Cox风险模型。第五节Cox（ProportionalHazardsModel,Coxh(t)limP(ttitt|tit)ft f(tS(t的负导数，即生存时间ti的密度函数，而h(t就是人们常说的0S(t)expt0方式，建立风险函数h(t)和解释性变量t之间的关系，就可以获得一个关系生存数著名统计学家DavidCox(1972提出了革命性的Cox等比例风险模型（ProportionalHazardsModel-CoxModelCox模型，以对参数模型(如Exponential、WeibullGompertz模型)等进行一般化(generalized)处理。该模型h(t)h0(t)exp了exp'x倍。因此，我们称此模型为等比例风险模型（ProportionalHazardsModel,CoxMoe理论上人们可以验证，通过恰当的选取基准风险函数h0(t)的函数形式，Cox模型就包含了前形式更加灵活。而且，Cox模型的最大优点是可以在不对基准风险函数h0(t作任何假设的情况下，能获得对回归系数的可靠估计及相关的统计推断方法，在此意义说，可以说Cox模型是生存分析技术史上的性突破，至今依然极具生命力。假设ti是一个观测到的时间，也就是说，有那么一个，他的生存时间是ti，并且没首先我们知道他至少存活了ti时间，并且“立刻”了。但是，存活时间至少为ti有多大呢？在Cox模型下，此概率为：0 ih(t)exp'x exp0 ih(t)exp'xexp'x0jR(ti

jR(ti其中，R(t)是这样一个集合，它包含了所有存活时间至少为ti的，而不再依赖基准风险函数的分布，即误差项的分布情况可以function时间点汇一起，就构成了似然函数，或称之为偏似然函数(partiallikelihood)，形式④这也就是为什么说“Cox模型对参数模型(Exponential、WeibullGompertz模型)等进行一般。若=h0(t=αlogt，它就转化Weibull exp'xti expxj'jR(tih0。其中，h0。一、Cox模型的STATA数据：数据（RecidivismData.stcoxagewexpmar,nologfailure_d:arrestysistime_t:Coxregression--BreslowmethodforNo.ofsubjects=432 Numberofobs=432No.offailures=114Timeatrisk=LRchi2(3)=Loglikelihood=- Prob>chi2=_t|Haz.RatioStd. P>|z|[95%Conf.+age| -0.009 wexp| -0.130 mar| -0.297 Cox模型中回归系数的理解与普通模型相同。如，每增加一岁，在其他变量当然，我们也可以在这一命令中添加选项nohr而得到模型中的回归系数：⑤但部分提及Cox模型的假定为：风险比不随时间变化(参见LawrenceC.Hamilton.StatisticswithStata.Schoenfeld残差的检验)。有关命令与选项的说明，参见helpstcox。.stcoxagewexpmar,nolognohrfailure_d:arrestysistime_t:Coxregression--BreslowmethodforNo.ofsubjects=432 Numberofobs=432No.offailures=114Timeatrisk=LRchi2(3)=_t|+ageStd.zP>|z|[95%Conf._t|+ageStd.zP>|z|[95%Conf..---wexp|mar

- -- -1.137891--例中没有出现)。假如Cox模型所拟合的mar(在婚状态)的系数是负的，而加速模型Weibull模型的拟合结果是正的。这是为什么呢？加速模型中mar系数为正说明在婚状态的人存活时间长(即没再的持续时间长)，这正好表明此类人群的风险函数小，因此在Cox模型中该变量系数的符号为负。它已经被吸收到了基准风险函数h0(t)中去了。附：Cox模型与加速模型系数符号相反的案例.use"D:\教学篇\中财教学\备课\回归分析在社会科学中的应用\课件\第七讲_史分析\第七讲社会史Stata实现\dataofpatient.dta",clear.stcoxhgbetagelogwbclogpbmproteinscalc,failure_d:vstatus==1ysistime_t:timeCoxregression--BreslowmethodforNo.ofsubjects Numberof No.offailuresTimeat LR=Loglikelihood-Prob>=_t|Haz. Std. [95%Conf.+hgb..- et..- age..- logwbc logpbm. scalc...stcoxhgbetagelogwbclogpbmproteinscalc,nolognohrfailure_d:vstatus==1ysistime_t:Coxregression--BreslowmethodforNo.of=Numberof=No.of=Timeat=LR=Log=-Prob>=_t+Std.z[95%hgb- .--.et-.--.age-.--.logwbc.-logpbm - -.scalc -.streg etagelogwbclogpbmproteinscalc,dis(weibull)nologfailure_d:vstatus==1ysistime_t:timeWeibullregression--acceleratedfailure-timeNo.of=Numberof=No.of=Timeat=LR=Log=-Prob>=_t Std. [95%Conf.+hgb|et|age|logwbc|

..-

- . - . - - - logpbm|-.protein|-.scalc|-.

- - - - - - _cons +/ln_p - +p 1/p 请学员试着解释符号不同的原因。二、 模型的变量选择与模型筛.estat用偏似然估计法(partialumlikelihoodmethod)进行估计的，其如下：PL

e n[ n Yexj其中，Yij=1如果tj≥ti，Yij=0如果tj<ti；δi=1如果发生，否则为0（被删注意：偏似然函数(partiallikelihood)的值取决于时间（eventtime）的序列，三、Cox模型的不足续的时间，会导致系数估计有偏倚。中的Efron方法等，能校正连续数据的系数和标准误的估计。命令如下：.stcoxagewexpmar,efronnologfailure_d:arrestysistime_t:Coxregression--EfronmethodforNo.ofsubjects=432 Numberofobs=432No.offailures=114Timeatrisk=LRchi2(3)=Loglikelihood=- Prob>chi2=_t|Haz.RatioStd. P>|z|[95%Conf.+age| -0.009 wexp| -0.130 mar| -0.297 四、 模型等比例风险假设的检Cox比例风险模型假设任两个间的风险率之比不随时间变化(proportionalStatastphplot,byvar)var的每个类别画出-ln(-ln生存))ln(时间)Cox.stphplot,by(wexp)adjust(age)failure_d:arrestysistime_t:若比例风险假设 hazards，PH）不成立，则意味着自变量与时间event第六节离散时间风险模型（Discrete-timeHazard值得注意的是，发生的时间出现（tiedeventtimes）的情况很常见，如同时发生的可能性就大为增加。(1)只发生在离散时点上（比如，的（比如，同期群中的职业转换因此，将这类数据视作是离散，比视为连续更有意义理方法是将某一发生的条件概率进行logit转换，并将其视作自变量的函数log(Pit) 1 i

k其中，Pit是某一在时点t之前没发生的条件下，第i个人在时点t发生该事为了估计该模型，我们的数据通常会转化为“人-记录”（-records）数据（即换句话说，在时间按年计算的情况下，假设A在数据中被观察了3年，3案例：数据(Recidivism.gen.listin |weekarrestagewexpmarid 1 2 3 4 5 6 7 8 9| 10 .expand(19377observations.sort.listidweekarrestagewexp |idweekarrestage mar | 0| 0| 0| 0| 0 | 0| 0| 0| 0| 0 | 0| 0| 0| 0| 0 | 0| 0| 0| 0| 0 | 0| 0| 0| 0| 0 | 0| 0| 1| 1| 1 | 1| 1| 1| 1| 1 | 1| 1| 1| 1| 1 第三步：对每个人标志“spell.listidweekarrestagewexpmartin||||t|||11001||11002||11003|| 4| 5 | 6| 7| 8| 9| 10 | 11| 12| 13| 14| 15 | 16| 17| 18| 19| 20 | 1| 2| 3| 4| 5 | 6| 7| 8| 9| 10|||2100||2100||2100||2100||2100|||2100|||2100||31101||31102||31103| | 4| 5| 6| 7| 8 | 9| 10| 11| 12| 13 .gen.receevent=1ifweek==t&arrest==1(114realchangesmade).listidweekarrestagewexpmart | weekarrest wexpmartevent | 0| 0| 0| 0| 0 | 0| 0| 0| 0| 0 | 0| 0| 0| 0| 0 | 0| 0| 0| 0| 1 | 0| 0| 0| 0| 0 | 0| 0| 0| 0| 0 | 0| 0| 0| 0| 0 | 0| 1| 0| 0| 0 | 0| 0| 0| 0| 0 | 0| 0| 0| 0| 0 | 0| 0| 0| 0| 0 | 0| 0| 0| 0| 0 | 0| 1| 0| 0| 0 | 0| 0| 0| 0| 0 | 0| 0| 0| 0| 0 | 0| 0| 0| 0| 0 | 0| 0| 0| 0| 0 | 0| 0| 0| 0| 0 | 0| 0| 0| 0| 0 | 0| 0| 0| 0| 0 | 0| 0| 0| 0| 0 | 0| 0| 0| 0| 0 | 0| 0| 0| 0.xi:logiteventagewexpmart,Logit Numberofobs=LRchi2(4)=Prob>chi2=Loglikelihood=- PseudoR2=event| Std. P>|z|[95%Conf.+age- --wexp- -.mar- --.t| _cons|-4.074031 -8.060.000-5.064355-.xi:logiteventagewexpmart,nologLogit Numberofobs=LRchi2(4)=Prob>chi2=Loglikelihood=- PseudoR2=event|OddsRatioStd. P>|z|[95%Conf.+age -..wexp -.mar -.t1.017161 logit模型所得到的回归系数等与我们通过其他模型(Cox模型)得到的系第七节PiecewiseExponential尽管Cox模型用处很广，目前已成为生存分析的主流方法，但对加速模型exponential模型的改进方法，即分段指数回归模型(Piecewiseexponentialmodels)。(Walder,LiandTreiman,1999,ASRarticle)。Cox模型，piecewiseexponential模型的一个优点是：我们可以利用它来风险率是如何随着时间例如我们想研究结7年之内7年之后的离数据：RecidivismData.***source:P.H.RossiandR.A.BerkandK.J.Lenihan(1980).Money,WorkandCrime:SomeExperimentalResults,AcademicPress:NewYork.*toestimatethepiecewiseexponentialmodelweneedtocreatea*“stacked”dataset.Thisiscommonlyused,powerful,andtrickinseveral*ofysis.STATAprovidesuswithaslickwaytocreateastackeddata*particularlysuitedtoestimationofapiecewiseexponential⑥该数据未能提供，类似的案例可以参见.use tobeginwith,weneedtocreateanIDvariable.Thereasonforthisisthatweareactuallygoingtocreateduplicates..gen.stsetweek,id(id)failure(arrest)id:idfailureevent:arrest!=0&arrest<.obs.timeinterval:(week[_n-1],week]exitonorbefore:432total 432obs.remaining,432114failuresinsinglefailure-per-subjectdata19809totalysistimeatrisk,atriskfromt=0earliestobservedentryt=0lastobservedexitt=nowwearereadytoissuetheamazing–stsplit-command.Whenwedoso,we*toindicatehowwearegoingtocreateastackingofthedata.Inthe*data,wewilldivide52weeksinto4equallengthsegmentsof13weekseach.*variableindicatingthedivisionisnamedas.stsplitJ,at(13,26,(1141observations(episodes)*letlookatwhatthedatahavebeen.listid_d_t_t0 |id_d_t_t0J||1.10002.11133.20004.21135.|301300 |3|31||4000||40|