rvlenth / emmeans

Estimated marginal means
https://rvlenth.github.io/emmeans/
341 stars 30 forks source link

Regridding outcome of interrupted time series #420

Closed jrahilly14 closed 1 year ago

jrahilly14 commented 1 year ago

Hi,

I've reached a point where I just can't quite figure out how to do what I need with the package and though my best option would be going to the source. I've read the vignettes, but can't equate to the model we have developed.

I have conducted an Interrupted Time Series Analysis, which entails fitting a segmented regression. The data therefore consists of an outcome (in this case number of new retail outlets), time (in this instance quarters), a dummy variable identifying whether the time period is prior to or after intervention (0 or 1) and an elapsed time variable (0 for all time points prior to intervention and 1:8 for all after).

I fit the model:

fit1a <- glm(outcome ~ time + intervention + elapsed, family = "poisson") 

and can derive a simple summary.

However, it is our intention to supply a difference between the point estimate of the model and the predicted counterfactual (continuation of prior trend) at time 20. I think have figured out how to do so with the emmeans package...however tis may be incorrect....

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")|>
  confint()

But I want to convert the confidence interval back to the original count. I assume I do so using the regrid function, but cannot get to grips with how I would do so in my particualr circumstance. I apologise for such a basic question, but am sort of losing hope with tbh

rvlenth commented 1 year ago

(I edited the formatting of your question for clarity.)

I'm going to save an intermediate step to make it simpler to talk about:

EMM <- 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)
    )
) 

So the question is what to do with EMM, knowing that estimates are on the log scale because of the default log link with the Poisson model.

If you really want a difference of counts, you need to regrid before doing the contrasts, because then we will have means already on the count scale before they are contrasted:

EMM |> regrid(transform = "response")  |>    # now we have counts
       contrast(method = "revpairwise")  |>
       confint()

But an additional approach is

EMM |> contrast(method = "revpairwise")  |>    # yielding differences of log(mean)s
       confint(type = "response")

In this code, we note that the difference of logs is the log of a ratio, and type = "response" converts these to ratios of mean counts. This latter approach only works with log (and logit) links, because for other than log transformations, there is no transformation family that arises from differences of transformed means.

rvlenth commented 1 year ago

I think this is resolved, so am closing this issue.