Объединить подогнанные значения plm в набор данных

Объединение plm подогнанных значений обратно в исходный набор данных требует некоторых промежуточных шагов — plm отбрасывает все строки с отсутствующими данными, и, насколько я могу судить, объект plm не содержит информации об индексе. Порядок данных не сохраняется — см. комментарий Джованни Милло, одного из авторов plm, в этой темы:

«... порядок ввода не всегда сохраняется: наблюдения всегда переупорядочиваются (индивидуально, по времени) внутренне, так что выходные данные, которые вы получаете, упорядочиваются соответствующим образом ...»

Шаги вкратце:

1. Получите подходящие значения из предполагаемого объекта plm. Это один вектор, но записи названы. Имена соответствуют позиции в указателе.
2. Получите индекс, используя функцию index(). Он может возвращать как индивидуальные, так и временные индексы. Обратите внимание, что индекс может содержать больше строк, чем подобранные значения, если строки были удалены из-за отсутствующих данных. (Также можно сгенерировать индекс непосредственно из исходных данных, но я не видел обещания, что исходный порядок данных сохраняется в том, что возвращает plm.)
3. Слияние с исходными данными, поиск значений идентификатора и времени из индекса.

Пример кода приведен ниже. Немного длинно, но я попытался прокомментировать. Код не оптимизирован, я хотел явно перечислить шаги. Кроме того, я использую data.tables, а не data.frames.

library(data.table); library(plm)

### Generate dummy data. This way we know the "true" coefficients
set.seed(100)
n <- 500 # Run with more data if you want to get closer to the "true" coefficients
DT <- data.table(CJ(id = c("a","b","c","d","e"), time = c(1:(n / 5))))
DT[, x1 := rnorm(n)]
DT[, x2 := rnorm(n)]
DT[, y  := x1 + 2 * x2 + rnorm(n) / 10]

setkey(DT, id, time)
# # Make it an unbalanced panel & put in some NAs
DT <- DT[!(id == "a" & time == 4)]
DT[.("a", 3), x2 := as.numeric(NA)]
DT[.("d", 2), x2 := as.numeric(NA)]

str(DT)

### Run the model -- both individual and time effects; "within" model
summary(PLM <- plm(data = DT, id = c("id", "time"), formula = y ~ x1 + x2, model = "within", effect = "twoways", na.action = "na.omit"))

### Merge the fitted values back into the data.table DT
# Note that PLM$model$y is shorter than the data, i.e. the row(s) with NA have been dropped
cat("\nRows omitted (due to NA): ", nrow(DT) - length(PLM$model$y))

# Since the objects returned by plm() do not contain the index, need to generate it from the data
# The object returned by plm(), i.e. PLM$model$y, has names that point to the place in the index
# Note: The index can also be done as INDEX <- DT[, j = .(id, time)], but use the longer way with index() in case plm does not preserve the order
INDEX <- data.table(index(x = pdata.frame(x = DT, index = c("id", "time")), which = NULL)) # which = NULL extracts both the individual and time indexes
INDEX[, id := as.character(id)]
INDEX[, time := as.integer(time)] # it is returned as a factor, convert back to integer to match the variable type in DT

# Generate the fitted values as the difference between the y values and the residuals
if (all(names(PLM$residuals) == names(PLM$model$y))) { # this should not be needed, but just in case...
    FIT <- data.table(
        index   = as.integer(names(PLM$model$y)), # this index corresponds to the position in the INDEX, from where we get the "id" and "time" below
        fit.plm = as.numeric(PLM$model$y) - as.numeric(PLM$residuals)
    )
}

FIT[, id   := INDEX[index]$id]
FIT[, time := INDEX[index]$time]
# Now FIT has both the id and time variables, can match it back into the original dataset (i.e. we have the missing data accounted for)
DT <- merge(x = DT, y = FIT[, j = .(id, time, fit.plm)], by = c("id", "time"), all = TRUE) # Need all = TRUE, or some data from DT will be dropped!

Остатки - это отклонение модели от значения на LHS формулы .... которую вы нам не показали. В пакете 'plm' есть функция fitted.panelmodel, но, похоже, ожидается, что будет значение fitted, которое функция plm не возвращает по умолчанию, и это не задокументировано, а также способ, который я вижу чтобы заставить его кашлять.

