r - plm:使用 fixef() 手动计算固定效应双向模型的拟合值

Question

请注意:我试图让代码同时处理时间和个人固定效应以及不平衡的数据集。下面的示例代码适用于平衡数据集。

也请参阅下面的编辑，请

我正在尝试使用 plm 手动计算固定效应模型(具有个体效应和时间效应)的拟合值包裹。这更多是为了确认我了解模型和包的机制，我知道我可以从 plm 中获得拟合值。对象，来自两个相关问题( here 和 here )。

来自 plm小插图(第 2 页)，底层模型是:

y _it = 阿尔法 + 测试版 _转置 * x _it + ( mu _i + lambda _t + epsilon _it)

其中 mu_i 是误差项的单个分量(又名“个体效应”)，而 lambda_t 是“时间效应”。

可以使用 fixef() 提取固定效应。我想我可以使用它们(与自变量一起)来计算模型的拟合值，使用(带有两个自变量)以这种方式:

适合 _it = 阿尔法 + 测试版 _1 * x1 + 测试版 _2 * x2 + 亩 _i + lambda _t

这就是我失败的地方——我得到的值远不及拟合值(我得到的是模型对象中实际值和残差之间的差异)。一方面，我看不到 alpha任何地方。我尝试将固定效果显示为与第一个、平均值等的差异，但没有成功。

我缺少什么？这很可能是对模型的误解，或者是代码中的错误，恐怕……提前致谢。

PS:其中一个相关问题暗示pmodel.response()应该与我的问题有关(以及没有 plm.fit 函数的原因)，但它的帮助页面并不能帮助我理解该函数的实际作用，而且我找不到任何示例来解释它产生的结果。

谢谢!

我所做的示例代码:

library(data.table); library(plm)

set.seed(100)
DT <- data.table(CJ(id=c("a","b","c","d"), time=c(1:10)))
DT[, x1:=rnorm(40)]
DT[, x2:=rnorm(40)]
DT[, y:=x1 + 2*x2 + rnorm(40)/10]
DT <- DT[!(id=="a" & time==4)] # just to make it an unbalanced panel
setkey(DT, id, time)

summary(plmFEit <- plm(data=DT, id=c("id","time"), formula=y ~ x1 + x2, model="within", effect="twoways"))

# Extract the fitted values from the plm object
FV <- data.table(plmFEit$model, residuals=as.numeric(plmFEit$residuals))
FV[, y := as.numeric(y)]
FV[, x1 := as.numeric(x1)]
FV[, x2 := as.numeric(x2)]

DT <- merge(x=DT, y=FV, by=c("y","x1","x2"), all=TRUE)
DT[, fitted.plm := as.numeric(y) - as.numeric(residuals)]

FEI <- data.table(as.matrix(fixef(object=plmFEit, effect="individual", type="level")), keep.rownames=TRUE) # as.matrix needed to preserve the names?
setnames(FEI, c("id","fei"))
setkey(FEI, id)
setkey(DT, id)
DT <- DT[FEI] # merge the fei into the data, each id gets a single number for every row

FET <- data.table(as.matrix(fixef(object=plmFEit, effect="time", type="level")), keep.rownames=TRUE) # as.matrix needed to preserve the names?
setnames(FET, c("time","fet"))
FET[, time := as.integer(time)] # fixef returns time as character
setkey(FET, time)
setkey(DT, time)
DT <- DT[FET] # merge the fet into the data, each time gets a single number for every row

# calculate the fitted values (called calc to distinguish from those from plm)
DT[, fitted.calc := as.numeric(coef(plmFEit)[1] * x1 + coef(plmFEit)[2]*x2 + fei + fet)]
DT[, diff := as.numeric(fitted.plm - fitted.calc)]

all.equal(DT$fitted.plm, DT$fitted.calc)

我的 session 如下:

R version 3.2.2 (2015-08-14)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 8 x64 (build 9200)

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252    LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] plm_1.4-0           Formula_1.2-1       RJSONIO_1.3-0       jsonlite_0.9.17     readxl_0.1.0.9000   data.table_1.9.7    bit64_0.9-5         bit_1.1-12          RevoUtilsMath_3.2.2

