Stata:中断时间序列分析-itsa

发布时间:2022-07-15 阅读 1826

Stata连享会   主页 || 视频 || 推文 || 知乎 || Bilibili 站

温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。

New! lianxh 命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc, ihelp, rdbalance, gitee, installpkg

课程详情 https://gitee.com/lianxh/Course

课程主页 https://gitee.com/lianxh/Course

⛳ Stata 系列推文:

PDF下载 - 推文合集

作者:李亭 (上海财经大学)
邮箱2930175934@qq.com

编者按:本文主要摘译自下文,特此致谢!
Source:Linden A. Conducting interrupted time-series analysis for single-and multiple-group comparisons[J]. The Stata Journal, 2015, 15(2): 480-500. -PDF-


目录


中断时间序列 (ITS) 是公共卫生、公共政策和卫生服务领域的一种常见研究方法。该方法通过控制结果变量干预前的上升或下降趋势,来检验干预措施的效果,包括中断前后的水平变化和趋势变化。此外,中断可以是有计划的,如推出新的卫生政策,或非计划的,如意外的环境暴露。

1. 理论分析

中断时间序列 (ITS) 分为单组的 ITS 设计 (无控制组) 和多组的 ITS 设计 (存在一个或多个控制组)。在进行统计分析时,中断时间序列 (ITS) 对干预实施前后的两个时间段进行线性回归拟合,以检验政策实施后水平下降或升高的幅度是否有统计学意义。

1.1 公式和方法

当研究中只有单个组 (无控制组) 时,标准 ITSA 回归模型采用以下形式:

其中,Tt 是研究开始后的时间,Xt 是干预的虚拟变量 (干预前为 0,否则为 1) 。在单组设计中,β0 表示结果变量的截距或起始水平,β1 表示在引入干预之前结果变量的斜率或轨迹,β2 表示在引入干预之后结果变量的水平变化,β3 表示干预前和干预后回归斜率之间的差异。因此,β2 和 β3 的显著程度是我们关注的重点,分别表示干预时的瞬时效果和随着时间推移所造成的长期影响。

当存在控制组时,标准 ITSA 回归模型将扩展为以下形式:

相比较单组分析而言,多组分析中包含了 β4-β7 四个新的系数。其中,Z 为描述组别的虚拟变量 (干预组为 0,否则为 1)。β4 表示干预前处理与对照初始水平的差异,β5 表示干预前处理与对照斜率的差异,β6 表示处理与对照在干预开始后瞬时变化量的差异,β7 表示处理与对照在干预开始后和干预前斜率的差异。

1.2 适用的情形

中断时间序列 (ITS) 要求干预前后结果变量的变化趋势为线性,此时可用前述的分段回归分析。但时间序列数据可能存在季节性和周期性,因此必须考虑数据存在自相关的可能。此外,ITS 应用的样本含量,一般认为使用 ITS 模型干预前后不少于 10 年,尤其当样本量 (观察点) 少时,更多地反映了结果变量 “短期改变” 而不是 “长期趋势”。

ITS 是准实验设计,并不适合作因果关系的推断。ITS 分析可确定结果变量在干预后是否有水平和趋势改变,但不能确定这种改变与干预措施有明确的因果关系,因伴随的其他解释变量可直接或间接地影响结果变量。

ITS 设计还可在单组 (处理组) ITS 的基础上,再增加一个平行对照组,试验组在干预前后比较的同时,还可与对照组比较,以控制同期非实验因素的干扰,得到干预措施与结果联系的更可靠结论。

2. 命令介绍

本文将介绍 itsa 命令,它基于两种普通最小二乘 (OLS) 回归方法执行中断时间序列分析。该命令不仅可以对单组 (无控制组) 和多组 (有控制组) 进行分析,而且还可以估计存在多个处理期的情形。

使用中断时间序列 (ITS) 进行统计分析时必须考虑数据存在自相关的可能性。ITSA 方法使用的两种通用方法是自回归移动平均模型 (ARMA) 和旨在调整自相关的 OLS 回归模型。itsa 命令依赖于普通最小二乘 (OLS) 回归而不是基于自回归移动平均模型 (ARMA) 的回归方法,因为前者通常更灵活,适用范围更广泛。

命令安装:

ssc install itsa, replace

