Support side-by-side outputs from multi-models

[lmwang9527]

In academic papers, it is common to show estimation outputs from similar models side-by-side to facilitate comparison. stargazer supports this:

fit1 <- lm(mpg ~ qsec + am + wt + gear + factor(vs), data = mtcars)
fit2 <- lm(mpg ~ am + wt + gear + factor(vs), data = mtcars)
    
library(stargazer)
stargazer(fit1, fit2, 
          column.labels = c("Model1", "Model2"), 
          dep.var.caption = "", 
          dep.var.labels.include=F, 
          type="text", no.space=TRUE)

which produces:

=================================================================
                            Model1                 Model2        
                             (1)                    (2)          
-----------------------------------------------------------------
qsec                       1.211**                               
                           (0.473)                               
am                          3.088                  2.485         
                           (1.823)                (1.985)        
wt                        -3.918***              -3.786***       
                           (0.826)                (0.905)        
gear                        -0.149                 -0.862        
                           (1.066)                (1.130)        
factor(vs)1                 0.052                 3.708***       
                           (1.851)                (1.291)        
Constant                    10.345               32.817***       
                           (9.983)                (5.211)        
-----------------------------------------------------------------
Observations                  32                     32          
R2                          0.850                  0.812         
Adjusted R2                 0.821                  0.784         
Residual Std. Error    2.551 (df = 26)        2.800 (df = 27)    
F Statistic         29.415*** (df = 5; 26) 29.145*** (df = 4; 27)
=================================================================
Note:                                 *p<0.1; **p<0.05; ***p<0.01

With pixieddust, what I came up is something like this:

library(dplyr)
library(broom)
library(pixiedust)
library(tibble)
library(stringr)

tf1 <- broom::tidy(fit1) %>% select(-statistic) 
tf2 <- broom::tidy(fit2) %>% select(-statistic)

tf <- tf1 %>% 
  full_join(tf2, by="term", suffix=c("_Model1", "_Model2"))

custom_head <- (str_split(names(tf), "_", simplify=TRUE)) %>% 
  as_data_frame() %>% 
  select(V2, V1) %>% 
  mutate(V2=ifelse(V1=="std.error", V2, "")) %>% 
  t() %>% 
  as_data_frame()

glance_stats <- c("adj.r.squared", "df")
custom_foot <- bind_rows(broom::glance(fit1), broom::glance(fit2)) %>% 
  t() %>% 
  as.data.frame() %>% 
  rownames_to_column() %>%
  filter(rowname %in% glance_stats) %>% 
  transmute(term=rowname, X1="", X2=V1, X3="", X4="", X5=V2, X6="") %>% 
  as_data_frame()

dust(tf) %>% 
  redust(custom_head, part = "head") %>% 
  redust(custom_foot, part = "foot") %>% 
  sprinkle(rows = 1, border = "top") %>%
  sprinkle(cols = 5, border = "left", part="head") %>%
  sprinkle(cols = 5, border = "left", part="body") %>%
  sprinkle(cols = 5, border = "left", part="foot") %>%
  sprinkle_table(round=3, na_string="") %>% 
  sprinkle(rows = 1, border = "top", part = "foot") %>%
  sprinkle(rows = 2, border = "bottom", part = "foot") %>%
  sprinkle_print_method("html")

It works, but what I don't like is that there is too much messing around (and too little magic) and apparently not scalable (what if I want to add one more model). Can pixiedust add support of this (or is there a better way that I am not aware of)?

BTW, thanks for the great work creating the package!

[nutterb]

The current devel branch includes a gaze function. It doesn't actually apply dust, but it does format into the stargazer like format. I haven't worked out significance stars yet, but it has been on my mind.

devtools::install_github("nutterb/pixiedust", ref = "devel")

fit1 <- lm(mpg ~ qsec + am + wt + gear + factor(vs), data = mtcars)
fit2 <- lm(mpg ~ am + wt + gear + factor(vs), data = mtcars)

gaze(fit1, fit2)
gaze(with_qsec = fit1, 
     without_qsec = fit2)
gaze(fit1, fit2, include_glance = FALSE)
gaze(fit1, fit2, glance_vars = c("AIC", "BIC"))

Still to do

• [ ] Requirement
• [ ] Tests
• [ ] Review documentation
• [ ] Finalize defaults for glance_vars