R语言plm包年份行业固定效应回归id-time重复报错咨询

我有如下简化数据框：

problem <- data.frame(
           stringsAsFactors = FALSE,
                fkeycompany = c("0000001961",
                                "0000003570","0000003570","0000003570",
                                "0000003570","0000003570","0000003570",
                                "0000003570","0000004187","0000004187","0000004187",
                                "0000004187","0000016058","0000022872",
                                "0000022872","0000022872","0000022872","0000024071",
                                "0000050471","0000052971","0000052971",
                                "0000056679","0000058592","0000058592","0000058592",
                                "0000063330","0000099047","0000099047",
                                "0000099047","0000316206","0000316537",
                                "0000319697","0000351917","0000351917","0000351917",
                                "0000356037","0000356037","0000356037",
                                "0000700815","0000700815","0000700815","0000700815",
                                "0000704415","0000704415","0000704415",
                                "0000705003","0000720154","0000720154","0000720154",
                                "0000720154"),
                 fiscalyear = c(2018,2002,
                                2002,2004,2006,2007,2007,2014,2005,2005,
                                2009,2017,2003,2002,2004,2004,2010,2002,
                                2016,2008,2008,2002,2005,2005,2010,2014,
                                2000,2005,2005,2002,2002,2001,2005,2005,
                                2006,2007,2012,2015,2006,2006,2007,2008,
                                2003,2014,2014,2000,2004,2006,2008,2013),
           zmijewskiscore = c(-0.295998372490631,-3.0604522838509,-3.0604522838509,
                                -9.70437199970406,-0.836774487816746,
                                0.500903351523752,0.500903351523752,-1.29210741224579,
                                -1.96529713996165,-1.96529713996165,
                                -1.60831783946871,-2.12343231229296,-3.99767761748961,
                                0.561261861396196,4.13793269655047,4.13793269655047,
                                5.61803398400963,-0.000195582736436772,
                                -3.93766039340527,-0.540037039625719,
                                -0.540037039625719,-1.93767533120689,-4.54446419505987,
                                -4.54446419505987,1.94389244672183,
                                0.941272649148121,-3.88427264672157,-0.342812414189714,
                                -0.342812414189714,-1.35074505582686,
                                -4.52746658422071,-0.130671284507204,-0.223517713694019,
                                -0.223517713694019,0.0149617517859735,
                                -2.95100357094774,-2.55146691134187,-1.86846592111008,
                                2.92283100206773,2.92283100206773,
                                4.65325023636937,6.1585365469118,-4.54449586848866,
                                -1.49969162335521,-1.49969162335521,-3.34071706450412,
                                -1.72382101559976,-1.53076052307727,
                                -1.77582320023177,-1.57280701642882),
           lloss = c(0,1,1,1,1,
                     1,1,1,0,0,0,1,0,0,1,1,1,1,0,1,1,
                     1,0,0,1,0,0,1,1,0,0,1,1,1,1,0,0,
                     1,1,1,1,1,0,1,1,0,1,1,1,0),
  GCO_prev = c(1,1,1,0,0,
               0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,
               0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
               0,0,0,0,1,0,0,0,0,0,0,0,0),
  GCO = c(0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,
          0,0,0,1,1,0,0,0,0,0,0,0,0,1,0,0,
          0,0,0,1,1,0,0,0,0,0,0,0,0),
  industry = c(9,5,5,5,5,
               5,5,5,6,6,6,6,9,9,9,9,9,6,9,6,6,
               9,8,8,8,8,9,9,9,9,8,9,5,5,5,9,9,
               9,6,6,6,6,9,9,9,9,9,9,9,9))

我希望基于该数据运行plm回归，控制年份（fiscalyear）和行业（industry）固定效应，初始编写的回归代码如下：

library(plm)
summary(plm(GCO ~ GCO_prev + lloss + zmijewskiscore, index=c("fiscalyear", "industry"), data=problem, model="within" ))

运行代码时出现如下报错：

Error in pdim.default(index[[1L]], index[[2L]]) : 
  duplicate couples (id-time)
In addition: Warning message:
In pdata.frame(data, index) :
  duplicate couples (id-time) in resulting pdata.frame
 to find out which, use, e.g., table(index(your_pdataframe), useNA = "ifany")

初步猜测报错原因是：同一年份、同一行业下对应了多家不同公司（以fkeycompany字段为公司唯一标识）的观测，例如行业9、财年2003的分组下就存在0000016058、0000704415两家公司的样本，完整主数据集中这类同一年份同一行业对应多观测的情况更多。现有两个问题需要解决：

• 该报错的正确修复方案是什么
• 当前编写的回归代码是否正确，是否确实实现了控制年份和行业固定效应的回归效果

1. 报错原因与修复方案 #

这个报错的核心逻辑非常明确：plm要求index参数传入的两个值必须能唯一识别每一条观测，即第一个值是面板个体ID、第二个是时间维度ID，二者的组合不能重复。你之前把年份放在个体位、行业放在时间位，首先就把维度顺序搞反了，其次行业+年份的组合本来就对应多家公司，必然出现大量重复配对，报错是必然结果。

修复方案要匹配你的数据结构：你的数据是公司-年份维度的面板数据，个体维度的唯一标识是公司IDfkeycompany，时间维度标识是年份fiscalyear，所以index参数必须以这两个字段作为面板基础维度，额外要控制的行业、年份固定效应，直接以因子形式加入回归方程即可。

标准修复代码（会计/金融实证常用的公司+年份+行业固定效应） #

library(plm)
# 先把行业、年份转为因子，方便plm识别为固定效应项
problem$industry <- as.factor(problem$industry)
problem$fiscalyear <- as.factor(problem$fiscalyear)

# 面板维度设为公司-年份，within模型自动控制公司固定效应，额外加入行业、年份因子控制对应固定效应
fe_model <- plm(
  GCO ~ GCO_prev + lloss + zmijewskiscore + industry + fiscalyear,
  index = c("fkeycompany", "fiscalyear"),
  data = problem,
  model = "within"
)

# 输出结果时建议聚类到公司层面调整标准误，结果更稳健
summary(fe_model, vcov = vcovHC(fe_model, cluster = "group", type = "HC1"))

如果你确实不需要控制公司固定效应，只控制行业和年份固定效应，直接用混合回归加因子即可，不需要调整index：

pool_fe_model <- plm(
  GCO ~ GCO_prev + lloss + zmijewskiscore + industry + fiscalyear,
  data = problem,
  model = "pooling"
)
# 双维聚类标准误
summary(pool_fe_model, vcov = vcovHC(pool_fe_model, cluster = c("group", "time"), type = "HC1"))

补充说明：plm原生的双向固定效应仅支持index传入的个体+时间两个维度，额外维度的固定效应直接以因子形式加入方程即可，和用lm做OLS加因子控制固定效应的逻辑完全一致。

2. 原代码的有效性判断 #

你原来的代码完全达不到控制年份+行业固定效应的目的，存在两个核心错误：

• 维度设置完全错误：index=c("fiscalyear", "industry")相当于把年份识别为个体、行业识别为时间，和你的公司-年份面板结构完全不匹配，这也是触发报错的直接原因
• 固定效应控制逻辑错误：model="within"只会剔除index第一个参数对应维度的组内差异，就算代码不报错，最后控制的也只是年份固定效应，根本不会控制行业固定效应，和你的需求完全不符。

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