如何基于PLM固定效应模型系数绘制时间维度X-Y关系图？

完全可以实现，而且不用手动收集系数，有更高效的自动化方法，以下是具体步骤：

核心思路调整 #

你当前的模型仅加入了as.factor(Time)作为控制项，这只能控制时间固定效应，无法直接得到父母财富（X）对子女财富（Y）的边际效应随时间的变化。要实现这个目标，需要给X和时间虚拟变量添加交互项，让模型估计每个年度X的单独系数。

步骤1：拟合带交互项的PLM模型 #

直接用公式X*as.factor(Time)自动生成X的主效应、时间虚拟变量主效应，以及二者的交互项：

library(plm)
# 拟合模型
model <- plm(Y ~ X*as.factor(Time) + a + b, 
             data = df, 
             index = c("ID", "Time"), 
             model = "within", 
             effect = "individual")
summary(model)

模型输出中，形如X:as.factor(Time)XXXX的系数就是对应年度下X对Y的边际效应。

步骤2：提取并整理年度系数 #

从模型结果中筛选出交互项系数，整理成适合绘图的数据框：

# 提取所有系数
coefs <- coef(model)
# 筛选X与时间交互的系数
time_specific_coefs <- coefs[grepl("X:as.factor\\(Time\\)", names(coefs))]
# 整理成绘图数据：提取年度数值，对应系数
plot_df <- data.frame(
  Year = as.numeric(gsub("X:as.factor\\(Time\\)", "", names(time_specific_coefs))),
  Marginal_Effect = as.numeric(time_specific_coefs)
)

步骤3：绘制边际效应随时间变化图 #

用ggplot2绘制折线图，直观展示年度变化趋势：

library(ggplot2)
ggplot(plot_df, aes(x = Year, y = Marginal_Effect)) +
  geom_line(color = "#2E86AB", linewidth = 1.2) +
  geom_point(size = 3, fill = "white", shape = 21, color = "#2E86AB") +
  labs(title = "父母财富对子女财富的边际效应年度变化",
       x = "年度", y = "边际效应系数") +
  theme_minimal()

可选：绘制分年度X-Y拟合线 #

如果需要直接展示每个年度X与Y的散点+拟合线，可以按年度分组拟合线性模型后绘图：

library(dplyr)
# 按年度分组拟合控制a、b后的线性模型
yearly_fits <- df %>%
  group_by(Time) %>%
  do(fit = lm(Y ~ X + a + b, data = .)) %>%
  mutate(
    slope_X = coef(fit)[["X"]],
    intercept = coef(fit)[["(Intercept)"]]
  )
# 绘图
ggplot(df, aes(x = X, y = Y)) +
  geom_point(alpha = 0.2, color = "#999999") +
  geom_abline(data = yearly_fits, 
              aes(intercept = intercept, slope = slope_X, color = as.factor(Time)),
              linewidth = 1) +
  labs(title = "分年度父母财富与子女财富的关系",
       x = "父母财富", y = "子女财富", color = "年度") +
  theme_minimal()

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