R: Plotting panel model predictions using plm & pglm

Question

I've created two regression models using a linear panel model with plm, and a generalized panel model using poisson with the pglm package.

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

I now want to compare the two fits graphically by plotting the fitted values in a set of scatterplots. Preferably along these lines using ggplot2:

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

I considered simply using ggplot2's stat_smooth(), but (perhaps unsurprisingly) it doesn't seem to recognize the panel format of my data. Manually extracting the predicted values with predict also does not seem to work for the pglm-model.

How do I overlay the predicted values for my two panel models in this plot?

Answer 1

Similar to @mtoto, I am also not familiar with either library(plm) or library(gplm). But the predict method for plm is available, it's just not exported. pglm does not have a predict method.

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

Of note, I do not understand why you are using a Poisson model for the wage data. It's clearly not a Poisson distribution since it takes non-integer values (below). You could try a negative binomial if you wish, though I'm not sure that's available with random effects. But you could use MASS::glm.nb for instance.

> 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

Solution 1: use 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) 

Solution 2 - lme4

Alternatively, you can get similar fits from the lme4 package, which does have a predict method defined:

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 

I haven't run the poisson models because it's clearly incorrectly specified. You could do some sort of variable transformation to handle it or perhaps a negative binomial. In any case, let's finish the example:

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