library(plm)
data("Produc", package = "plm")
zz <- plm(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp,
          data = Produc, index = c("state","year"))
summary(zz)  # the example on the plm page:
> str(fitted(zz))
 NULL
> names(zz$model)
[1] "log(gsp)"  "log(pcap)" "log(pc)"   "log(emp)"  "unemp"    
> Produc[ , c("Yvar", "Fitted")] <- cbind( zz$model[ ,"log(gsp)", drop=FALSE], zz$residuals)
> str(Produc)
'data.frame':   816 obs. of  12 variables:
 $ state : Factor w/ 48 levels "ALABAMA","ARIZONA",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ year  : int  1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 ...
 $ pcap  : num  15033 15502 15972 16406 16763 ...
 $ hwy   : num  7326 7526 7765 7908 8026 ...
 $ water : num  1656 1721 1765 1742 1735 ...
 $ util  : num  6051 6255 6442 6756 7002 ...
 $ pc    : num  35794 37300 38670 40084 42057 ...
 $ gsp   : int  28418 29375 31303 33430 33749 33604 35764 37463 39964 40979 ...
 $ emp   : num  1010 1022 1072 1136 1170 ...
 $ unemp : num  4.7 5.2 4.7 3.9 5.5 7.7 6.8 7.4 6.3 7.1 ...
 $ Yvar  :Classes 'pseries', 'pseries', 'integer'  atomic [1:816] 10.3 10.3 10.4 10.4 10.4 ...
  .. ..- attr(*, "index")='data.frame': 816 obs. of  2 variables:
  .. .. ..$ state: Factor w/ 48 levels "ALABAMA","ARIZONA",..: 1 1 1 1 1 1 1 1 1 1 ...
  .. .. ..$ year : Factor w/ 17 levels "1970","1971",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ Fitted: num  -0.04656 -0.03064 -0.01645 -0.00873 -0.02708 ...

У меня упрощенный метод. Основная проблема здесь двоякая:

1) pdata.frames сортирует ввод в алфавитном порядке по имени, а затем по году. Это можно решить, сначала отсортировав фрейм данных перед запуском plm.

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

library(plm)
set.seed(100)
n <- 10 # Run with more data if you want to get closer to the "true" coefficients
DT <- data.frame(id = c("a","c","b","d","e"), time = c(1:(n / 5)),x1 = rnorm(n),x2= rnorm(n),x3=rnorm(n))
DT$Y = DT$x2 + 2 * DT$x3 + rnorm(n) / 10 # make x1 a function of other variables
DT$x3[3]=NA  # add an NA to show this works with missing data 
DT  

# now can add drop.index = F, but note that DT is now sorted by order(id,time)
pdata.frame(DT,index=c('id','time'),drop.index = F)

# order DT to match pdata.frame that will be used for plm
DT=DT[order(DT$id,DT$time),]

# formulas
formulas =Y~x1+x2+x3 
formulas_dataframe = Y~x1+x2+x3 +id+time # add id and time for model.frame

# estimate
random <- plm(formulas, data=DT, index=c("id", "time"), model="random",na.action = 'na.omit')
summary(random) 

# merge prediction and and model.frame 
fitted = data.frame(fitted = random$model[[1]] - random$residuals)
model_data = cbind(as.data.frame(as.matrix(random$model)),fitted)  # this isn't really needed but shows that input and model.frame are same
model_data = cbind(model_data,na.omit(model.frame(formulas_dataframe,DT)))  
model_data

Я написал функцию (predict.out.plm), чтобы делать прогнозы из выборки после оценки моделей первых различий или фиксированных эффектов с помощью plm.

Далее функция добавляет предсказанные значения к индексам исходных данных. Это делается с помощью rownames, сохраненных в plm - attributes(plmobject)$index, и rownames в model.matrix.

для получения более подробной информации см. функцию, размещенную здесь:

https://stackoverflow.com/a/44185441/2409896

Для этого поста прошло некоторое время, но я считаю, что самый простой способ сделать это сейчас:

Fe.model.fitted <- cbind(FE.model$model, 
                         resid=FE.model$residuals, 
                         fitted=plm:::fitted_exp.plm(FE.model))

Функция fitted_exp.plm не экспортируется пакетом plm, но мы можем использовать ::: для ее извлечения.