命令语法:

itsa depvar [indepvars] [if] [in] [weight], trperiod(numlist) 
  [single treatid(#) contid(numlist) prais lag(#) figure[(twoway_options)] 
  posttrend replace prefix(string) model_options]

其中,

  • trperiod(*numlist*):指定干预开始的时间,也可以指定多个时间段,该选项为必选项;
  • single:指定 itsa 将用于单组分析。相反,省略该项表示是多组分析;
  • treatid(*#*):指定当数据集包含多个样本时,正在研究的处理样本标识符。当数据集只包含一个样本数据时,该选项必须舍弃;
  • prais:指定拟合 prais 模型,如果没有指定该选项,itsa 将默认使用 newey 进行回归;
  • lag(*#*):指定选择 newey 模型时,应考虑的最大滞后期数。如果选择滞后 0 期,结果和使用 regress , vce(robust) 一致。此外 prais 模型和 lag() 不能同时使用,因为 prais 默认执行的是 AR(1) 过程;
  • figure:生成中断时间序列分析图命令,在多组分析中,则将控制组的观测值取平均后绘制;
  • posttrend:使用 lincom 模型产生处理后趋势估计值。在单组 ITSA 的情况下,会产生一个估计数。在多组 ITSA 的情况下,则对处理组、控制组和两组的差异进行估计。在有多个处理期的情况下,将为每个处理期制作一个单独的表格。

3. Stata 操作

3.1 单组 ITS 设计

本文使用的案例数据来自于 Abadie (2010) 对加州烟草控制计划效果的研究。1988 年 11 月美国加州通过了当代美国最大规模的控烟法-99 号提案 (anti-tobacco legislation),并于 1989 年 1 月开始生效。该法将加州的香烟消费税 (cigarette excise tax) 提高了每包 25 美分。我们使用 itsa 来评估 99 号提案在减少加州人均香烟销售方面的影响。更具体地说,我们评估了 99 号提案的引入是否导致人均香烟销售水平和趋势发生了变化。

首先,我们加载数据,并将数据集声明为面板。

. lxhuse cigsales.dta, clear
. tsset state year

Panel variable: state (strongly balanced)
 Time variable: year, 1970 to 2000
         Delta: 1 unit

接下来,我们指定一个单组 ITSA,加利福尼亚 (研究中第3个州) 作为处理组,1989 作为干预的开始,采用滞后一期的 newey 模型进行估计:

. itsa cigsale, single treat(3) trperiod(1989) lag(1) posttrend figure

Regression with Newey–West standard errors      Number of obs     =         31
Maximum lag = 1                                 F(  3,        27) =     331.45
                                                Prob > F          =     0.0000
------------------------------------------------------------------------------
             |             Newey–West
    _cigsale | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
          _t |     -1.779      0.383    -4.64   0.000       -2.566      -0.993
      _x1989 |    -20.058      4.724    -4.25   0.000      -29.752     -10.364
    _x_t1989 |     -1.495      0.437    -3.42   0.002       -2.391      -0.598
       _cons |    132.226      4.253    31.09   0.000      123.499     140.952
------------------------------------------------------------------------------

Postintervention Linear Trend: 1989
Treated: _b[_t]+_b[_x_t1989]
------------------------------------------------------------------------------
Linear Trend |      Coef.   Std. Err.      t              [95% Conf. Interval]
-------------+----------------------------------------------------------------
     Treated |  -3.274126   .2688039   -12.18   0.000    -3.825666   -2.722586
------------------------------------------------------------------------------

可以看出,人均香烟销售的起始水平估计为 132 包,在 1989 年之前,销量似乎每年都显著下降 1.78 包。在干预的第一年 (1989 年),人均香烟销量似乎显著下降了 20.06 包。随后,年销售趋势 (相对于干预前的趋势) 显著下降约人均每年 1.49 包。此外,lincom 模型估计的结果表明在 1989 年干预后,加州的人均香烟销售量以每年约 3.28 包的速度下降。

为了确保拟合一个正确的自相关结构的模型,我们使用 Actaum 和 Schaffer (2013) 的方法来检验数据的自相关性。

. actest, lag(6)

Cumby-Huizinga test for autocorrelation
  H0: variable is MA process up to order q
  HA: serial correlation present at specified lags >q
-----------------------------------------------------------------------------
  H0: q=0 (serially uncorrelated)        |  H0: q=specified lag-1
  HA: s.c. present at range specified    |  HA: s.c. present at lag specified
-----------------------------------------+-----------------------------------
    lags   |      chi2      df     p-val | lag |      chi2      df     p-val
-----------+-----------------------------+-----+-----------------------------
   1 -  1  |     15.242      1    0.0001 |   1 |     15.242      1    0.0001
   1 -  2  |     15.255      2    0.0005 |   2 |      3.300      1    0.0693
   1 -  3  |     15.325      3    0.0016 |   3 |      1.192      1    0.2749
   1 -  4  |     15.896      4    0.0032 |   4 |      0.000      1    0.9880
   1 -  5  |     16.057      5    0.0067 |   5 |      1.113      1    0.2914
   1 -  6  |     16.078      6    0.0133 |   6 |      2.051      1    0.1521
-----------------------------------------------------------------------------
  Test allows predetermined regressors/instruments
  Test requires conditional homoskedasticity

如表的右侧结果所示,自相关只出现在滞后 1 期时。因此,我们指定的滞后 1 期的初始模型应该正确地解释了这种自相关。

3.2 多组 ITS 设计

在这个例子中,我们使用多组设计来评估 99 号提案在减少加州人均香烟销售 (包装包装) 方面的影响。更具体地说,我们现在将加州的经验与数据文件中其他 38 个州的经验进行了比较。

. itsa cigsale, treat(3) trperiod(1989) lag(1) replace figure

Regression with Newey–West standard errors      Number of obs     =      1,209
Maximum lag = 1                                 F(  7,      1201) =     364.04
                                                Prob > F          =     0.0000
------------------------------------------------------------------------------
             |             Newey–West
    _cigsale | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
          _t |     -0.548      0.294    -1.86   0.063       -1.125       0.029
          _z |     -3.274      5.338    -0.61   0.540      -13.746       7.198
        _z_t |     -1.232      0.464    -2.65   0.008       -2.142      -0.321
      _x1989 |    -17.252      3.815    -4.52   0.000      -24.737      -9.766
    _x_t1989 |     -0.504      0.525    -0.96   0.338       -1.534       0.527
    _z_x1989 |     -2.806      5.842    -0.48   0.631      -14.268       8.655
  _z_x_t1989 |     -0.991      0.666    -1.49   0.137       -2.297       0.315
       _cons |    135.499      3.554    38.13   0.000      128.527     142.472
------------------------------------------------------------------------------

可以看出,加州与其他州之间的初始平均水平 (z) 差异并不显著,但干预前的平均趋势斜率 (_z_t) 差异具有显著性。下图也证实了这一点:在政策干预前,38 个州的平均香烟销售轨迹似乎一直高于加州。鉴于干预前趋势变化的这种差异,我们可以说其他 38 个州与加州没有可比性,从而对 (z_x1989) 和 (z_x_t1989) 的效果估计可能有偏差。因此,该模型可以通过将控制组的选择限制在这两个变量上值相似的对照组来改进。

接下来,我们将控制组限制在那些在起始水平和干预前变化趋势上与加州可比的州。作者将这种可比性定义为在初始水平 (z) 和干预前平均趋势斜率 (z_t)上的 p 值都大于 10%。作者经过迭代分析后,发现有三个州比较符合这个标准:科罗拉多州、爱达荷州和蒙大拿州。

. itsa cigsale, treat(3) trperiod(1989) contid(4 8 19) lag(1) replace posttrend figure

Regression with Newey–West standard errors      Number of obs     =        124
Maximum lag = 1                                 F(  7,       116) =     251.48
                                                Prob > F          =     0.0000
------------------------------------------------------------------------------
             |             Newey–West
    _cigsale | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
          _t |     -1.465      0.384    -3.82   0.000       -2.225      -0.704
          _z |      1.731      5.737     0.30   0.763       -9.632      13.095
        _z_t |     -0.315      0.533    -0.59   0.556       -1.371       0.741
      _x1989 |    -13.589      4.180    -3.25   0.002      -21.869      -5.309
    _x_t1989 |      0.475      0.499     0.95   0.344       -0.514       1.463
    _z_x1989 |     -6.469      6.185    -1.05   0.298      -18.720       5.781
  _z_x_t1989 |     -1.969      0.653    -3.01   0.003       -3.263      -0.675
       _cons |    130.495      4.009    32.55   0.000      122.553     138.436
------------------------------------------------------------------------------

Comparison of Linear Postintervention Trends: 1989
Treated    : _b[_t] + _b[_z_t] + _b[_x_t1989] + _b[_z_x_t1989]
Controls   : _b[_t] + _b[_x_t1989]
Difference : _b[_z_t] +  _b[_z_x_t1989]
------------------------------------------------------------------------------
Linear Trend |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     Treated |  -3.274126   .2593693   -12.62   0.000     -3.78784   -2.760413
    Controls |  -.9898601   .2882862    -3.43   0.001    -1.560847    -.418873
-------------+----------------------------------------------------------------
  Difference |  -2.284266   .3877903    -5.89   0.000    -3.052334   -1.516198
------------------------------------------------------------------------------

可以看出,处理组在结果变量起始水平 (z) 和干预前平均趋势斜率 (z_t) 上均与控制组具有可比性。虽然在干预的第一年 (z_x_1989) 没有统计学上显著的处理效果,这可能是政策具有一定的滞后性,但与控制组相比,干预后趋势有显著的下降,为每年 1.97 支人均香烟销量。此外,从 posttrend 分析结果也可以看出,处理组在干预后香烟销量下降了 3.27 包,控制组同期仅下降了 1 包,两者的差异为人均每年 2.28 包。

上述结果强调了在进行多组 ITSA 时,确保处理组和控制组在起始水平和干预前变化趋势上具有可比性的重要性。在选择合适的控制组时,可以使用一个迭代过程,将每个控制组单位分别与处理组进行比较。当 β4 和 β5 的 p 值大于指定阈值的组可以保留作为对照,纳入最终模型,这种方法也可以很容易地扩展到其他协变量。

4. 相关结论

虽然随机对照试验仍然是金标准的研究设计,但在某些情况下,例如当大规模的干预或政策变化针对整个人群时,这种设计不可行或不实用。当干预前和干预后的多个时间点的数据可用时,中断时间序列设计为评估政策处理效果提供了一个稳健的准自然实验替代方案。

在本文中,我们演示了 itsa 命令在 Stata 中的基本操作,并对单组和多组 ITS 设计进行了案例分析。实际上,itsa 命令可以分析存在多个处理期的情况。此外,itsa 命令还解决了其他重要问题,如选择和指定模型的标准、数据自相关检验、干预开始之前的稳健性测试,以及选择可比对照组等。

5. 参考文献

  • Abadie A, Diamond A, Hainmueller J. Synthetic control methods for comparative case studies: Estimating the effect of California’s tobacco control program[J]. Journal of the American statistical Association, 2010, 105(490): 493-505. -PDF-
  • 邵华, 邵华, 王琦琦, 等. 中断时间序列分析及其在公共卫生领域中的应用[J]. 中华流行病学杂志, 2015, 36(009): 1015-1017. -PDF-

6. 相关推文

Note:产生如下推文列表的 Stata 命令为:
lianxh 因果推断, m
安装最新版 lianxh 命令:
ssc install lianxh, replace

相关课程

免费公开课

最新课程-直播课

专题 嘉宾 直播/回看视频
最新专题 文本分析、机器学习、效率专题、生存分析等
研究设计 连玉君 我的特斯拉-实证研究设计-幻灯片-
面板模型 连玉君 动态面板模型-幻灯片-
面板模型 连玉君 直击面板数据模型 [免费公开课,2小时]
  • Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。

课程主页

课程主页

关于我们

  • Stata连享会 由中山大学连玉君老师团队创办,定期分享实证分析经验。
  • 连享会-主页知乎专栏,700+ 推文,实证分析不再抓狂。直播间 有很多视频课程,可以随时观看。
  • 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标,输入简要关键词,以便快速呈现历史推文,获取工具软件和数据下载。常见关键词:课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法

连享会小程序:扫一扫,看推文,看视频……

扫码加入连享会微信群,提问交流更方便

✏ 连享会-常见问题解答:
https://gitee.com/lianxh/Course/wikis

New! lianxhsongbl 命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh