Search code examples
rstargazer

how to output a gee model using stargazer (not working for me)


Built some models that I want to present using stargazer but I observe the following error with my code:

Error in if (object.name$family$family == "gaussian") { : argument is of length zero

My code is:

library(drgee)
library(stargazer)

sat_mod01 <- gee(sat~GAD*Ratg + whtb + GRatio + Control_sec + Op_Ratio + Time_Cat + Yr, family = guassian, data = cd_df,
                 clusterid = "id", corstr = "AR-1")
mat_mod01 <- gee(mat~GAD*Ratg + whtb + GRatio + Control_sec + Op_Ratio + Time_Cat + Yr, family = guassian, data = cd_df,
                 clusterid = "id", corstr = "AR-1")
lat_mod01 <- gee(lat~GAD*Ratg + whtb + GRatio + Control_sec + Op_Ratio + Time_Cat + Yr, family = guassian, data = cd_df,
                 clusterid = "id", corstr = "AR-1")
ag_att_mod01 <- gee(ag_att~GAD*Ratg + whtb + GRatio + Control_sec + Op_Ratio + Time_Cat + Yr, family = guassian, 
            data = cd_df, clusterid = "id", corstr = "AR-1")
stargazer(sat_mod01, mat_mod01, lat_mod01, ag_att_mod01, title = "Regression Reults", align = T)

A minimum reproducible dataset is included below (data is sensitive but have reproduced variables, the code yields the same error.

structure(list(id = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), GAD = structure(c(2L, 2L, 1L, 
2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 2L), .Label = c("0", 
"1"), class = "factor"), whtb = structure(c(2L, 2L, 2L, 1L, 1L, 
2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 
2L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 
2L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 2L), .Label = c("0", 
"1"), class = "factor"), Yr = structure(c(4L, 5L, 7L, 7L, 7L, 
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L), .Label = c("2015", 
"2016", "2017", "2018", "2019", "2020", "2021"), class = "factor"), 
    GRatio = c(0.0523917995444191, 0.0523917995444191, 0.0523917995444191, 
    0.0523917995444191, 0.0523917995444191, 0.0523917995444191, 
    0.0523917995444191, 0.0523917995444191, 0.0523917995444191, 
    0.0523917995444191, 0.0523917995444191, 0.0523917995444191, 
    0.0523917995444191, 0.0411764705882353, 0.0411764705882353, 
    0.0411764705882353, 0.0373467112597547, 0.0373467112597547, 
    0.0373467112597547, 0.0373467112597547, 0.0373467112597547, 
    0.0373467112597547, 0.0373467112597547, 0.0373467112597547, 
    0.0373467112597547, 0.0373467112597547, 0.0373467112597547, 
    0.0373467112597547, 0.0373467112597547, 0.0373467112597547, 
    0.0373467112597547, 0.0373467112597547, 0.0373467112597547, 
    0.0373467112597547, 0.0373467112597547, 0.0373467112597547, 
    0.0373467112597547, 0.0373467112597547, 0.0373467112597547, 
    0.0373467112597547, 0.0373467112597547, 0.0373467112597547, 
    0.0373467112597547, 0.0373467112597547, 0.0373467112597547, 
    0.0373467112597547, 0.0373467112597547, 0.0373467112597547, 
    0.0373467112597547, 0.0373467112597547), Time_Cat = structure(c(1L, 
    3L, 1L, 2L, 1L, 2L, 3L, 1L, 2L, 2L, 1L, 2L, 2L, 3L, 1L, 1L, 
    1L, 3L, 2L, 2L, 3L, 3L, 3L, 2L, 3L, 3L, 3L, 2L, 3L, 1L, 1L, 
    3L, 3L, 2L, 1L, 3L, 3L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
    3L, 3L, 3L, 3L), .Label = c("Afternoon", "Evening", "Morning"
    ), class = "factor"), Control_sec = c(600, 180, 600, 600, 
    600, 1800, 1800, 1800, 1800, 900, 900, 900, 900, 600, 600, 
    60, 300, 600, 300, 300, 600, 600, 600, 300, 300, 300, 300, 
    300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 
    300, 300, 180, 180, 300, 300, 300, 300, 300, 300, 300), Ratg = c(722, 
    506, 529, 623, 856, 618, 585, 578, 570, 553, 539, 533, 541, 
    267, 342, 462, 597, 594, 565, 558, 580, 676, 729, 744, 766, 
    854, 836, 859, 767, 803, 776, 762, 742, 821, 820, 822, 811, 
    803, 808, 783, 620, 729, 753, 732, 730, 612, 611, 653, 766, 
    759), sat = c(2.35, 1.82, 2.88, 1.82, 2.53, 1.82, 2.35, 2.35, 
    2.35, 2.53, 1.82, 1.82, 1.82, 1.82, 1.82, 2.71, 2, 2.53, 
    2.71, 2.7, 2, 2.35, 1.82, 1.82, 2.35, 1.64, 2.53, 2.53, 1.82, 
    0.76, 2.53, 2.53, 1.82, 2.7, 2.35, 1.82, 1.64, 2, 2.35, 2.17, 
    1.64, 1.64, 2, 1.64, 1.82, 2.35, 2.35, 1.82, 2, 1.47), mat = c(0.83, 
    0.87, 2.88, 1.3, 1.35, 0.76, 0.94, 2.53, NA, 1.23, 1.05, 
    2.7, 0.76, 2.46, 2.7, 0.58, 1.82, 3.24, 0.87, 3.06, 1.64, 
    1.11, 2.64, 2.99, 1.65, 0.94, 1.82, 1.64, 1.35, 1.23, 1, 
    0.76, 2.11, 2.18, 0.65, 1.65, 2.35, 2.17, 1.82, 1.29, 0.94, 
    0.94, 2.6, 1.64, 0.76, 2.71, 2.17, 2, 2.06, 2.35), lat = c(0.83, 
    0, 0, 0, 0.36, 0.87, 0.76, 0.29, 2.59, 4.24, 3.71, 0.82, 
    1.59, 2.12, 1.64, 0.58, 0.65, 3.24, 0.87, 2.35, 2.24, 1.29, 
    2.75, 2.17, 0, 0.29, 0.76, 0.87, 0.58, 0.65, 0.58, 1.88, 
    1.64, 0.47, 0.65, 2.71, 0.76, 2.17, 1.05, 1.29, 0.76, 0, 
    2.6, 0.87, 1.88, 2.99, 2.88, 1.05, 2.24, 0.58), ag_att = c(1.33666666666667, 
    0.896666666666667, 1.92, 1.04, 1.41333333333333, 1.15, 1.35, 
    1.72333333333333, NA, 2.66666666666667, 2.19333333333333, 
    1.78, 1.39, 2.13333333333333, 2.05333333333333, 1.29, 1.49, 
    3.00333333333333, 1.48333333333333, 2.70333333333333, 1.96, 
    1.58333333333333, 2.40333333333333, 2.32666666666667, 1.33333333333333, 
    0.956666666666667, 1.70333333333333, 1.68, 1.25, 0.88, 1.37, 
    1.72333333333333, 1.85666666666667, 1.78333333333333, 1.21666666666667, 
    2.06, 1.58333333333333, 2.11333333333333, 1.74, 1.58333333333333, 
    1.11333333333333, 0.86, 2.4, 1.38333333333333, 1.48666666666667, 
    2.68333333333333, 2.46666666666667, 1.62333333333333, 2.1, 
    1.46666666666667), Op_Ratio = c(928, 603, 507, 685, 1084, 
    459, 498, 592, 636, 656, 513, 555, 569, 138, 355, 479, 621, 
    1299, 660, 573, 644, 696, 840, 689, 830, 838, 867, 866, 775, 
    897, 740, 563, 609, 883, 800, 946, 757, 826, 832, 838, 625, 
    729, 711, 578, 671, 667, 573, 694, 769, 852)), row.names = c(NA, 
-50L), class = c("tbl_df", "tbl", "data.frame"))

Without stargazer, they all run fine and output the desired result. I have searched thoroughly, the documentation for stargazer says it supports gee but have not seen any implementation in gee. Please any ideas?


Solution

  • I was able to replicate your problem and this appears to be a bug in stargazer. Unfortunately, it is very hard to diagnose the actual problem in that package because it consists of a single function that is 7000 lines long! Also, I’m not sure we should expect a fix soon since stargazer has not been updated since 2018.

    In the meantime, you may want to consider alternative regression table packages. For example, it is pretty easy to summarize your models with the modelsummary package (disclaimer: I am the author, so biased).

    modelsummary supports over one hundred models out of the box, but unfortunately not the ones produced by the drgee package. However, it is very easy to add support for these models by defining two simple functions: tidy.gee and glance.gee. Click here for details on the general strategy.

    library(modelsummary)
    library(drgee)
    
    tidy.gee <- function(x, ...) {
      out <- data.frame(
        term = names(coef(x)),
        estimate = coef(x),
        std.error = sqrt(diag(vcov(x)))
      )
      return(out)
    }
    
    glance.gee <- function(x, ...) {
      out <- data.frame(
        nobs = x$gee.data$n.obs
      )
      return(out)
    }
    
    models <- list(
      "SAT" = gee(sat ~ GAD * Ratg + whtb + GRatio + Control_sec + Op_Ratio + Time_Cat + Yr,
                  family = guassian, data = cd_df, clusterid = "id", corstr = "AR-1"),
      "MAT" = gee(mat ~ GAD * Ratg + whtb + GRatio + Control_sec + Op_Ratio + Time_Cat + Yr,
                  family = guassian, data = cd_df, clusterid = "id", corstr = "AR-1"),
      "LAT" = gee(lat ~ GAD * Ratg + whtb + GRatio + Control_sec + Op_Ratio + Time_Cat + Yr,
                  family = guassian, data = cd_df, clusterid = "id", corstr = "AR-1"),
      "AG ATT" = gee(ag_att ~ GAD * Ratg + whtb + GRatio + Control_sec + Op_Ratio + Time_Cat + Yr,
                     family = guassian, data = cd_df, clusterid = "id", corstr = "AR-1")
    )
    
    modelsummary(models)
    

    enter image description here