Plm: используя fixef(), чтобы вручную вычислить установленные значения для модели с фиксированными эффектами twoways

Обратите внимание: я пытаюсь заставить код работать как с временными, так и с отдельными фиксированными эффектами и с несбалансированным набором данных. Пример кода ниже работает со сбалансированным набором данных.

См. также редактирование ниже, пожалуйста

Я пытаюсь вручную вычислить установленные значения модели фиксированных эффектов (с индивидуальными и временными эффектами) с помощью пакета plm. Это скорее упражнение, чтобы подтвердить, что я понимаю механику модели и пакета, я знаю, что могу получить установленные значения из объекта plm из двух связанных вопросов (here и здесь).

Из виньетки plm (стр. 2) базовая модель:

y _it = alpha + betastrong > _перемещенный * x _it + ( mu _i + lambda _t + epsilon _it)

где mu_i - индивидуальная составляющая члена ошибки (a.k.a. "индивидуальный эффект" ), а lambda_t - "эффект времени".

Фиксированные эффекты можно извлечь, используя fixef(), и я подумал, что могу использовать их (вместе с независимыми переменными) для вычисления установленных значений для модели, используя (с двумя независимыми переменными) следующим образом:

fit _it = alpha + betastrong > _1 * x1 + betastrong > _2 * x2 + mu _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)

Мой сеанс выглядит следующим образом:

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, верхний абзац) - тот же наклон (ы), различные перехваты.

Это то же самое, что запустить обычную регрессию OLS после добавления манекенов для time и id. Используя приведенный ниже код, я могу дублировать установленные значения из пакета plm с помощью базы 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, но без эффекта.

Ответы

Ответ 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 из отсутствующих значений из-за задержки. Для меня это именно то, что мне нужно.

Ответ 2

Это то, что вы хотели? Извлеките фиксированные эффекты на fixef и сопоставьте их с отдельным индексом. Вот пример данных Grunfeld:

data(Grunfeld, package = "plm")
fe <- plm(inv ~ value + capital, data=Grunfeld, model = "within")

temp <- merge(Grunfeld, data.frame(fixef_firm = names(fixef(fe)), fixef = as.numeric(fixef(fe))), all.x =T, by.x = c("firm"), by.y=c("fixef_firm"))
fitted_by_hand <- temp$fixef + fe$coefficients[1] * Grunfeld$value +  fe$coefficients[2] * Grunfeld$capital

fitted <- fe$model[ , 1] - fe$residuals

# just to remove attributs and specific classes 
fitted_by_hand <- as.numeric(fitted_by_hand)
fitted <- as.numeric(fitted)

all.equal(fitted, fitted_by_hand) # TRUE
cbind(fitted, fitted_by_hand) # see yourself

Ответ 3

Это работает для несбалансированных данных с effect="individual" и временными манерами y ~ x +factor(year):

fitted <- pmodel.response(plm.model)-residuals(plm.model)

Ответ 4

Я очень близко отношусь к предложению Helix123, чтобы вычесть within_intercept (он включается в каждый из двух фиксированных эффектов, поэтому вам нужно исправить это).

В моих ошибках реконструкции есть очень наводящий пример: индивидуальный a всегда выключен на -0.004858712 (за каждый период времени). Лица b, c, d всегда отключены на 0.002839703 за каждый период времени, за исключением периода 4 (где нет наблюдения за a), где они отключены на -0.010981192.

Любые идеи? Похоже, что отдельные фиксированные эффекты отбрасываются дисбалансом. Правильная балансировка работает правильно.

Полный код:

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)

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

summary(plmFEit)

DT[, resids := residuals(plmFEit)]

FEI <- data.table(as.matrix(fixef(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(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

DT[, fitted.calc := plmFEit$coefficients[[1]] * x1 + plmFEit$coefficients[[2]] * x2 +
     fei + fet - within_intercept(plmFEit)]

DT[, myresids := y - fitted.calc]
DT[, myerr := resids - myresids]