如何绘制含交互项的双向固定效应回归的边际效应？——feols回归后iplot绘图异常问题咨询

It looks like the issue with iplot(tal_lpm4) is that it's defaulting to showing the interaction term coefficients rather than the total marginal effect of provtariff for each sex group. Let's break down why this happens and how to fix it properly.

Why Your Current iplot() Isn't Working #

Your model includes both provtariff*i(sex) and an explicit i(sex) term, which creates redundant parameterization (the interaction already includes the sex main effect). More importantly, iplot() in fixest by default visualizes the coefficients of the interaction terms themselves—not the combined effect of the provtariff main effect plus the interaction term for each group.

Looking at your regression results:

• For males, the true marginal effect of provtariff is just the main provtariff coefficient (0.200287).
• For females, it's the main provtariff coefficient plus the provtariff:sex::Female interaction term (0.200287 + (-0.216309) = -0.016022).

iplot() isn't combining these values automatically, hence the incorrect plot.

Solution 1: Use the marginaleffects Package (Recommended) #

The marginaleffects package is designed explicitly for calculating and visualizing marginal effects, and it works seamlessly with fixest models. Here's how to use it:

1. Install and load the package:

install.packages("marginaleffects")
library(marginaleffects)

2. Calculate sex-specific marginal effects for provtariff:

# Compute marginal effects grouped by sex
mfx <- marginaleffects(tal_lpm4, 
                       variables = "provtariff", 
                       by = "sex")

3. Plot the results (this will show the total marginal effect + confidence intervals for each sex):

plot(mfx)

This will automatically handle the coefficient combination and correct standard error calculation (no manual math needed!).

Solution 2: Simplify Your Model and Use fixest's Built-in Tools #

First, simplify your model to remove redundant terms—provtariff*sex already includes both the main effects of sex and the interaction, so you don't need i(sex) separately:

tal_lpm4_clean <- feols(
  tal ~ provtariff*sex + as.numeric(educ) + age + I(age^2) | year + tinh,
  data = employment0204,
  vcov_cluster(~year + tinh),
  weights = ~hhwt
)

Then, use fixest's get_margins() to extract the group-specific marginal effects, followed by iplot():

# Extract marginal effects for provtariff by sex
margins <- get_margins(tal_lpm4_clean, vars = "provtariff", by = "sex")

# Plot the marginal effects
iplot(margins)

Solution 3: Manual Calculation + ggplot (For Full Control) #

If you want complete control over the plot, you can manually compute the marginal effects and their standard errors, then use ggplot2:

library(ggplot2)

# Extract coefficients and variance-covariance matrix
coefs <- coef(tal_lpm4)
vcov_mat <- vcov(tal_lpm4)

# Calculate marginal effects
male_effect <- coefs["provtariff"]
female_effect <- coefs["provtariff"] + coefs["provtariff:sex::Female"]

# Calculate standard errors (using variance rules for sums)
male_se <- sqrt(vcov_mat["provtariff", "provtariff"])
female_se <- sqrt(
  vcov_mat["provtariff", "provtariff"] +
  vcov_mat["provtariff:sex::Female", "provtariff:sex::Female"] +
  2*vcov_mat["provtariff", "provtariff:sex::Female"]
)

# Create a data frame for plotting
mfx_data <- data.frame(
  sex = c("Male", "Female"),
  estimate = c(male_effect, female_effect),
  se = c(male_se, female_se)
)

# Plot with ggplot2
ggplot(mfx_data, aes(x = sex, y = estimate)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = estimate - 1.96*se, ymax = estimate + 1.96*se), width = 0.2) +
  labs(
    x = "Sex",
    y = "Marginal Effect of provtariff",
    title = "Sex-Specific Marginal Effects of provtariff on tal"
  ) +
  theme_minimal()

All three methods will give you the correct sex-specific marginal effects you're looking for. The marginaleffects approach is the most robust and least error-prone.

内容的提问来源于stack exchange，提问作者anrisakaki96