通常，时间序列分析中对前后影响的评估称为中断时间序列。这是一种非常通用的建模方法，用于检验强假设：

$\mathcal{H}_0: \mu_{ijt} = f_i(t)$与$\mathcal{H}_1 : \mu_{ijt} = f_i(t) + \beta(t)X_{ijt}$

其中是个体在时间的治疗分配。最简单的例子是将视为一个常数函数，将视为 0,1 指标 0：干预前 1：干预前后。即使干预的实际“效果”与此不同，该测试也可以检测多种场景中的差异，例如，如果是任何非零函数，则工作常数参数将估计对干预的时间平均积极响应，并且是非零的。$X_{ijt}$$i$$t$$\beta$$X_{ijt}$$\beta(t)$$\beta$

事前干预的时间序列分析的一个挑战是使用参数建模方法进行自相关。通过多次重复时间和函数，可以将趋势分解为滞后效应、季节性效应等。这将消除误差项中自相关的需要。因此没有必要使用预测，但模型本身直接预测在干预后时间段内将观察到的情况。

datasets考虑R 包中著名的 Air Passengers 数据。

## construct an analytic dataset to predict time trend using auto-regressive and seasonal components
AirPassengers <- data.frame('flights'=as.numeric(AirPassengers))
AirPassengers$month <- factor(month.name, levels=month.name)
AirPassengers$year <- rep(1949:1960, each=12)
AirPassengers$lag <- c(NA, AirPassengers$flights[-nrow(AirPassengers)])

plot(AirPassengers$flights, type='l')

AirPassengers$fitted <- exp(predict(lm(log(flights) ~ month + year, data=AirPassengers)))
lines(AirPassengers$fitted, col='red')

很明显，这可以很好地预测基于时间的趋势。但是，如果您对关于“飞行增加”是否在 1955 年发布的假设检验感兴趣，您可以更新数据集以包含一个 0/​​1 指示符，以指示时间段是否在该点之后并进行测试它在线性模型中的意义。

例如：

library(lmtest)
library(sandwich)
AirPassengers$post <- AirPassengers$year >= 1955
fit <- lm(log(flights) ~ month + year + post, data=AirPassengers)
coeftest(fit, vcov. = vcovHC)['postTRUE', ]

给我：

> coeftest(fit, vcov. = vcovHC)['postTRUE', ]
  Estimate Std. Error    t value   Pr(>|t|) 
0.03720327 0.01783242 2.08627126 0.03890842 

这是虚假发现的一个很好的例子，以及实际上并不显着的统计显着效果。通过允许月份特定效应之间的异质性，可以进行更一般的测试。

nullmodel <- lm(log(flights) ~ month + year, data=AirPassengers)
fullmodel <- lm(log(flights) ~ post*month + year, data=AirPassengers)
waldtest(nullmodel, fullmodel, vcov=vcovHC, test='Chisq')

这两个都是用于分段回归的“中断时间序列”的一般方法的示例。这是一个定义松散的术语，我对作者在大多数情况下描述他们的确切方法时使用的细节很少感到有点失望。

为方便起见重复您的数据生成代码...

set.seed(101)  ## don't to forget to set the seed for reproducibility
##generate a comparable time-series ts() 
base <- rnbinom(n = 120, size = 1400, prob = 0.5)
season <- rep(c(600, 400, 150, 0, -50, -80, -300, -600, 50, 100, 200, 300), 10)

## dangerous to name your time-series the same as the ts() function
ts0 <- ts(base + season, start = c(2000,1), end = c(2009,12), frequency = 12)
##generate the independent variables
month.f <- factor(rep(1:12, 10))
dummy.months <- model.matrix(~month.f +0)
lag1 <- lag(ts0); lag2 <- lag(lag(ts0))
time.pre <- 72
time.post <- 48

我正在稍微改变你的虚拟变量。为了获得最大的可解释性，我认为您希望将截距和干预前斜率一直应用到整个数据集，以便“干预指标”和“干预后斜率”代表与干预前行为的差异。

pre.interv.trend <- 1:(time.pre+time.post)
interv.indicator <- c(rep(0, time.pre), rep(1, time.post))
post.interv.trend <- c(rep(0, time.pre), 1:time.post)
df <- data.frame(ts0, dummy.months,lag1,lag2,interv.indicator,
                 pre.interv.trend,post.interv.trend)

The fitted model:

    fit <- MASS::glm.nb(ts0 ~. + 0, data = df)

我们要分离的术语（唯一代表干预效果的术语）：

    params <- c("interv.indicator","post.interv.trend")

使用predict(.,type="terms") 几乎但并不完全符合您的要求。它将得到正确的预测，但置信区间未能考虑项之间的协方差......

    pp <- predict(fit,type="terms",terms=params,se.fit=TRUE)
    pp.fit <- rowSums(pp$fit)
    pp.se <- sqrt(rowSums(pp$se.fit^2))

这有点乏味，但实际上是正确的......

    m <- model.matrix(ts0~.+0,data=df)
    m2 <- m[,params]
    pred <- drop(m2 %*% coef(fit)[params])
    predsd <- sqrt(diag(m2 %*% vcov(fit)[params,params] %*% t(m2)))

绘制结果（细线忽略协方差，粗线包括它）：

    par(las=1,bty="l")
    matplot(cbind(pp.fit-1.96*pp.se,pp.fit,pp.fit+1.96*pp.se),
            type="l",lty=c(2,1,2),col=c(1,2,1),
            ylab="intervention effect (link scale)")
    matlines(cbind(pred-1.96*predsd,pred,pred+1.96*predsd),
        lty=c(2,1,2),type="l",col=c(1,2,1),lwd=2)