roffsetpoissonbrmsr-marginaleffects

Calculating the average marginal effect of a single day of exposure to a treatment from a Poisson GLMM with an offset


I am trying to calculate the average marginal effect of a treatment (i.e., the average change in the response variable when participants are exposed to the treatment) per day of treatment exposure (which varies between participants).

I'm fitting the following Poisson regression with an offset in brms:

library(brms)
library(marginaleffects)
set.seed(42) # For reproducibility

# Number of observations
n <- 100

# Generate the data
DF_anly_1 <- data.frame(
  catch_count = rpois(n, lambda = 2), # Simulated count data using a Poisson distribution
  Treatment = sample(c('Control', 'Treatment'), n, replace = TRUE), # Treatment groups
  village = sample(c('Village1', 'Village2', 'Village3'), n, replace = TRUE), # Village categories
  phase_exp = sample(c('Phase1', 'Phase2'), n, replace = TRUE), # Experimental phases
  boat = sample(c('Boat1', 'Boat2', 'Boat3'), n, replace = TRUE), # Random effects for boats
  day_count = sample(1:30, n, replace = TRUE) # Days of count/exposure
)

# Model 
m1.WF_brms <- brm(catch_count ~ Treatment + village + phase_exp + (1|boat) + offset(log(day_count)),
              data = DF_anly_1, family = poisson(),
              chains = 4, core = 4, iter = 4000)

I am calculating the average marginal effect of the treatment using the marginaleffects package on the response scale:

# Marginal treatment effect 
AME <- comparisons(
  m1.WF_brms, 
  variables = "Treatment", 
  re_formula = NULL,
  type = "response") %>% 
  tidy()

My understanding is that this function calculates the treatment effect, holding constant the average effect of all other variables (i.e., the population-level effect of the treatment). Does this also apply to offsets? In other words, can I get the average marginal effect per day of treatment exposure by dividing the AME by the mean number of exposure days:

AME$estimate / mean(DF_anly_1$day_count)   

I have implemented the advice here (see below): https://discourse.mc-stan.org/t/estimating-marginal-effects-of-offset-in-brms-poisson-regression-model/4463/4, but this appears to give the conditional rather than average marginal effect:

marginal_effects(m1.WF_brms, 
                 conditions = data.frame(day_count = 1), 
                 effects  = "Treatment", 
                 type = "response")

