Uwaga: próbuję uzyskać kod do pracyobie czas i indywidualne efekty stałe orazniezrównoważony zestaw danych. Poniższy przykładowy kod działa ze zbilansowanym zestawem danych.

Zobacz także edycję poniżej

Próbuję ręcznie obliczyć dopasowane wartości modelu efektów stałych (z efektami indywidualnymi i czasowymi) za pomocąplm pakiet. Jest to bardziej ćwiczenie potwierdzające, że rozumiem mechanikę modelu i pakietu, wiem, że mogę sam dopasować wartości zplm obiekt, z dwóch powiązanych pytań (tutaj itutaj).

Odplm winieta (str. 2), podstawowym modelem jest:

y_it =alfa + beta_przezroczony *x_it + (mu_i +lambda_t +epsilon_to)

gdzie mu_i jest indywidualnym składnikiem terminu błędu (a.k.a. „efekt indywidualny”), a lambda_t to „efekt czasu”.

Ustalone efekty można wyodrębnić za pomocąfixef() i pomyślałem, że mogę je wykorzystać (razem z niezależnymi zmiennymi) do obliczenia dopasowanych wartości dla modelu, używając (z dwiema niezależnymi zmiennymi) w ten sposób:

dopasowanie_it =alfa + beta_1 *x1 + beta_2 *x2 + mu_i +lambda_t

To jest miejsce, w którym zawodzę - wartości, które otrzymuję, nie znajdują się w pobliżu dopasowanych wartości (które otrzymuję jako różnicę między wartościami rzeczywistymi a resztami w obiekcie modelu). Po pierwsze, nie widzęalpha gdziekolwiek. Próbowałem grać ze stałymi efektami pokazywanymi jako różnice od pierwszego, od średniego itd., Bez powodzenia.

Czego mi brakuje Może to być niezrozumienie modelu lub błąd w kodzie, boję się ... Z góry dziękuję.

PS: Jedno z powiązanych pytań to sugerujepmodel.response() powinien być związany z moją sprawą (a powodem nie maplm.fit funkcji), ale jej strona pomocy nie pomaga mi zrozumieć, co ta funkcja faktycznie robi, i nie mogę znaleźć żadnych przykładów, jak zinterpretować wynik, który ona wywołuje.

Dzięki!

Przykładowy kod tego, co zrobiłem:

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)

Moja sesja jest następująca:

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  

Edytuj: (2015-02-22) Ponieważ wzbudziło to pewne zainteresowanie, postaram się wyjaśnić dalej. Próbowałem dopasować model „stałych efektów” (a.k.a. „w obrębie” lub „najmniejszych kwadratów atrapowych zmiennych”, ponieważplm winieta wywołuje go na str. 3, górny akapit) - takie same nachylenie (-a), różne przechwyty.

To samo, co po uruchomieniu zwykłej regresji OLS po dodaniu manekinówtime iid. Używając poniższego kodu, mogę zduplikować dopasowane wartości zplm pakiet za pomocą bazylm(). W przypadku manekinów jest jasne, że pierwsze elementy zarówno id, jak i time są grupą do porównania. To, czego wciąż nie mogę zrobić, to korzystać z udogodnieńplm pakiet, aby zrobić to samo, mogę łatwo wykonać używająclm().

# 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)

Mam nadzieję, że jest to przydatne dla innych, którzy mogą być zainteresowani. Problem może dotyczyć tego, jakplm ifixef radzić sobie z brakującą wartością, którą celowo stworzyłem. Próbowałem grać ztype= parametrfixef ale bez skutku.