13

注意:時間と個々の固定効果、および不均衡なデータセットの両方でコードを動作させようとしています。以下のサンプル コードは、バランスの取れたデータセットで動作します。

以下の編集もご覧ください。

plmパッケージを使用して、固定効果モデル (個人効果と時間効果の両方) の適合値を手動で計算しようとしています。plmこれは、モデルとパッケージの仕組みを理解していることを確認するための演習です。2 つの関連する質問 (こことここ)から、オブジェクトから適合値自体を取得できることを知っています。

plmビネット (p.2) から、基になるモデルは次のとおりです。

y _it = alpha + beta _transposed * x _it + ( mu _i + lambda _t + epsilon _it)

ここで、mu_i は誤差項の個別の要素 (別名「個別効果」) であり、lambda_t は「時間効果」です。

固定効果は次のように使用して抽出できます。fixef()これらを (独立変数と一緒に) 使用して、(2 つの独立変数を使用して) モデルの適合値を計算できると思いました。

フィット_it =アルファ+ベータ_1 * x1 +ベータ_2 * x2 +ミュー_i +ラムダ_t

これは私が失敗するところです - 私が得た値は適合値にどこにも近くありません(モデルオブジェクトの実際の値と残差の差として得られます)。一つには、私はどこにも見えませんalpha。固定効果が最初からの差、平均値からの差などとして表示されるようにしてみましたが、成功しませんでした。

私は何が欠けていますか？モデルの誤解、またはコードのエラーである可能性が高いと思います...よろしくお願いします。

pmodel.response()PS: 関連する質問の 1 つは、私の問題 (および機能がない理由) に関連するはずのヒントです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)

私のセッションは次のとおりです。

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パッケージのビネットがp.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

4

4 に答える 4

1

私の場合、 lm() ソリューションが機能しなかったため、これが役立つことがわかりました（plmパッケージと比較して異なる係数が得られました）

したがって、ここで plm パッケージの作成者による提案を適用するだけですhttp://r.789695.n4.nabble.com/fitted-from-plm-td3003924.html

だから私がしたことはただ適用することです

plm.object <- plm(y ~ lag(y, 1) + z +z2, data = mdt, model= "within", effect="twoways")
fitted <- as.numeric(plm.object$model[[1]] - plm.object$residuals)

as.numeric 関数が必要なのは、さらに操作するためにプラグインするためのベクトルとして使用する必要があるためです。また、モデルの右側に遅延従属変数がある場合、上記の as.numeric を使用したソリューションは、遅延のために欠損値の NET のベクトルを提供することも指摘したいと思います。私にとって、これはまさに私が必要としていたものです。

于 2015-04-20T16:08:51.660 に答える
0

編集: 双方向不平衡モデルに適合、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"その構成要素に分割する場合は、参照を選択する必要があり、以下に示す 2 つが自然に思い浮かびます。

# 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 の man ページに基づいています?fixef。そこには、(平衡および不平衡) 一方向モデルの処理方法も示されています)。

于 2015-11-12T14:28:06.040 に答える