Solution

  • Sorry, but I don’t understand exactly what quantity you are looking for. But I’ll give you three options with marginaleffects and show manual computations so you know exactly what happens under the hood.

    First, notice that if we call brms::posterior_epred() we get 8000 draws from the posterior distribution of predicted values (there are 100 points in your dataset):

    library(brms)
    library(marginaleffects)
    
    p <- posterior_epred(m1.WF_brms, re_formula = NULL)
    dim(p)
    # [1] 8000  100
    

    Now, for each row in the original data, we make predictions in the counterfactual worlds where Treatment is either “Control” or “Treatment”, and we take the difference between those predictions, and then take the median of the posterior draws. This gives us 100 estimates, one for each row:

    d0 <- transform(DF_anly_1, Treatment = "Control")
    d1 <- transform(DF_anly_1, Treatment = "Treatment")
    p0 <- posterior_epred(m1.WF_brms, d0, re_formula = NULL)
    p1 <- posterior_epred(m1.WF_brms, d1, re_formula = NULL)
    
    effects <- apply(p1 - p0, 2, median) 
    head(effects)
    # [1] -0.24480083 -0.14519547 -0.01207146 -0.09687572 -0.27200092 -0.29490911
    

    This is equivalent to calling comparisons():

    comparisons(
      m1.WF_brms, 
      variables = "Treatment", 
      re_formula = NULL,
      type = "response") |> head()
    # 
    #       Term            Contrast Estimate  2.5 % 97.5 %
    #  Treatment Treatment - Control  -0.2448 -1.081 0.4142
    #  Treatment Treatment - Control  -0.1452 -0.575 0.2697
    #  Treatment Treatment - Control  -0.0121 -0.054 0.0215
    #  Treatment Treatment - Control  -0.0969 -0.389 0.1745
    #  Treatment Treatment - Control  -0.2720 -1.201 0.4602
    #  Treatment Treatment - Control  -0.2949 -1.293 0.5282
    # 
    # Columns: rowid, term, contrast, estimate, conf.low, conf.high, predicted_lo, predicted_hi, predicted, tmp_idx, catch_count, Treatment, village, phase_exp, day_count, boat 
    # Type:  response
    

    Now we can take the average of those 100 estimates. This is equivalent to calling avg_comparisons(). Warning: I use the avg_ prefix because the tidy() feature that your code relies on is now deprecated in marginaleffects 0.18.0:

    median(apply(p1 - p0, 1, mean))
    # [1] -0.2159907
    
    avg_comparisons(
      m1.WF_brms, 
      variables = "Treatment", 
      re_formula = NULL,
      type = "response")
    # 
    #       Term            Contrast Estimate  2.5 % 97.5 %
    #  Treatment Treatment - Control   -0.216 -0.884  0.373
    # 
    # Columns: term, contrast, estimate, conf.low, conf.high 
    # Type:  response
    

    Finally, you can also report aggregated estimates for each day of exposure using the by argument:

    avg_comparisons(
      m1.WF_brms, 
      by = "day_count",
      variables = "Treatment", 
      re_formula = NULL,
      type = "response")
    # 
    #       Term                        Contrast day_count Estimate   2.5 % 97.5 %
    #  Treatment mean(Treatment) - mean(Control)         1  -0.0143 -0.0582 0.0253
    #  Treatment mean(Treatment) - mean(Control)         2  -0.0253 -0.1075 0.0437
    #  Treatment mean(Treatment) - mean(Control)         3  -0.0388 -0.1676 0.0663
    #  Treatment mean(Treatment) - mean(Control)         4  -0.0536 -0.2318 0.0923
    #  Treatment mean(Treatment) - mean(Control)         5  -0.0702 -0.3078 0.1258
    #  Treatment mean(Treatment) - mean(Control)         6  -0.0822 -0.3423 0.1414
    #  Treatment mean(Treatment) - mean(Control)         7  -0.1010 -0.4138 0.1789
    #  Treatment mean(Treatment) - mean(Control)         8  -0.0966 -0.3994 0.1689
    #  Treatment mean(Treatment) - mean(Control)         9  -0.1457 -0.6396 0.2570
    #  Treatment mean(Treatment) - mean(Control)        10  -0.1409 -0.5689 0.2489
    #  Treatment mean(Treatment) - mean(Control)        11  -0.1428 -0.6072 0.2453
    #  Treatment mean(Treatment) - mean(Control)        12  -0.1599 -0.6642 0.2801
    #  Treatment mean(Treatment) - mean(Control)        13  -0.1775 -0.7256 0.3081
    #  Treatment mean(Treatment) - mean(Control)        14  -0.1822 -0.7467 0.3196
    #  Treatment mean(Treatment) - mean(Control)        15  -0.1954 -0.7978 0.3444
    #  Treatment mean(Treatment) - mean(Control)        16  -0.1968 -0.8214 0.3402
    #  Treatment mean(Treatment) - mean(Control)        17  -0.2415 -0.9836 0.4246
    #  Treatment mean(Treatment) - mean(Control)        18  -0.2448 -1.0810 0.4142
    #  Treatment mean(Treatment) - mean(Control)        19  -0.2328 -0.9672 0.3940
    #  Treatment mean(Treatment) - mean(Control)        20  -0.2533 -1.0799 0.4302
    #  Treatment mean(Treatment) - mean(Control)        21  -0.2893 -1.1689 0.5012
    #  Treatment mean(Treatment) - mean(Control)        22  -0.3194 -1.2644 0.5934
    #  Treatment mean(Treatment) - mean(Control)        23  -0.3380 -1.4055 0.5845
    #  Treatment mean(Treatment) - mean(Control)        24  -0.3073 -1.2686 0.5329
    #  Treatment mean(Treatment) - mean(Control)        25  -0.3229 -1.2969 0.5640
    #  Treatment mean(Treatment) - mean(Control)        26  -0.3215 -1.3020 0.5658
    #  Treatment mean(Treatment) - mean(Control)        27  -0.3411 -1.3967 0.5950
    #  Treatment mean(Treatment) - mean(Control)        28  -0.3539 -1.4158 0.6239
    #  Treatment mean(Treatment) - mean(Control)        29  -0.3764 -1.5477 0.6613
    #  Treatment mean(Treatment) - mean(Control)        30  -0.4241 -1.7160 0.7473
    # 
    # Columns: term, contrast, day_count, estimate, conf.low, conf.high, predicted_lo, predicted_hi, predicted, tmp_idx 
    # Type:  response