R:使用plm和pglm绘制面板模型预测

Question

我使用带有plm的线性面板模型创建了两个回归模型,并使用带有pglm包的泊松创建了一个广义面板模型.

library(plm); library(pglm)
data(Unions)  # from pglm-package
punions <- pdata.frame(Unions, c("id", "year"))

fit1 <- plm(wage ~ exper + rural + married, data=punions, model="random")
fit2 <- pglm(wage ~ exper + rural + married, data=punions, model="random", family="poisson")

我现在想通过绘制一组散点图中的拟合值来图形化地比较两个拟合.最好沿着这些行使用ggplot2:

library(ggplot2)
ggplot(punions, aes(x=exper, y=wage)) +
    geom_point() +
    facet_wrap(rural ~ married)

我考虑过简单地使用ggplot2 stat_smooth(),但(也许不足为奇)它似乎不能识别我数据的面板格式.手动提取预测值predict似乎也不适用于pglm模型.

如何在此图中叠加我的两个面板模型的预测值？

Answer 1

与@mtoto类似,我也不熟悉library(plm)或者library(gplm).但是预测方法plm是可用的,它只是没有输出.pglm没有预测方法.

R> methods(class= "plm")
[1] ercomp          fixef           has.intercept   model.matrix    pFtest          plmtest         plot            pmodel.response
 [9] pooltest        predict         residuals       summary         vcovBK          vcovHC          vcovSCC        
R> methods(class= "pglm")
no methods found

值得注意的是,我不明白你为什么使用泊松模型来计算工资数据.它显然不是泊松分布,因为它采用非整数值(下图).如果你愿意的话,你可以尝试使用负二项式,但我不确定它是否具有随机效果.但你可以使用MASS::glm.nb例如.

> quantile(Unions$wage, seq(0,1,.1))
         0%         10%         20%         30%         40%         50%         60%         70%         80%         90%        100% 
 0.02790139  2.87570334  3.54965422  4.14864865  4.71605855  5.31824370  6.01422463  6.87414349  7.88514525  9.59904809 57.50431282

解决方案1:使用 plm

punions$p <- plm:::predict.plm(fit1, punions)
# From examining the source code, predict.plm does not incorporate 
# the random effects, so you do not get appropriate predictions. 
# You just get the FE predictions.

ggplot(punions, aes(x=exper, y=p)) +
  geom_point() +
  facet_wrap(rural ~ married) 

解决方案2 - lme4

或者,您可以从lme4包中获得类似的拟合,它具有定义的预测方法:

library(lme4)
Unions$id <- factor(Unions$id)
fit3 <- lmer(wage ~ exper + rural + married + (1|id), data= Unions)
# not run:
fit4 <- glmer(wage ~ exper + rural + married + (1|id), data= Unions, family= poisson(link= "log"))

R> fit1$coefficients
(Intercept)       exper    ruralyes  marriedyes 
  3.7467469   0.3088949  -0.2442846   0.4781113 
R>  fixef(fit3)
(Intercept)       exper    ruralyes  marriedyes 
  3.7150302   0.3134898  -0.1950361   0.4592975 

我没有运行泊松模型,因为它显然是错误指定的.您可以进行某种变量处理以处理它或者可能是负二项式.无论如何,让我们完成这个例子:

# this has RE for individuals, so you do see dispersion based on the RE
Unions$p <- predict(fit3, Unions)
ggplot(Unions, aes(x=exper, y=p)) +
    geom_point() +
    facet_wrap(rural ~ married)

Answer 2

我不熟悉 pglm 包，但似乎没有类似的函数predict()可以从数据帧生成未来值的向量。

在所有其他情况下（应该都是tbh），您可以轻松地在 ggplot 中绘制此图，甚至使用构面换行。您只需将预测作为新列添加到数据框中：

punions$pred1 <- predict(fit1,punions,class="lm")

然后将其添加为附加项geom_line()：

    ggplot() + geom_point(data=punions, aes(x=exper, y=wage)) +
    geom_line(data=punions,aes(x=exper, y= pred1), color = "red") +
    facet_wrap(rural ~ married)