R Stargazer with pglm model - convert binominal pglm model in plm model

Question

I am using stargazer to create my plm summary tables.

library(plm)
library(pglm)
data("Unions", package = "pglm")
anb1 <- plm(wage ~ union + exper + rural, Unions, model = "random", method = "bfgs")
stargazer(anb1)

Unfortunately stargazer does not support pglm models. I am looking for a solution on how to plot the results of a pglm model with binary dependent variable, as the following stargazer call does not work with pglm models.

anb2 <- pglm(union ~ wage + exper + rural, Unions, family = "binomial",
            model = "random", method = "bfgs")
stargazer(anb2)

Any alternative to just extract each summary item and than format it respectively? The class of the outcome is:

[1] "maxLik" "maxim" "list"

I was surprised to see attempts to "plot" a model using `stargazer`. I thought it only formatted model summaries. Can you give details about what sort of "plot" you do desire? — IRTFM, Feb 10 '17 at 15:12
Hi 42, sorry for the confusion. I just meant to create nice looking output tables. — Nils, Feb 10 '17 at 15:46
Since `pglm` isn't supported by `stargazer`, and `stargazer` is not extensible, check out [`texreg`](https://cran.r-project.org/web/packages/texreg/vignettes/texreg.pdf) which is extensible for models not included in the package. — paqmo, Feb 13 '17 at 15:46

score 4 · Accepted Answer · answered May 19 '17 at 10:33

Here's a simple extract function to make texreg work with pglm:

extract.pglm <- function (model, include.nobs = TRUE, include.loglik = TRUE, ...) {
   s <- summary(model, ...)
   coefficient.names <- rownames(s$estimate)
   coefficients <- s$estimate[, 1]
   standard.errors <- s$estimate[, 2]
   significance <- s$estimate[, 4]
   loglik.value <- s$loglik
   n <- nrow(model$model)
   gof <- numeric()
   gof.names <- character()
   gof.decimal <- logical()
   if (include.loglik == TRUE) {
      gof <- c(gof, loglik.value)
      gof.names <- c(gof.names, "Log-Likelihood")
      gof.decimal <- c(gof.decimal, TRUE)
   }
   if (include.nobs == TRUE) {
      gof <- c(gof, n)
      gof.names <- c(gof.names, "Num. obs.")
      gof.decimal <- c(gof.decimal, FALSE)
   }
   tr <- createTexreg(coef.names = coefficient.names, coef = coefficients, 
                 se = standard.errors, pvalues = significance, gof.names = gof.names, 
                 gof = gof, gof.decimal = gof.decimal)
   return(tr)
}

In order for this code to work, you should also register the function so that it handles the pglm maxLik objects by default when extract is called:

setMethod("extract", signature = className("maxLik", "maxLik"), 
      definition = extract.pglm)

After that, you can use texreg with pglm just like with plm or other models supported by texreg.

score 4 · Answer 2 · answered Jun 25 '19 at 17:55

Here's a much simpler solution. Stargazer still doesn't support pglm as of 6/25/2019 but coeftest does so just pass the model to stargazer through coeftest.

(Also notice the change to the name of the data object in pglm since @giamcomo)

library(plm)
library(pglm)
library(lmtest)
library(stargazer)
data("UnionWage", package = "pglm")

anb2 <- pglm(union ~ wage + exper + rural, UnionWage, family = "binomial",
             model = "random", method = "bfgs")

stargazer(anb2)

summary(anb2)

stargazer(coeftest(anb2), type="text")

Here are the output

> stargazer(anb2)

% Error: Unrecognized object type.
> 
> summary(anb2)
--------------------------------------------
Maximum Likelihood estimation
BFGS maximization, 35 iterations
Return code 0: successful convergence 
Log-Likelihood: -1655.034 
5  free parameters
Estimates:
            Estimate Std. error t value  Pr(> t)    
