Suppose I have a panel dataset with monthly observations over 10 years. I have a simple dummy intervention, where some policy is put in place around the Spring in every year and only affects some people or groups. In other words, the units over time experience multiple shocks. Examples include something like emergency room re-entry or a crime policy affecting districts during warm weather. Note: the dummy is turning 'on' and 'off' over time and there is a seasonal component observed in the raw trends over time.
Let's assume I'm observing aggregate hospital admissions over time. To include hospital-specific linear time trends, I would multiply each hospital dummy by a continuous linear time trend variable, where $t = 1, 2, 3, ... , T$. To assess the robustness of the findings, it is common to estimate the following equation:
$$ y_{it} = \alpha_{0i} + \beta_i (t \times \alpha_{1i}) + \lambda_t + \delta D_{it} + u_{it}, $$
where $D_{it}$ is a treatment dummy equal to 1 if the facility was treated and is in a post-treatment time period, 0 otherwise. The parameters $\alpha_{0,i}$ and $\lambda_t$ denote hospital and month fixed effects, respectively. The interaction of the continuous time trend variable $t$ with each $\alpha_{1,i}$ gives each $i$ it's own unique linear time trend. This works well and is quite common in the literature with one intervention. But suppose the policy is removed and then reinstated in subsequent years. In other words, the intervention repeats every year around Spring time.
Here is what I am observing in practice. The toy dataset shows visits to a facility over 24 months. Every 12 months I see this inverted U-shape. The shaded regions represent the shocks observed over time.
# Loading the required packages
library(tidyr)
library(dplyr)
# Creating some fake data
set.seed(1987)
y <- c(8.4, 10.0, 11.8, 12.2, 13.1, 13.3, 12.0, 12.4, 12.0, 10.3, 10.0, 10.5, 9.3, 8.0, 8.1, 10.1 , 11.5, 12.1, 12.5, 12.1, 10.7, 8.8, 8.7, 7.0)
shocks <- c(5, 6, 18, 19) # on/off shocks
df <-
tibble(unit = rep(1:10, each = 24),
month = rep(1:24, 10),
y = rep(y, 10) + rnorm(240, 1, 5)
) %>%
mutate(group = ifelse(unit > 5, "Exposed", "Unexposed"),
intervention = ifelse(unit > 5 & month %in% shocks, 1, 0),
time = rep(1:10, 24), # linear time trend
sine = sin(2*pi*time/12), # sin()
cosine = cos(2*pi*time/12) # cos()
)
# Producing a fitted trend line for each group
# Shaded regions show transient shocks over time
# Some units become exposed to a policy and others do not
df %>%
ggplot(., aes(x = factor(month), y = y, color = factor(group))) +
geom_smooth(aes(group = group)) +
labs(x = "Month (24 Periods)", y = "Mean Outcome") +
annotate("rect", xmin = 5, xmax = 7, ymin = -Inf, ymax = Inf, fill = "maroon", alpha = .2) +
annotate("rect", xmin = 18, xmax = 20, ymin = -Inf, ymax = Inf, fill = "maroon", alpha = .2) +
scale_color_manual(name = "Group:", values = c("Exposed" = "maroon",
"Unexposed" = "royalblue")) +
theme_classic() +
theme(
legend.position = "bottom"
)
Questions:
(1) Is it valid to allow each unit to have their own sinusoidal trend, which is allowed to flexibly model the cyclical pattern over time between the two groups?
(2) Is there a better approach? I don't see how adding a quadratic time trend would help given that the intervention repeats. I would argue that including a higher order polynomial term would be too demanding. Or maybe not?
To be clear, this is more of a robustness check—not the main specification. I don't think it has been done before, in part because most quasi-experimental evaluations involve units with one treatment history. In applied work, the linear time trend is a way to adjust for the possibility that the treatment group and the control group were on somewhat different growth trajectories before the shock. But with multiple interventions, I wanted to adjust for any differential "cyclical" patterns observed over time, as they seem to repeat every year with the reintroduction of the policy.
Please review my R code below. I am somewhat new to including sinusoidal trends. I welcome any insight or criticism of this approach.
# I multiplied each unit by sin() and cos() separately
> summary(mod <- lm(y ~ as.factor(unit)*sine + as.factor(unit)*cosine + as.factor(month) + intervention, data = df))
Call:
lm(formula = y ~ as.factor(unit) * sine + as.factor(unit) * cosine +
as.factor(month) + intervention, data = df)
Residuals:
Min 1Q Median 3Q Max
-11.4567 -3.1386 -0.0271 3.0985 10.8461
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.8550 2.8408 3.821 0.000181 ***
as.factor(unit)2 -0.6000 1.9757 -0.304 0.761700
as.factor(unit)3 0.8304 1.7146 0.484 0.628729
as.factor(unit)4 -1.5308 1.7592 -0.870 0.385329
as.factor(unit)5 1.4748 1.9535 0.755 0.451217
as.factor(unit)6 1.9856 1.6350 1.214 0.226108
as.factor(unit)7 2.5895 2.0160 1.284 0.200584
as.factor(unit)8 1.2193 1.7436 0.699 0.485217
as.factor(unit)9 0.5477 1.8026 0.304 0.761602
as.factor(unit)10 3.3756 1.9880 1.698 0.091192 .
sine 2.6086 2.8070 0.929 0.353920
cosine -4.1625 3.3955 -1.226 0.221802
as.factor(month)2 0.2706 2.7212 0.099 0.920883
as.factor(month)3 3.8197 3.8195 1.000 0.318591
as.factor(month)4 -3.7247 4.6494 -0.801 0.424083
as.factor(month)5 -1.4893 5.5054 -0.271 0.787066
as.factor(month)6 -1.3394 5.7128 -0.234 0.814891
as.factor(month)7 2.3822 5.3802 0.443 0.658440
as.factor(month)8 2.3526 4.8477 0.485 0.628043
as.factor(month)9 3.2434 3.7984 0.854 0.394272
as.factor(month)10 3.3325 2.9592 1.126 0.261546
as.factor(month)11 2.8721 2.3530 1.221 0.223786
as.factor(month)12 0.3874 2.7212 0.142 0.886934
as.factor(month)13 -0.7625 3.8195 -0.200 0.841997
as.factor(month)14 -7.9720 4.6494 -1.715 0.088077 .
as.factor(month)15 -7.7604 5.4126 -1.434 0.153318
as.factor(month)16 -2.0334 5.6234 -0.362 0.718065
as.factor(month)17 -1.1144 5.3802 -0.207 0.836129
as.factor(month)18 2.9657 4.9512 0.599 0.549912
as.factor(month)19 1.2505 3.9296 0.318 0.750660
as.factor(month)20 4.6338 2.9592 1.566 0.119072
as.factor(month)21 1.9352 2.3530 0.822 0.411883
as.factor(month)22 -3.5087 2.7212 -1.289 0.198861
as.factor(month)23 -4.7884 3.8195 -1.254 0.211542
as.factor(month)24 -6.3957 4.6494 -1.376 0.170601
intervention -0.7227 2.0136 -0.359 0.720054
as.factor(unit)2:sine -2.5287 5.0387 -0.502 0.616358
as.factor(unit)3:sine -0.1849 3.4525 -0.054 0.957344
as.factor(unit)4:sine -0.7943 3.5224 -0.226 0.821831
as.factor(unit)5:sine -2.1947 4.9591 -0.443 0.658593
as.factor(unit)6:sine -1.7032 2.0679 -0.824 0.411209
as.factor(unit)7:sine -6.4227 5.0401 -1.274 0.204138
as.factor(unit)8:sine -0.3358 3.4440 -0.098 0.922423
as.factor(unit)9:sine -4.5655 3.5362 -1.291 0.198271
as.factor(unit)10:sine -2.6118 4.9475 -0.528 0.598196
as.factor(unit)2:cosine 1.4540 6.0715 0.239 0.810997
as.factor(unit)3:cosine 6.8385 4.2671 1.603 0.110720
as.factor(unit)4:cosine -1.1452 4.1782 -0.274 0.784312
as.factor(unit)5:cosine 7.7818 6.1802 1.259 0.209551
as.factor(unit)6:cosine 1.2134 2.4830 0.489 0.625636
as.factor(unit)7:cosine 5.3385 6.0575 0.881 0.379291
as.factor(unit)8:cosine 3.0046 4.2863 0.701 0.484185
as.factor(unit)9:cosine 0.2422 4.1670 0.058 0.953707
as.factor(unit)10:cosine 6.4106 6.1764 1.038 0.300653
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 5.262 on 186 degrees of freedom
Multiple R-squared: 0.2956, Adjusted R-squared: 0.0949
F-statistic: 1.473 on 53 and 186 DF, p-value: 0.03175
All that matters is the estimate on intervention
, at least for reporting purposes. I'm curious to see if effects hold even after the inclusion of such trends. I understand this is a very demanding procedure. To offer some perspective, I'm working with about 60 units and each unit is observed over 120 months.