p-Value from Z-ratio

Question

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.

Answer 1

In general, since emmeans is a widely used package you should probably assume that whatever computation it's doing is being done correctly: the danger is that it might not be doing what you want, or you might not understand what it's doing (or both). The most obvious thing that occurs to me is that computing pairwise contrasts automatically invokes a multiple-comparisons correction. From the vignette:

In its out-of-the-box configuration, pairs() sets two defaults for summary(): adjust = "tukey" (multiplicity adjustment), and infer = c(FALSE, TRUE) (test statistics, not confidence intervals). You may override these, of course, by calling summary() on the result with different values for these.

So if the contrast you're showing us is one of several pairwise contrasts, the p-value will be larger than that computed from a single Z-statistic (google provides lots of good hits for "tukey pairwise contrasts correction" , I wasn't sure which one to pick).