(Intercept) -3.43651    0.29175 -11.779  < 2e-16 ***
wage         0.82896    0.15014   5.521 3.37e-08 ***
exper       -0.06590    0.02318  -2.843  0.00447 ** 
ruralyes     0.07558    0.24866   0.304  0.76116    
sigma        4.26050    0.22606  18.847  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
--------------------------------------------
> 
> stargazer(coeftest(anb2), type="text")

====================================
             Dependent variable:    
         ---------------------------

------------------------------------
wage               0.83***          
                   (0.15)           
exper             -0.07***          
                   (0.02)           
ruralyes            0.08            
                   (0.25)           
sigma              4.26***          
                   (0.23)           
Constant          -3.44***          
                   (0.29)           
====================================
====================================
Note:    *p<0.1; **p<0.05; ***p<0.01
>

This is clever, is there any way to retain R and N, etc? Or do these need to be added manually? — wolfsatthedoor, Aug 14 '19 at 20:04

score 1 · Answer 3 · answered May 01 '17 at 17:13

One possible solution is to do the following.

anb2 <- pglm(union ~ wage + exper + rural, Unions, family = "binomial",
         model = "random", method = "bfgs")
model = summary(anb2)

Load or install the following libraries

 library(dplyr)
 library(xtable)
 library('gtools')

Create a vector with the name of the co-variates

var = c('Intercept', 'wage', 'exper', 'ruralyes', 'sigma')

Then do

 model_summary = model$estimate %>% as.data.frame() %>% 
 mutate(term = var, Estimate = round(Estimate, 2), SE = round(`Std. error`, 2), p.value = stars.pval(`Pr(> t)`)) %>% 
 select(term, Estimate, SE, p.value)

 > model_summary
        term Estimate   SE p.value
 1 Intercept    -2.86 0.23     ***
 2      wage     0.12 0.02     ***
 3     exper    -0.06 0.02       *
 4  ruralyes     0.09 0.25        
 5     sigma     4.30 0.23     ***

Then you can use xtable on the data.frame

 library(xtable)
 xtable(model_summary)

score 1 · Answer 4 · answered May 01 '17 at 21:07

Another possibile (but not completely satisfactory) solution.

library(plm)
library(pglm)
library(stargazer)
data("Unions", package = "pglm")

anb2 <- pglm(union ~ wage + exper + rural, Unions, family = "binomial",
            model = "random", method = "bfgs")

# A "fake" model
anb0 <- plm(union ~ wage + exper + rural, Unions, family = "binomial",
            model = "random", method = "bfgs")

tstats <- summary(anb2)$estimate[,3][-5]
pvs <- summary(anb2)$estimate[,4][-5]
SEs <- summary(anb2)$estimate[,2][-5]
coefs <- summary(anb2)$estimate[,1][-5]

stargazer(anb0, type="text", coef=list(coefs), se=list(SEs),
 p = list(pvs), omit.stat="all")

Here is the table generated by stargazer:

====================================
             Dependent variable:    
         ---------------------------
                    union           
------------------------------------
wage              0.122***          
                   (0.024)          

exper             -0.058**          
                   (0.023)          

ruralyes            0.092           
                   (0.249)          

Constant          -2.857***         
                   (0.235)          

====================================
====================================
Note:    *p<0.1; **p<0.05; ***p<0.01

CZuC · Answer 5 · 2018-10-22T17:03:38.507

0

the new version of texreg (1.36.24) is available on GitHub (soon on CRAN) and pglm class has been added.

edited Oct 22 '18 at 17:03

answered Oct 18 '18 at 11:38

CZuC

69
6

1

I disagree with the reviewer. I think this does provide an answer. – Håken Lid Oct 18 '18 at 18:49
@HåkenLid I also agree. If reviewers believe that this answer is wrong or not useful, the correct action is to downvote. – Robert Columbia Oct 19 '18 at 03:40

R Stargazer with pglm model - convert binominal pglm model in plm model

5 Answers5