使用R中fepois进行泊松事件研究时仅得单个系数的问题

我尝试用R的fepois函数做泊松事件研究，分析自行车道安装对汽车事故的影响，但模型结果始终只返回事件时间的单个系数，而非每个时间单位对应一个系数。

模型运行代码及输出：

> model_poisson <- fepois(
+     accident_count ~ i(rel_time_binned, ref = -1) | segment_id + month,
+     data = df_panel_reg,
+     cluster = ~segment_id  # Cluster standard errors at segment level
+ )
> 
> summary(model_poisson)
Coefficients: 
                 Estimate Std. Error t value Pr(>|t|)    
rel_time_binned -0.017631   0.001332  -13.23   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

可复现的示例代码：

library(tidyverse)
library(lfe)

set.seed(123)
n_segments <- 100
n_periods <- 60  # months

# Create panel data
df_panel <- expand.grid(
  segment_id = 1:n_segments,
  month = 1:n_periods
) %>%
  mutate(
    # Random installation dates (some segments never treated)
    installation_month = sample(c(NA, 20:50), n_segments, replace = TRUE)[segment_id],
    # Calculate relative time to treatment
    rel_time = month - installation_month,
    # Generate accident counts (Poisson with treatment effect)
    treated = !is.na(installation_month) & month >= installation_month,
    lambda = exp(2.5 - 0.3 * treated + rnorm(n(), 0, 0.1)),
    accident_count = rpois(n(), lambda),
    # Convert segment_id to character to avoid issues
    segment_id = as.character(segment_id)
  )

# Define event time bins
create_event_time_dummies <- function(data, max_lag = 10, max_lead = 10) {
  data %>%
    mutate(
      # Create relative time variable
      rel_time = ifelse(is.na(installation_month), NA, month - installation_month),
      
      # Bin endpoints
      rel_time_binned = case_when(
        is.na(rel_time) ~ NA_real_,  # Never treated
        rel_time < -max_lead ~ -max_lead,  # Bin distant pre-periods
        rel_time > max_lag ~ max_lag,      # Bin distant post-periods
        TRUE ~ rel_time
      ),
      
      # Convert to integer for fepois
      rel_time_binned = as.integer(rel_time_binned),
  
    )
}

df_panel <- create_event_time_dummies(df_panel, max_lag = 10, max_lead = 10)

df_panel_reg <- df_panel %>%
  filter(!is.na(rel_time_binned))

model_poisson <- fepois(
  accident_count ~ i(rel_time_binned, ref = -1) | segment_id + month,
  data = df_panel_reg,
  cluster = ~segment_id  # Cluster standard errors at segment level
)

summary(model_poisson)

问题原因及解决方法 #

问题出在rel_time_binned的变量类型上：你将其转换为整数，但lfe包中的i()函数需要变量是因子类型，才能为每个时间区间生成对应的虚拟变量。当变量是整数时，i()会将其视为连续变量处理，因此仅返回一个系数。

修改步骤 #

1. 在create_event_time_dummies函数中，将rel_time_binned转换为因子，并明确指定所有可能的水平（从-max_lead到max_lag），确保每个时间区间都被识别为独立类别：

rel_time_binned = factor(rel_time_binned, levels = -max_lead:max_lag)

2. 重新生成数据并运行模型，此时i()函数会自动为每个因子水平生成虚拟变量，排除指定的参考水平-1后，就能得到每个时间单位对应的系数。

修改后的完整create_event_time_dummies函数：

create_event_time_dummies <- function(data, max_lag = 10, max_lead = 10) {
  data %>%
    mutate(
      rel_time = ifelse(is.na(installation_month), NA, month - installation_month),
      rel_time_binned = case_when(
        is.na(rel_time) ~ NA_real_,
        rel_time < -max_lead ~ -max_lead,
        rel_time > max_lag ~ max_lag,
        TRUE ~ rel_time
      ),
      # 转换为因子并指定完整水平范围
      rel_time_binned = factor(rel_time_binned, levels = -max_lead:max_lag)
    )
}

重新运行模型后，summary(model_poisson)会输出每个相对时间区间的系数，清晰展示不同时间点的处理效应。

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