loaded via a namespace (and not attached):
 [1] bdsmatrix_1.3-2  Rcpp_0.12.1      lattice_0.20-33  zoo_1.7-12       MASS_7.3-44      grid_3.2.2       chron_2.3-47     nlme_3.1-122     curl_0.9.3       rstudioapi_0.3.1 sandwich_2.3-4  
[12] tools_3.2.2  

编辑:(2015-02-22)
由于这引起了一些兴趣，我将尝试进一步澄清。我试图拟合一个“固定效应”模型(又名“内部”或“最小二乘虚拟变量”，正如 plm package vignette 在第 3 页，顶部段落中所称的那样)——相同的斜率，不同的截距。

这与在为 time 添加虚拟变量后运行普通 OLS 回归相同。和 id .使用下面的代码，我可以复制 plm 中的拟合值使用 base 封装 lm() .对于假人，很明显 id 和 time 的第一个元素是要比较的组。我还做不到的是如何使用 plm的设施。包做同样的事情我可以使用 lm() 轻松完成.

# fit the same with lm() and match the fitted values to those from plm()
lmF <- lm(data = DT, formula = y ~ x1 + x2 + factor(time) + factor(id))
time.lm <- coef(lmF)[grep(x = names(coef(lmF)), pattern = "time", fixed = TRUE)]
time.lm <- c(0, unname(time.lm)) # no need for names, the position index corresponds to time

id.lm <- coef(lmF)[grep(x = names(coef(lmF)), pattern = "id", fixed = TRUE)]
id.lm <- c(0, unname(id.lm))
names(id.lm) <- c("a","b","c","d") # set names so that individual values can be looked up below when generating the fit

DT[, by=list(id, time), fitted.lm := coef(lmF)[["(Intercept)"]]  +  coef(lmF)[["x1"]] * x1  +  coef(lmF)[["x2"]] * x2  +  time.lm[[time]]  +  id.lm[[id]]]
all.equal(DT$fitted.plm, DT$fitted.lm)

希望这对其他可能感兴趣的人有用。问题可能是关于如何 plm和 fixef处理我故意创建的缺失值。我试着玩 type= fixef 的参数但没有效果。

Answer 1

编辑:适应双向不平衡模型，需要plm版本>= 2.4-0
这是你想要的吗？通过 fixef 提取固定效应.以下是非平衡双向模型上的 Grunfeld 数据示例(对平衡双向模型的工作方式相同):

gtw_u <- plm(inv ~ value + capital, data = Grunfeld[-200, ], effect = "twoways")
yhat <- as.numeric(gtw_u$model[ , 1] - gtw_u$residuals) # reference
pred_beta <- as.numeric(tcrossprod(coef(gtw_u), as.matrix(gtw_u$model[ , -1])))
pred_effs <- as.numeric(fixef(gtw_u, "twoways")) # sum of ind and time effects

all.equal(pred_effs + pred_beta, yhat) # TRUE -> matches fitted values (yhat)

如果您想在其组件中拆分个体效应和时间效应的总和(由 effect = "twoways" 给出)，您需要选择一个引用，并且自然会想到两个，它们都在下面给出:

# Splits of summed up individual and time effects:
# use one "level" and one "dfirst"
ii <- index(gtw_u)[[1L]]; it <- index(gtw_u)[[2L]]
eff_id_dfirst <- c(0, as.numeric(fixef(gtw_u, "individual", "dfirst")))[ii]
eff_ti_dfirst <- c(0, as.numeric(fixef(gtw_u, "time",       "dfirst")))[it]
eff_id_level <- as.numeric(fixef(gtw_u, "individual"))[ii]
eff_ti_level <- as.numeric(fixef(gtw_u, "time"))[it]

all.equal(pred_effs, eff_id_level  + eff_ti_dfirst) # TRUE
all.equal(pred_effs, eff_id_dfirst + eff_ti_level)  # TRUE

(这是基于 fixef 的手册页， ?fixef 。那里还显示了如何处理(平衡和非平衡)单向模型)。