1

I viewed lots of relevant bloggers who explained how to generate marginal effects, but somehow most of them did not take interaction effects as sample.

When I use Stata to do it, the margin effects are quite easy to compute, but R seems in another way and show different results.

I want to know how to compute proper marginal effects with interaction properly,

  1. why function margins() does not show the interaction factor?
  2. if I compute this by hand coef(glm)*mean(dlogis(predict(glm, type = "link")),na.rm=T) ,it looks more close to the Stata result, not with the function margins()?
  3. what's more, how can I compute the standard error of marginal effects handly, will it like this sqrt(diag(vcov(glm1)))*mean(dlogis(predict(glm, type = "link")),na.rm=T) ?
trydata=read.csv("https://vincentarelbundock.github.io/Rdatasets/csv/palmerpenguins/penguins.csv")
trydata[c("species","sex")]=lapply(trydata[c("species","sex" )], as.factor)
contrasts(trydata$species)=contr.treatment(3) 
contrasts(trydata$sex)=contr.treatment(2)

try=glm(bill_length_mm ~sex+species+bill_depth_mm + sex*species, data=trydata)
summary(try)
summary(margins(try))

some links I saw before and thought useful:

code: marginal effects with interaction

Marginal effects calculated by hand, but standard error doesn't show

Compare R vs STATA with various packages and functions: Interpreting Regression Results using Average Marginal Effects with R’s margins

Probit/Logit Marginal Effects in R

An Introduction to ‘margins’

Nick Cox
  • 35,529
  • 6
  • 31
  • 47
vlamenco
  • 37
  • 6

2 Answers2

4

I think margins() is trying to prevent you from doing something that doesn't make sense. There are a couple of reasons for the difference.

The most important is that when your x-variable (e.g., sex) is discrete (as it is in this case), the partial first derivative with respect to x is not meaningful because the idea of inducing an infinitesimal change to x doesn't make sense. In your data, sex can be either male or female. So, what margins() will do is calculate a first difference rather than a first derivative.

The other reason is that it is similarly not meaningful to think about the partial first derivative of the interaction term - the effect of changing the interaction term while holding the other terms constant. You can't change the interaction term while holding the constitutive terms fixed.

So, I am not sure how you would want to interpret a marginal effect for just the interaction term itself.

As for the third question, you would have to use the delta method or a parametric bootstrap to calculate the standard error of the partial first derivative since the partial first derivative is a non-linear combination of the coefficients.

Nick Cox
  • 35,529
  • 6
  • 31
  • 47
DaveArmstrong
  • 18,377
  • 2
  • 13
  • 25
  • I upvoted this answer, but I wonder if it shouldn't be on https://stats.stackexchange.com ? – IRTFM Feb 10 '22 at 21:27
  • Great thanks for you! but to be honest, I'm not good at math, too hard for me to understand the whole thing~ – vlamenco Feb 10 '22 at 21:36
  • 1
    @IRTFM it could be there, but this is also about why R is doing something in particular. Perhaps it could go either way. – DaveArmstrong Feb 10 '22 at 22:29
  • 1
    Sometimes the builders of R omit calculations they think are of questionable validity. Perhaps the real statisticians would have perspective. – IRTFM Feb 11 '22 at 04:28
3

You can use the marginaleffects package to do this (disclaimer: I am the author). For discrete regressors, the function will report pairwise contrasts.

Load data and package, and estimate the model:

library(marginaleffects)

trydata = read.csv("https://vincentarelbundock.github.io/Rdatasets/csv/palmerpenguins/penguins.csv")
trydata[c("species","sex")] = lapply(trydata[c("species","sex" )], as.factor)
contrasts(trydata$species) = contr.treatment(3)
contrasts(trydata$sex) = contr.treatment(2)

try <- glm(bill_length_mm ~ sex + species + bill_depth_mm + sex * species, data = trydata)

Average marginaleffects/contrasts:

mfx <- marginaleffects(try)
summary(mfx)

## Average marginal effects 
##            Term           Contrast  Effect Std. Error z value   Pr(>|z|) 2.5 %
## 1 bill_depth_mm                     0.5189     0.1510   3.436 0.00058933 0.223
## 2           sex      male - female  2.9130     0.3376   8.628 < 2.22e-16 2.251
## 3       species Chinstrap - Adelie  9.9775     0.3347  29.812 < 2.22e-16 9.322
## 4       species    Gentoo - Adelie 10.4467     0.5814  17.968 < 2.22e-16 9.307
##    97.5 %
## 1  0.8149
## 2  3.5748
## 3 10.6335
## 4 11.5862
## 
## Model type:  glm 
## Prediction type:  response

Marginal Effects at Representative values, using the newdata argument and the datagrid function:

marginaleffects(try,
                variables = "species",
                newdata = datagrid(sex = c("male", "female"), species = trydata$species))

##    rowid     type    term           contrast      dydx std.error bill_depth_mm
## 1      1 response species Chinstrap - Adelie 10.610124 0.4738125      17.16486
## 2      2 response species Chinstrap - Adelie  9.333474 0.4730566      17.16486
## 3      3 response species Chinstrap - Adelie 10.610124 0.4738125      17.16486
## 4      4 response species Chinstrap - Adelie  9.333474 0.4730566      17.16486
## 5      5 response species Chinstrap - Adelie 10.610124 0.4738125      17.16486
## 6      6 response species Chinstrap - Adelie  9.333474 0.4730566      17.16486
## 7      1 response species    Gentoo - Adelie 10.824147 0.6424885      17.16486
## 8      2 response species    Gentoo - Adelie 10.062312 0.6493940      17.16486
## 9      3 response species    Gentoo - Adelie 10.824147 0.6424885      17.16486
## 10     4 response species    Gentoo - Adelie 10.062312 0.6493940      17.16486
## 11     5 response species    Gentoo - Adelie 10.824147 0.6424885      17.16486
## 12     6 response species    Gentoo - Adelie 10.062312 0.6493940      17.16486
Vincent
  • 15,809
  • 7
  • 37
  • 39