R/plm 按索引提取残差

Question

我有一个使用以下方法创建的 plm 对象:

require(plm)
plm1 <- plm(Sepal.Length ~ Petal.Length + Petal.Width, data = iris, index = "Species")

我正在尝试提取残差以手动计算 Species 的 r 平方，似乎无法将 pseries 对象操作为像矩阵或 data.frame 这样可用的东西。

> data.frame(resid(plm1))
Error in as.data.frame.default(x[[i]], optional = TRUE) : 
  cannot coerce class '"pseries"' into a data.frame

如果我有类似的东西会很好:

> df1 <- data.frame(time = rep(1:10,15), Species = iris$Species, resid1 = runif(150))
> head(df1)
  time Species    resid1
1    1  setosa 0.7038776
2    2  setosa 0.2164597
3    3  setosa 0.1988884
4    4  setosa 0.9311872
5    5  setosa 0.7087211
6    6  setosa 0.9914357

我可以使用 ddply 或聚合来找到每个物种的 rsquared。

有什么建议？

Answer 1

这是一个老问题，但我想指出一些容易遗漏并且可能导致严重错误的内容。 previous answer by dickoa是正确的，但我想我会澄清为什么需要这样的解决方法，因为它可能并不明显。

阅读时another thread我学到了以下内容:如前所述 here , plm 不一定按照提供给函数的顺序保存数据。这意味着只需使用 residuals()如果您不小心，在 plm 对象上使用函数然后将其加入您的数据可能会导致错误的残差分组到错误的数据行!作为示例，请考虑以下内容:

require(plm)
data("Gasoline") # The Gasoline dataset from the plm package

plm1 <- plm(lgaspcar ~ lincomep + lrpmg + lcarpcap, data=Gasoline, method = "within", index = c("country", "year"))

coef(plm1)
  lincomep      lrpmg   lcarpcap 
 0.6622497 -0.3217025 -0.6404829 

head(residuals(plm1))
          1           2           3           4           5           6 
-0.18814207 -0.19642727 -0.14874420 -0.12476346 -0.12114060 -0.08684045 

注意我们得到的残差。现在让我们更改数据集的排序顺序。这不应该改变分析中的任何内容。

set.seed(1234)
Gasoline2 <- Gasoline[order(runif(nrow(Gasoline))), ] # We just change the order of the rows.

plm2 <- plm(lgaspcar ~ lincomep + lrpmg + lcarpcap, data=Gasoline2, method = "within", index = c("country", "year"))

coef(plm2)
  lincomep      lrpmg   lcarpcap 
 0.6622497 -0.3217025 -0.6404829 

head(residuals(plm2))
        258           7          64          73         268         186 
-0.18814207 -0.19642727 -0.14874420 -0.12476346 -0.12114060 -0.08684045 

乍一看，这似乎很好；估计系数与以前相同。但是，请注意残差的显示顺序与我们移动行之前的顺序相同。唯一改变的是与残差关联的名称现在反射(reflect)了它们在数据中的新位置。因此，在数据中的第 1 行进行后重新排序的观察结果是在第 258 行上进行了预重新排序。

Gasoline2[1, ]
    country year lgaspcar lincomep     lrpmg  lcarpcap
258  SWEDEN 1970 3.989372 -7.73261 -2.733592 -8.164506

Gasoline[258, ]
    country year lgaspcar lincomep     lrpmg  lcarpcap
258  SWEDEN 1970 3.989372 -7.73261 -2.733592 -8.164506

这意味着，如果我们使用 Gasoline2 作为我们正在使用的数据集，那么使用类似 cbind() 的函数在 Gasoline2和 residuals(plm2)将导致错误的残差与观测值相关联。

head(cbind(Gasoline, residuals(plm1)))
  country year lgaspcar  lincomep      lrpmg  lcarpcap residuals(plm1)
1 AUSTRIA 1960 4.173244 -6.474277 -0.3345476 -9.766840     -0.18814207
2 AUSTRIA 1961 4.100989 -6.426006 -0.3513276 -9.608622     -0.19642727
3 AUSTRIA 1962 4.073177 -6.407308 -0.3795177 -9.457257     -0.14874420
4 AUSTRIA 1963 4.059509 -6.370679 -0.4142514 -9.343155     -0.12476346
5 AUSTRIA 1964 4.037689 -6.322247 -0.4453354 -9.237739     -0.12114060
6 AUSTRIA 1965 4.033983 -6.294668 -0.4970607 -9.123903     -0.08684045

head(cbind(Gasoline2, residuals(plm2)))
     country year lgaspcar  lincomep      lrpmg  lcarpcap residuals(plm2)
258   SWEDEN 1970 3.989372 -7.732610 -2.7335921 -8.164506     -0.18814207
7    AUSTRIA 1966 4.047537 -6.252545 -0.4668377 -9.019822     -0.19642727
64   DENMARK 1966 4.233643 -5.851866 -0.3961885 -8.681541     -0.14874420
73   DENMARK 1975 4.033015 -5.612967 -0.3939543 -8.274632     -0.12476346
268 SWITZERL 1961 4.441330 -6.111640 -0.8655847 -9.158229     -0.12114060
186    JAPAN 1974 4.007964 -5.852553 -0.1909064 -8.846520     -0.08684045

正如我们在上面看到的，残差被分配到 Gasoline2 示例中的错误行。

发生什么了？好吧，正如前面提到的， plm不保留观察的顺序。使用 attr()函数 dickoa 在之前的回答中指出，我们可以看到 plm按国家和年份重新组织数据。

head( attr(residuals(plm2), "index") )
  country year
1 AUSTRIA 1960
2 AUSTRIA 1961
3 AUSTRIA 1962
4 AUSTRIA 1963
5 AUSTRIA 1964
6 AUSTRIA 1965

这就是原始汽油数据的结构方式，这就是残差以相同顺序显示的原因。

因此，我们可以使用 attr(residuals(plm2), "index") 的事实。为我们提供残差及其相应的国家和年份指标，以便将残差添加到原始数据中。正如所指出的 here , plyr包对此非常有帮助。

require(plyr)
resids2 <- data.frame(residual = residuals(plm2), attr(residuals(plm2), "index"))
Gasoline2$year <- factor(Gasoline2$year) # Needed since resids2$year is a factor, and Gasoline2$years was an integer. plyr does not accept them to be of different types.
Gasoline2 <- join(Gasoline2, resids2, by = c("country", "year"))

head(Gasoline2)
   country year lgaspcar  lincomep      lrpmg  lcarpcap    residual
1   SWEDEN 1970 3.989372 -7.732610 -2.7335921 -8.164506 -0.02468148
2  AUSTRIA 1966 4.047537 -6.252545 -0.4668377 -9.019822 -0.02479759
3  DENMARK 1966 4.233643 -5.851866 -0.3961885 -8.681541  0.03175032
4  DENMARK 1975 4.033015 -5.612967 -0.3939543 -8.274632 -0.06575219
5 SWITZERL 1961 4.441330 -6.111640 -0.8655847 -9.158229 -0.05789130
6    JAPAN 1974 4.007964 -5.852553 -0.1909064 -8.846520 -0.21957156

这给了我们正确的结果。