R中非平衡面板数据绘制滞后变量散点图及年份固定效应回归方法

你需要先把长面板数据转换为匹配2014年tfp和2015年c_sales的截面数据集，再调用ggplot绘图，完整代码如下：

# 加载依赖包
library(tidyverse)

# 步骤1：按企业ID匹配两年变量
plot_data <- df %>%
  # 提取2014年tfp数据
  filter(year == 2014) %>% select(id, tfp) %>% rename(tfp_2014 = tfp) %>%
  # 关联2015年销售额数据
  inner_join(df %>% filter(year == 2015) %>% select(id, c_sales) %>% rename(c_sales_2015 = c_sales),
             by = "id") %>%
  na.omit() # 剔除仅存在单年数据的企业

# 步骤2：绘制散点图
ggplot(plot_data, aes(x = tfp_2014, y = c_sales_2015)) +
  # 散点设置，alpha调整透明度避免样本重叠
  geom_point(alpha = 0.5, size = 1.8, color = "#2c3e50") +
  # 可选添加线性拟合线与95%置信区间
  geom_smooth(method = "lm", se = TRUE, color = "#e74c3c", linewidth = 1) +
  # 轴标签与主题设置
  labs(x = "2014年企业全要素生产率（tfp）",
       y = "2015年企业销售额（c_sales）",
       title = "2014年生产率与次年销售额相关性") +
  theme_bw()

如果样本量过大重叠严重，可以把geom_point替换为geom_jitter(width = 0.02, height = 0.02)添加微小随机扰动。

如果你要检验全样本范围内滞后1期生产率对当期销售额的影响并控制年份固定效应，先构造滞后变量再执行回归即可，两种常用实现方式如下：

2.1 数据预处理 #

# 生成企业层面滞后1期的tfp变量
df <- df %>%
  group_by(id) %>%
  arrange(year) %>%
  mutate(lag_tfp = lag(tfp, 1)) %>% # 取上一年的tfp值
  ungroup() %>%
  drop_na(lag_tfp, c_sales, year)

2.2 回归操作 #

方法1：lm手动控制年份虚拟变量 #

fe_reg1 <- lm(c_sales ~ lag_tfp + factor(year), data = df)
summary(fe_reg1)
# 输出结果中lag_tfp的系数就是滞后生产率对销售额的平均影响，factor(year)对应的各组为年份固定效应估计值

方法2：plm包实现面板固定效应 #

library(plm)
# 声明面板数据结构，个体维度为id，时间维度为year
p_df <- pdata.frame(df, index = c("id", "year"))
# effect参数设置为"time"即仅控制年份固定效应，需同时控制企业个体固定效应则改为"twoways"
fe_reg2 <- plm(c_sales ~ lag_tfp, data = p_df, effect = "time", model = "within")
summary(fe_reg2)

如果你仅需要检验2014年tfp对2015年销售额的单期关系，直接用第一部分的plot_data做普通OLS回归即可，单期截面数据无需额外控制年份固定效应。

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