How can the magnitude of an intervention be quantified in a segmented time series regression?
I am attempting to replicate the methodology of Decline in pneumonia admissions after routine childhood immunisation with pneumococcal conjugate vaccine in the USA: a time-series analysis. There are several other published papers with similar methodology where the methods chapter is not informative. All of these papers cite Wagner et al., however that paper describes a much simpler linear model.
I have a 120 month time series of count data that exhibits strong seasonal variation, has a baseline negative trend and a known public health intervention at $t = 72$ that would be expected to decrease the counts (increase the negative trend). I have fitted this with a negative binomial model with AR(2) errors and dummy variables for months, pre-intervention trend, an intervention indicator and post-intervention trend. I did this in R statistics with the function
glm.nb from the
#In R version 3.3.1 with packages dplyr, MASS #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) ts <- ts(base + season, start = c(2000,1), end = c(2009,12), frequency = 12) #generate the independant variables month.f <- factor(rep(1:12, 10)) dummy.months <- model.matrix(~month.f +0) require(dplyr); lag1 <- lag(ts); lag2 <- lag(lag(ts)) time.interv <- 72 pre.interv.trend <- c(1:time.interv, rep(0, 48)) interv.indicator <- c(rep(0, time.interv), rep(1, 48)) post.interv.trend <- c(rep(0, time.interv), 1:48) df <- cbind.data.frame(ts, dummy.months,lag1,lag2,interv.indicator,pre.interv.trend,post.interv.trend) #the fitted model require(MASS); fit <- glm.nb(ts ~. + 0, data = df)
I have attempted several approaches
- I have tried using the
forecastpackage to forecast the pre-intervention time series and then subtract the observed values from the predicted. However, the 95%CI intervals became so large that there would have been no theoretical way for the observed time series to fall outside them.
- I have refit the model without the intervention variables and subtracted
fit_nonintervention. The refitted model exhibits fairly similar fitted values with an overall decreased model fit.