Plotting dynamic treatment effects using marginaleffects package in R

Question

I am attempting to use the marginaleffects package in R to plot the effect of treatment on the outcome across time (i.e., dynamic treatment effects). My data consists of two groups (treated and control), seven time points, and a continuous outcome.

library(tibble)

d <- tibble(id = rep(1:10, each = 7),
            periods = factor(rep(0:6, times = 10)),
            treated = factor(ifelse(id<=5, 1, 0)),
            y = ifelse(treated==1, rnorm(100, 2, 0.5), rnorm(100, 0, 1)))

head(d, n=10)
# A tibble: 10 × 4
      id periods treated     y
   <int> <fct>   <fct>   <dbl>
 1     1 0       1        2.31
 2     1 1       1        2.32
 3     1 2       1        2.16
 4     1 3       1        1.74
 5     1 4       1        1.93
 6     1 5       1        2.17
 7     1 6       1        2.25
 8     2 0       1        2.02
 9     2 1       1        1.90
10     2 2       1        2.50

To fit the model, I used a standard two-way fixed-effects DiD approach where there are fixed effects for time (periods in my case) and group (treated in my case), and the treatment effect is represented by the interaction of the two (treated*periods).

In the model output (printed below), I am interested in plotting the effects for the treated*:periods* terms--which represent the effect of treatment across time after accounting for the group and time fixed effects. 

# fit model
m <- lm(y ~ treated*periods, data = d)

# print model output
tidy(m, conf.int = TRUE)

# A tibble: 14 × 7
   term              estimate std.error statistic   p.value conf.low conf.high
   <chr>                <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
 1 (Intercept)         0.269      0.308    0.872  0.387       -0.349     0.887
 2 treated1            2.11       0.436    4.84   0.0000107    1.24      2.98 
 3 periods1            0.475      0.436    1.09   0.281       -0.399     1.35 
 4 periods2           -0.0516     0.436   -0.118  0.906       -0.925     0.822
 5 periods3           -0.343      0.436   -0.787  0.435       -1.22      0.531
 6 periods4            0.0433     0.436    0.0994 0.921       -0.831     0.917
 7 periods5            0.578      0.436    1.33   0.190       -0.296     1.45 
 8 periods6           -0.220      0.436   -0.504  0.616       -1.09      0.654
 9 treated1:periods1  -1.07       0.617   -1.73   0.0897      -2.30      0.171
10 treated1:periods2   0.171      0.617    0.277  0.783       -1.07      1.41 
11 treated1:periods3   0.0223     0.617    0.0361 0.971       -1.21      1.26 
12 treated1:periods4  -0.635      0.617   -1.03   0.308       -1.87      0.601
13 treated1:periods5  -0.966      0.617   -1.57   0.123       -2.20      0.270
14 treated1:periods6  -0.170      0.617   -0.275  0.784       -1.41      1.07

However, when I attempt to plot the treatment effects using marginaleffects, the plotted values do not match the treatment effects from the model output (the treated*:periods* terms). I included the output of the call to avg_comparisons() rather than plot_comparisons() below for ease of reference:

library(marginaleffects)

avg_comparisons(m, variables = "treated", by = "periods")

    Term          Contrast periods Estimate Std. Error    z Pr(>|z|)    S 2.5 % 97.5 %
 treated mean(1) - mean(0)       0     2.11      0.436 4.84  < 0.001 19.5  1.26   2.97
 treated mean(1) - mean(0)       1     1.05      0.436 2.40  0.01657  5.9  0.19   1.90
 treated mean(1) - mean(0)       2     2.28      0.436 5.23  < 0.001 22.5  1.43   3.14
 treated mean(1) - mean(0)       3     2.13      0.436 4.89  < 0.001 19.9  1.28   2.99
 treated mean(1) - mean(0)       4     1.48      0.436 3.38  < 0.001 10.4  0.62   2.33
 treated mean(1) - mean(0)       5     1.14      0.436 2.62  0.00868  6.8  0.29   2.00
 treated mean(1) - mean(0)       6     1.94      0.436 4.45  < 0.001 16.8  1.09   2.80

It looks to me like avg_comparisons() is not taking the fixed effects into account when estimating the comparisons between the groups over time. How do I obtain avg_comparisons() or plot_comparisons() that match the treatment effects contained in the model output (treated*:periods*)?

I have attempted numerous variants of this call avg_comparisons(m, variables = c("periods", "treated"), by = "treated"), including switching which terms are included in the call to variables and which variables are included in the call to by. But none of them produce the output I am looking for.

I have also attempted avg_comparisons(m, variables = "treated:periods") to obtain the interaction effect, but this produced the error: Error: There is no valid predictor variable. Please change the variables argument or supply a new data frame to the 'newdata' argument.

Answer 1

I think you are misinterpreting your model coefficients and that avg_comparisons() correctly reports the quantities of interest.

Consider this code, where I set a seed for replicability:

library(tibble)
library(marginaleffects)
set.seed(1024)

d <- tibble(
  id = rep(1:10, each = 7),
  periods = rep(0:6, times = 10),
  treated = ifelse(id <= 5, 1, 0),
  y = ifelse(treated == 1, rnorm(100, 2, 0.5), rnorm(100, 0, 1)))

m <- lm(y ~ factor(treated) * factor(periods), data = d)

Estimate the average difference in predicted outcome between the treated and control groups, at different periods:

e <- avg_comparisons(m, variables = "treated", by = "periods") 

First, note that the “effect” of a change from 0 to 1 in the treated variable in period 2 is not equal to just the interaction term. It is equal to the sum of these two coefficients:

coef(m)["factor(treated)1"] + coef(m)["factor(treated)1:factor(periods)2"]
# factor(treated)1 
#         2.153658

This can be verified using the predict() function from base R. See what happens to the predicted outcome when treated goes from 0 to 1, for different periods: 0 and 2.

d0 <- data.frame(treated = 0, periods = 0)
d1 <- data.frame(treated = 1, periods = 0)
p0 <- predict(m, newdata = d0)
p1 <- predict(m, newdata = d1)
p1 - p0
#        1 
# 2.792739

d0 <- data.frame(treated = 0, periods = 2)
d1 <- data.frame(treated = 1, periods = 2)
p0 <- predict(m, newdata = d0)
p1 <- predict(m, newdata = d1)
p1 - p0
#        1 
# 2.153658

Finally, see that these are exactly the results that avg_comparisons() reports:

subset(e, periods %in% c(0, 2))
# 
#     Term          Contrast Estimate Std. Error    z Pr(>|z|)    S CI low
#  treated mean(1) - mean(0)     2.79      0.578 4.83   <0.001 19.5   1.66
#  treated mean(1) - mean(0)     2.15      0.578 3.72   <0.001 12.3   1.02
#  CI high
#     3.93
#     3.29
# 
# Columns: term, contrast, periods, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted