Grouping coefficients in plot R

Viewed 30

I was looking to group coefficient from a single model in a dot and whisker plot:

Here is a sample:

    ## Simulate data
    set.seed(42)
    n <- 1e3
    
    d <- data.frame(
      # Covariates
      projects_0_300_minus2 = rnorm(n),
      projects_0_300_minus1 = rnorm(n),
      projects_0_300_0 = rnorm(n),
      projects_0_300_1 = rnorm(n),
      projects_0_300_2 = rnorm(n),
      projects_300_600_minus2 = rnorm(n),
      projects_300_600_minus1 = rnorm(n),
      projects_300_600_0 = rnorm(n),
      projects_300_600_1 = rnorm(n),
      projects_300_600_2 = rnorm(n),
      projects_600_800_minus2 = rnorm(n),
      projects_600_800_minus1 = rnorm(n),
      projects_600_800_0 = rnorm(n),
      projects_600_800_1 = rnorm(n),
      projects_600_800_2 = rnorm(n),
      # Individuals and firms
      id = factor(sample(20, n, replace=TRUE)),
      firm = factor(sample(13, n, replace=TRUE)),
      # Noise
      u = rnorm(n)
    )
    
    
    # Effects for individuals and firms
    id.eff <- rnorm(nlevels(d$id))
    firm.eff <- rnorm(nlevels(d$firm))
    
    # Left hand side
    d$y <- d$projects_0_300_minus2 + 0.5*d$projects_0_300_minus1 + d$projects_0_300_0 + 0.5*d$projects_0_300_1 +d$projects_0_300_2  +d$projects_300_600_minus2 +d$projects_300_600_minus1 + 0.5*d$projects_300_600_0 +d$projects_300_600_1 + 0.5*d$projects_300_600_2  + 0.5*d$projects_600_800_0 +d$projects_600_800_minus2 + 0.5*d$projects_600_800_minus1 +d$projects_600_800_0 + 0.5*d$projects_600_800_1 +d$projects_600_800_2+ id.eff[d$id] + firm.eff[d$firm] + d$u
    
## extract predictor names
    d_pred <- subset(d, select = -c(id,firm,u,y))
    d_pred_names <- names(d_pred)
    
## model
    Formula <- formula(paste("y ~ ", paste(d_pred_names, collapse=" + "),"| id + firm | 0 | 0"))
    results <- felm(Formula, data = d)
        
    #extract names for grouping of coefficient s
    projects_minus2 <- names(d[ , grepl( "00_minus2" , names( d ) ) ])
    projects_minus1 <- names(d[ , grepl( '00_minus1' , names( d ) ) ])
    projects_0 <- names(d[ , grepl( "00_0" , names( d ) ) ])
    projects_1 <- names(d[ , grepl( "00_1" , names( d ) ) ])
    projects_2 <- names(d[ , grepl( "00_2" , names( d ) ) ])

## plot individual BUT WANT TO GROUP THESE INTO ONE PLOT    
    Turnover_projects_minus2 <- jtools::plot_coefs(results, coefs = projects_minus2) + labs(title = 'Turnover_projects -2')
    Turnover_projects_minus1 <- jtools::plot_coefs(results, coefs = projects_minus1) + labs(title = 'Turnover_projects -1')
    Turnover_projects_0 <- jtools::plot_coefs(results, coefs = projects_0) + labs(title = 'Turnover_projects 0')
    Turnover_projects_1 <- jtools::plot_coefs(results, coefs = projects_1) + labs(title = 'Turnover_projects 1')
    Turnover_projects_2 <- jtools::plot_coefs(results, coefs = projects_2) + labs(title = 'Turnover_projects 2')
    
    
    ggarrange(Turnover_projects_minus2, Turnover_projects_minus1, Turnover_projects_0,Turnover_projects_1,Turnover_projects_2,
                    ncol = 2, nrow = 5)

Using this image as an example:

enter image description here

I would like the variables that include seperately: 0-300, 300-600, 600-800, which will be on the Y axis where the above image has 'wt' , 'cyl', 'disp'.

Then, for each one of these, they would grouped by 5 points that include 'minus2', 'minus1' , '0' , '1' , '2' like the above image where 'model1', 'model2', 'model3' are grouped together.

Something like this i think but couldnt get it work..: https://www.cgoodman.com/blog/archives/2019/12/03/custom-coefficient-graphs/

1 Answers

To my mind, it would be easier to obtain the estmates with confidence intervals directly from your model and plot these. Taking your felm output results, we can do:

library(tidyverse)

plot_df <- as.data.frame(summary(results)$coef) %>%
  rownames_to_column() %>%
  mutate(var1 = sub("^.*_(\\d+_\\d+)_.*$", "\\1", rowname),
         var1 = sub("_", " - ", var1),
         var2 = sub("^.*_(.*)$", "\\1", rowname),
         var2 = factor(var2, c("minus1", "minus2", 0:2)),
         upper = Estimate + 1.96 * `Std. Error`,
         lower = Estimate - 1.96 * `Std. Error`) %>%
  select(var1, var2, Estimate, lower, upper)

This gives us the two variables required for grouping, the estimates, and the confidence intervals:

plot_df
#>         var1   var2  Estimate     lower     upper
#> 1    0 - 300 minus2 0.9377860 0.8725008 1.0030712
#> 2    0 - 300 minus1 0.4913298 0.4253458 0.5573139
#> 3    0 - 300      0 1.0522792 0.9884674 1.1160910
#> 4    0 - 300      1 0.4887432 0.4232850 0.5542015
#> 5    0 - 300      2 0.9327972 0.8685721 0.9970223
#> 6  300 - 600 minus2 0.9883478 0.9257882 1.0509073
#> 7  300 - 600 minus1 0.9717786 0.9052546 1.0383026
#> 8  300 - 600      0 0.5193846 0.4546181 0.5841512
#> 9  300 - 600      1 1.0052123 0.9400112 1.0704134
#> 10 300 - 600      2 0.4651625 0.4003738 0.5299512
#> 11 600 - 800 minus2 1.0270672 0.9627909 1.0913436
#> 12 600 - 800 minus1 0.5258786 0.4608404 0.5909167
#> 13 600 - 800      0 1.4825069 1.4190417 1.5459721
#> 14 600 - 800      1 0.5573935 0.4916004 0.6231867
#> 15 600 - 800      2 1.0396711 0.9741703 1.1051719

We can now plot this however we choose using ggplot, such as

ggplot(plot_df, aes(x = Estimate, y = var2)) +
  geom_linerange(aes(xmin = lower, xmax = upper), size = 1,
                 color = "deepskyblue4", alpha = 0.5) +
  geom_point(size = 3, color = "deepskyblue4") +
  scale_y_discrete(labels = seq(-2, 2), expand = c(0.25, 0)) +
  labs(y = NULL) +
  facet_grid(var1~., switch = "y") +
  theme_minimal(base_size = 20) +
  theme(strip.placement = "outside",
        axis.text.y = element_text(size = 12),
        panel.spacing.y = unit(0, "mm"))

enter image description here

Or

ggplot(plot_df, aes(x = var1, y = Estimate, color = var2)) +
  geom_linerange(aes(ymin = lower, ymax = upper), size = 1,
                 position = position_dodge(width = 0.5)) +
  geom_point(position = position_dodge(width = 0.5), size = 3) +
  coord_flip() +
  labs(x = NULL, color = NULL) +
  scale_color_brewer(palette = "Set1", labels = seq(-2, 2)) +
  theme_minimal(base_size = 20)

enter image description here

Created on 2022-09-10 with reprex v2.0.2

Related