I have run an Interrupted Time Series Analysis based upon the below code:
glm(`Subject Total` ~ Quarter + int2 + time_since_intervention2 ,
df, family = "poisson")
I have used the emmeans package to estimate the pairwise difference between the counterfactual and point estimate and get the below output:
contrast estimate SE df z.ratio p.value
Quarter20 int21 time_since_intervention24 - Quarter20 int20 time_since_intervention20 -0.341 0.160 Inf -2.140 0.1406
The above estimate (0.341) cross-checks against manually derived outcomes. However I had a question about the p-value included within the above output and the manually derived equivalent [undertaken as a means of checking]. The p-value in the above is 0.146 [non-significant]. However, when calculated directly from the z-ratio (2*pnorm(-2.140)) I get a p-value of 0.03. Is the emmeans output correct (am I doing something wrong by assuming pnorm?) or is the manually calculated more likely to be accurate?
UPDATE:
The data frame is as below. Quarters represent time. Subject Total the outcome. Int2 a dummy variable to identify the point of intervention (0 pre/1 post). Time_since_intervention2 another dummy variable 0 prior to intervention 1:8 after.
> df[,c(1,2,9,11)]
Quarter Subject Total int2 time_since_intervention2
1 1 33 0 0
2 2 32 0 0
3 3 35 0 0
4 4 34 0 0
5 5 23 0 0
6 6 34 0 0
7 7 33 0 0
8 8 24 0 0
9 9 31 0 0
10 10 32 0 0
11 11 21 0 0
12 12 26 0 0
13 13 22 0 0
14 14 28 0 0
15 15 27 0 0
16 16 22 0 0
17 17 14 1 1
18 18 16 1 2
19 19 20 1 3
20 20 19 1 4
21 21 13 1 5
22 22 15 1 6
23 23 16 1 7
24 24 8 1 8
Call:
glm(formula = `Subject Total` ~ Quarter + int2 + time_since_intervention2,
family = "poisson", data = df)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.4769 -0.5111 0.1240 0.6103 0.9128
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.54584 0.09396 37.737 <0.0000000000000002 ***
Quarter -0.02348 0.01018 -2.306 0.0211 *
int2 -0.23652 0.21356 -1.108 0.2681
time_since_intervention2 -0.02624 0.04112 -0.638 0.5234
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 63.602 on 23 degrees of freedom
Residual deviance: 13.368 on 20 degrees of freedom
AIC: 140.54
Number of Fisher Scoring iterations: 4
The above suggests that the level change at the intervention point was non-significant.
We want to report the difference between the point estimate at Quarter 20 with the corresponding counterfactual (extrapolation of pre-policy trends/patterns) at the same point. At the moment I have done that using emmeans pairwise comparison.
emmeans(
object = fit1a,
specs = c("Quarter", "int2", "time_since_intervention2"),
at = list(
Quarter = c(20),
int2 = c(0, 1),
time_since_intervention2 = c(0,4)
)
) |>
contrast(method = "revpairwise")
contrast(method = "revpairwise")
contrast estimate SE df
Quarter20 int21 time_since_intervention20 - Quarter20 int20 time_since_intervention20 -0.237 0.214 Inf
Quarter20 int20 time_since_intervention24 - Quarter20 int20 time_since_intervention20 -0.105 0.164 Inf
Quarter20 int20 time_since_intervention24 - Quarter20 int21 time_since_intervention20 0.132 0.346 Inf
Quarter20 int21 time_since_intervention24 - Quarter20 int20 time_since_intervention20 -0.341 0.160 Inf
Quarter20 int21 time_since_intervention24 - Quarter20 int21 time_since_intervention20 -0.105 0.164 Inf
Quarter20 int21 time_since_intervention24 - Quarter20 int20 time_since_intervention24 -0.237 0.214 Inf
z.ratio p.value
-1.108 0.6849
-0.638 0.9197
0.380 0.9813
-2.140 0.1406
-0.638 0.9197
-1.108 0.6849
The only outcome that can exist is Quarter20 int21 time_since_intervention24 - Quarter20 int20 time_since_intervention20.
Not sure I'm doing this the right way.