Interrupted time series data
kaz_yos
2015-08-15
```
References
Experimental and Quasi-Experimental Designs for Generalized Causal Inference
Introductory Time Series with R
The use of segmented regression in analysing interrupted time series studies: an example in pre-hospital ambulance care: http://www.implementationscience.com/content/9/1/77
Wats: Wrap Around Time Series graphics: https://cran.r-project.org/web/packages/Wats/
Load packages
library(Wats)
library(ggplot2)
library(nlme)Load data
## Monthly Growth Fertility Rates (GFR) for 12 urban Oklahoma counties
data(CountyMonthBirthRate2014Version)
## Extract canadian county as an example
canadian <- subset(CountyMonthBirthRate2014Version, CountyName == "canadian")
## Date of Oklahoma bombing
dateBombing <- as.Date("1995-04-19")
## Indicator for post-bombing dates
canadian$postBombing <- (canadian$Date >= dateBombing)
## Date as numeric
canadian$DateNum <- as.numeric(canadian$Date)
## Center the date at bombing
canadian$DateNumCtr <- canadian$DateNum - as.numeric(dateBombing)Plot
ggplot(data = canadian, mapping = aes(x = Date, y = BirthRateMonthly)) +
layer(geom = "line") +
theme_bw() + theme(legend.key = element_blank())
Linear model fit for a correlogram
The lag 1 autocorrelation is around 0.25.
lmFit1 <- lm(BirthRateMonthly ~ DateNum, data = canadian)
## Auto- and Cross- Covariance and -Correlation Function Estimation
acf(resid(lmFit1))
Crude interrupted time series
\(E[Y] = \beta_0 + \beta_1 Time + \beta_2 Period2 + \beta_3 Time \times Period2\)
where
- \(\beta_0\) is the pre period intercept
- \(\beta_1\) is the pre period slope (baseline time trend)
- \(\beta_2\) is the immediate effect of the event
- \(\beta_3\) is the slope change after the event
- \(\beta_0 + \beta_2\) is the post period intercept
- \(\beta_1 + \beta_3\) is the post period slope
- \(Time\) is the date centered at the event (necessary for interpretability of immediate effect)
- \(Period2\) is the indicator for the time period after the event
## Fit model
glsFit1 <- gls(model = BirthRateMonthly ~ DateNumCtr + postBombing + DateNumCtr:postBombing,
data = canadian,
correlation = corAR1(0.25))
summary(glsFit1)## Generalized least squares fit by REML
## Model: BirthRateMonthly ~ DateNumCtr + postBombing + DateNumCtr:postBombing
## Data: canadian
## AIC BIC logLik
## 234.8338 251.3553 -111.4169
##
## Correlation Structure: AR(1)
## Formula: ~1
## Parameter estimate(s):
## Phi
## 0.260391
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 4.399954 0.17586388 25.019089 0.0000
## DateNumCtr -0.000278 0.00015751 -1.765280 0.0801
## postBombingTRUE 0.059139 0.25593277 0.231074 0.8177
## DateNumCtr:postBombingTRUE 0.000439 0.00025098 1.747655 0.0832
##
## Correlation:
## (Intr) DtNmCt pBTRUE
## DateNumCtr 0.860
## postBombingTRUE -0.665 -0.572
## DateNumCtr:postBombingTRUE -0.559 -0.644 -0.122
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.27069072 -0.70619943 0.03102496 0.60955441 2.29540604
##
## Residual standard error: 0.5530618
## Degrees of freedom: 120 total; 116 residualOverplot prediction
## Create prediction dataset
newdata <- data.frame(DateNumCtr = seq(min(canadian$DateNumCtr), max(canadian$DateNumCtr), by = 1))
newdata$postBombing <- (newdata$DateNumCtr >= 0)
## Predict
newdata$BirthRateMonthly <- predict(glsFit1, newdata = newdata)
ggplot(data = canadian, mapping = aes(x = DateNumCtr, y = BirthRateMonthly)) +
layer(geom = "line") +
layer(data = subset(newdata, DateNumCtr < 0), geom = "line", color = "red", size = 1.5) +
layer(data = subset(newdata, DateNumCtr >= 0), geom = "line", color = "red", size = 1.5) +
theme_bw() + theme(legend.key = element_blank())
Explicit modeling of monthly fluctuations
The immediate change postBombingTRUE is non-significant. The trend change DateNumCtr:postBombingTRUE is positive and significant.
## Fit
glsFit2 <- gls(model = BirthRateMonthly ~ DateNumCtr + postBombing + DateNumCtr:postBombing + factor(Month),
data = canadian,
correlation = corAR1(0.25))
summary(glsFit2)## Generalized least squares fit by REML
## Model: BirthRateMonthly ~ DateNumCtr + postBombing + DateNumCtr:postBombing + factor(Month)
## Data: canadian
## AIC BIC logLik
## 227.1337 272.251 -96.56685
##
## Correlation Structure: AR(1)
## Formula: ~1
## Parameter estimate(s):
## Phi
## 0.2130256
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 4.056658 0.19580489 20.717860 0.0000
## DateNumCtr -0.000271 0.00012568 -2.152381 0.0337
## postBombingTRUE -0.010639 0.20603440 -0.051637 0.9589
## factor(Month)2 -0.215932 0.18251667 -1.183083 0.2394
## factor(Month)3 0.289306 0.20102227 1.439174 0.1531
## factor(Month)4 0.141948 0.20482160 0.693033 0.4898
## factor(Month)5 0.559483 0.20585345 2.717869 0.0077
## factor(Month)6 0.524528 0.20591945 2.547250 0.0123
## factor(Month)7 0.763552 0.20589578 3.708439 0.0003
## factor(Month)8 0.637137 0.20586202 3.094971 0.0025
## factor(Month)9 0.793359 0.20574894 3.855955 0.0002
## factor(Month)10 0.415827 0.20511644 2.027271 0.0452
## factor(Month)11 -0.062976 0.20186108 -0.311976 0.7557
## factor(Month)12 0.645521 0.18516481 3.486198 0.0007
## DateNumCtr:postBombingTRUE 0.000445 0.00019924 2.231004 0.0278
##
## Correlation:
## (Intr) DtNmCt pBTRUE fc(M)2 fc(M)3 fc(M)4 fc(M)5 fc(M)6 fc(M)7 fc(M)8 fc(M)9 f(M)10
## DateNumCtr 0.634
## postBombingTRUE -0.476 -0.578
## factor(Month)2 -0.474 -0.013 0.010
## factor(Month)3 -0.527 -0.023 0.020 0.550
## factor(Month)4 -0.544 -0.033 0.031 0.465 0.595
## factor(Month)5 -0.503 0.015 -0.056 0.446 0.508 0.601
## factor(Month)6 -0.510 0.006 -0.044 0.443 0.491 0.517 0.606
## factor(Month)7 -0.516 -0.003 -0.031 0.442 0.488 0.500 0.522 0.606
## factor(Month)8 -0.522 -0.013 -0.018 0.442 0.487 0.497 0.504 0.522 0.606
## factor(Month)9 -0.527 -0.022 -0.006 0.442 0.487 0.496 0.499 0.504 0.522 0.606
## factor(Month)10 -0.530 -0.031 0.007 0.440 0.485 0.494 0.496 0.498 0.502 0.520 0.605
## factor(Month)11 -0.527 -0.039 0.019 0.433 0.477 0.486 0.486 0.488 0.490 0.494 0.513 0.597
## factor(Month)12 -0.486 -0.041 0.030 0.397 0.438 0.447 0.446 0.447 0.449 0.450 0.455 0.473
## DateNumCtr:postBombingTRUE -0.404 -0.637 -0.125 0.009 0.010 0.010 0.020 0.017 0.014 0.011 0.008 0.005
## f(M)11 f(M)12
## DateNumCtr
## postBombingTRUE
## factor(Month)2
## factor(Month)3
## factor(Month)4
## factor(Month)5
## factor(Month)6
## factor(Month)7
## factor(Month)8
## factor(Month)9
## factor(Month)10
## factor(Month)11
## factor(Month)12 0.558
## DateNumCtr:postBombingTRUE 0.001 -0.009
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.16271810 -0.60854485 -0.01295044 0.61383062 2.38590679
##
## Residual standard error: 0.4603772
## Degrees of freedom: 120 total; 105 residual## Provide dataset for prediction
newdata2 <- data.frame(Date = seq(min(canadian$Date), max(canadian$Date), by = 1))
newdata2$DateNum <- as.numeric(newdata2$Date)
newdata2$postBombing <- (newdata2$Date >= dateBombing)
newdata2$DateNum <- as.numeric(newdata2$Date)
newdata2$DateNumCtr <- newdata2$DateNum - as.numeric(dateBombing)
newdata2$Month <- as.numeric(format(newdata2$Date, "%m"))
## Predict
newdata2$BirthRateMonthly <- predict(glsFit2, newdata = newdata2)
## Plot
ggplot(data = canadian, mapping = aes(x = DateNumCtr, y = BirthRateMonthly)) +
layer(geom = "line") +
layer(data = subset(newdata2, DateNumCtr < 0), geom = "line", color = "red", size = 1) +
layer(data = subset(newdata2, DateNumCtr >= 0), geom = "line", color = "red", size = 1) +
theme_bw() + theme(legend.key = element_blank())
Modeling of seasons using trigonometric functions
Reference: Introductory Time Series with R
Rescale the time variable to an interpretable scale (year here). sin1 and cos1 is represent a yearly fluctuation. sin2 and cos2 represent a twice-a-year fluctuation, and so on. A 0-to-1 change in the time variable create a cycle. To create a cycle twice more frequent, multiply the time variable by 2.
In this case the result was consistent with the month indicator method.
## Sine waves
## Yearly cycle
canadian$sin1 <- sin(2 * pi * 1 * (canadian$DateNumCtr / 365.25))
## 6-month cycle
canadian$sin2 <- sin(2 * pi * 2 * (canadian$DateNumCtr / 365.25))
## 4-month cycle
canadian$sin3 <- sin(2 * pi * 3 * (canadian$DateNumCtr / 365.25))
## 3-month cycle
canadian$sin4 <- sin(2 * pi * 4 * (canadian$DateNumCtr / 365.25))
## Cosine waves
## Yearly cycle
canadian$cos1 <- cos(2 * pi * 1 * (canadian$DateNumCtr / 365.25))
## 6-month cycle
canadian$cos2 <- cos(2 * pi * 2 * (canadian$DateNumCtr / 365.25))
## 4-month cycle
canadian$cos3 <- cos(2 * pi * 3 * (canadian$DateNumCtr / 365.25))
## 3-month cycle
canadian$cos4 <- cos(2 * pi * 4 * (canadian$DateNumCtr / 365.25))
## Fit
glsFit3 <- gls(model = BirthRateMonthly ~ DateNumCtr + postBombing + DateNumCtr:postBombing +
sin1 + sin2 + sin3 + sin4 + cos1 + cos2 + cos3 + cos4,
data = canadian,
correlation = corAR1(0.25))
summary(glsFit3)## Generalized least squares fit by REML
## Model: BirthRateMonthly ~ DateNumCtr + postBombing + DateNumCtr:postBombing + sin1 + sin2 + sin3 + sin4 + cos1 + cos2 + cos3 + cos4
## Data: canadian
## AIC BIC logLik
## 244.2226 281.7724 -108.1113
##
## Correlation Structure: AR(1)
## Formula: ~1
## Parameter estimate(s):
## Phi
## 0.1443314
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 4.420599 0.13561687 32.59624 0.0000
## DateNumCtr -0.000278 0.00012155 -2.28856 0.0241
## postBombingTRUE 0.016509 0.19997853 0.08255 0.9344
## sin1 0.330713 0.06937414 4.76710 0.0000
## sin2 0.025177 0.06471866 0.38903 0.6980
## sin3 0.042194 0.06041427 0.69842 0.4864
## sin4 -0.181620 0.05667975 -3.20433 0.0018
## cos1 -0.088323 0.06935125 -1.27356 0.2056
## cos2 -0.031353 0.06516030 -0.48116 0.6314
## cos3 0.126563 0.06006912 2.10696 0.0374
## cos4 -0.084163 0.05609992 -1.50023 0.1365
## DateNumCtr:postBombingTRUE 0.000441 0.00019223 2.29222 0.0238
##
## Correlation:
## (Intr) DtNmCt pBTRUE sin1 sin2 sin3 sin4 cos1 cos2 cos3 cos4
## DateNumCtr 0.862
## postBombingTRUE -0.673 -0.581
## sin1 0.049 0.029 -0.078
## sin2 0.023 0.008 -0.030 -0.004
## sin3 0.035 0.037 -0.047 -0.004 0.000
## sin4 0.020 0.018 -0.030 0.008 0.004 0.012
## cos1 0.010 0.035 -0.026 0.002 0.006 0.008 -0.011
## cos2 -0.026 -0.032 0.033 -0.003 0.004 -0.017 -0.011 -0.001
## cos3 -0.017 -0.022 0.031 -0.001 -0.019 -0.015 -0.023 0.004 -0.002
## cos4 -0.004 0.002 0.010 -0.019 -0.015 -0.022 0.012 0.003 -0.006 -0.015
## DateNumCtr:postBombingTRUE -0.549 -0.635 -0.126 0.019 -0.002 0.005 0.006 0.016 0.001 -0.007 0.000
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.17610167 -0.58007853 -0.05586155 0.68640698 2.16931067
##
## Residual standard error: 0.4755701
## Degrees of freedom: 120 total; 108 residual## Provide dataset for prediction
newdata3 <- data.frame(Date = seq(min(canadian$Date), max(canadian$Date), by = 1))
newdata3$DateNum <- as.numeric(newdata3$Date)
newdata3$postBombing <- (newdata3$Date >= dateBombing)
newdata3$DateNum <- as.numeric(newdata3$Date)
newdata3$DateNumCtr <- newdata3$DateNum - as.numeric(dateBombing)
newdata3$sin1 <- sin(2 * pi * 1 * (newdata3$DateNumCtr / 365.25))
newdata3$sin2 <- sin(2 * pi * 2 * (newdata3$DateNumCtr / 365.25))
newdata3$sin3 <- sin(2 * pi * 3 * (newdata3$DateNumCtr / 365.25))
newdata3$sin4 <- sin(2 * pi * 4 * (newdata3$DateNumCtr / 365.25))
newdata3$cos1 <- cos(2 * pi * 1 * (newdata3$DateNumCtr / 365.25))
newdata3$cos2 <- cos(2 * pi * 2 * (newdata3$DateNumCtr / 365.25))
newdata3$cos3 <- cos(2 * pi * 3 * (newdata3$DateNumCtr / 365.25))
newdata3$cos4 <- cos(2 * pi * 4 * (newdata3$DateNumCtr / 365.25))
## Predict
newdata3$BirthRateMonthly <- predict(glsFit3, newdata = newdata3)
## Plot
ggplot(data = canadian, mapping = aes(x = DateNumCtr, y = BirthRateMonthly)) +
layer(geom = "line") +
layer(data = subset(newdata3, DateNumCtr < 0), geom = "line", color = "red", size = 1) +
layer(data = subset(newdata3, DateNumCtr >= 0), geom = "line", color = "red", size = 1) +
theme_bw() + theme(legend.key = element_blank())
- Top Page: http://rpubs.com/kaz_yos/
- github: https://github.com/kaz-yos