如何为面板数据模型拟合/实际值数据框添加日期变量并绘制时序图

Hey there! Let's tackle your problem step by step. You've run fixed, random, and first-difference panel models with plm, generated result data frames with actual/fitted values, but you're missing the YEAR variable—and since the first-difference model drops the first observation per unit, the row counts don't match up. Here's a reusable, robust way to fix this and create your time-series plots:

核心思路 #

The key here is that plm models retain the original panel data's index information. We can extract the UNIT and YEAR values directly from each model object, which ensures we match the correct dates to each observation (even for the first-difference model, which automatically removes the first time period per unit).

完整可复用代码 #

library(plm)
library(dplyr)
library(ggplot2)

# --------------- 原始数据生成（保留你的代码并修正细节）---------------
YEAR <- c(2015, 2016, 2017, 2018, 2015, 2016, 2017, 2018, 2015, 2016, 2017, 2018, 2015, 2016, 2017, 2018)
# 修正日期格式：单独年份无法直接转Date，转成当年1月1日的标准日期格式（月度数据同理可调整）
YEAR <- as.Date(ISOdate(YEAR, 1, 1))
UNIT <- c("A", "A", "A", "A", "B", "B", "B", "B", "C", "C", "C", "C", "D", "D", "D", "D")
set.seed(123) # 设置种子保证结果可复现
Y <- sample(100:1000, 16)
X1 <- sample(10:50, 16)
X2 <- sample(20:60, 16)
data <- data.frame(YEAR, UNIT, Y, X1, X2) # 补充X1/X2到原始数据，避免模型报错
crime.p <- pdata.frame(data, index=c("UNIT","YEAR"))

# --------------- 模型估计（保留你的代码）---------------
fixedeff <- plm(log(Y)~X1 + X2, data=crime.p, model="within")
randomeff <- plm(log(Y)~X1 +X2, data=crime.p, model="random")
firstdiff <- plm(log(Y)~X1 + X2, data=crime.p, model="fd")

# --------------- 提取结果并关联日期/单位信息（关键通用化修改）---------------
# 定义通用函数：输入plm模型和模型名称，返回带日期的结果数据框
get_model_results <- function(model, model_name) {
  # 提取模型对应的面板索引（自动匹配模型保留的观测）
  model_index <- as.data.frame(index(model))
  # 生成结果数据框并合并索引
  results <- data.frame(
    fitted = predict(model),
    residuals = model$residuals
  ) %>%
    mutate(
      actual = fitted + residuals,
      model = model_name
    ) %>%
    bind_cols(model_index) %>% # 合并单位和日期变量
    select(UNIT, YEAR, actual, fitted, residuals, model)
  return(results)
}

# 批量生成三个模型的结果
fixx_results <- get_model_results(fixedeff, "fixed")
random_results <- get_model_results(randomeff, "random")
fd_results <- get_model_results(firstdiff, "fd")

# 合并所有结果
fitted_res_all <- bind_rows(fixx_results, random_results, fd_results)

# --------------- 绘制时序图（支持年度/月度数据）---------------
# 按单位和模型分组，绘制实际值vs拟合值的时序对比图
ggplot(fitted_res_all, aes(x = YEAR)) +
  geom_line(aes(y = actual, color = "Actual"), linewidth = 1) +
  geom_line(aes(y = fitted, color = "Fitted"), linetype = "dashed", linewidth = 1) +
  facet_grid(model ~ UNIT) + # 按模型和单位分面展示
  scale_color_manual(values = c("Actual" = "black", "Fitted" = "blue")) +
  labs(title = "Panel Model Actual vs Fitted Values Over Time",
       x = "Year", y = "Log(Y)", color = "Value Type") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

关键细节解释 #

• 日期格式修正：你的原始代码中as.Date(YEAR)会报错，因为单独的年份不是有效的日期格式。用ISOdate(YEAR, 1, 1)可以将年份转为当年1月1日的标准日期，这个逻辑也适用于月度数据（只需传入年月参数即可）。
• 通用函数get_model_results：这个函数可以复用在任何plm模型上，自动提取对应观测的UNIT和YEAR，完美解决不同模型观测数不一致的问题（比如fd模型自动排除了每个单位的第一期数据，索引会自动匹配剩余的有效观测）。
• 绘图灵活性：用ggplot2的分面功能可以同时查看不同单位、不同模型的时序对比，你也可以根据需求调整（比如去掉分面，只查看单个单位或单个模型的结果）。

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