Skip to content

Instantly share code, notes, and snippets.

@brshallo
Last active August 28, 2020 13:01
Show Gist options
  • Save brshallo/f923b9b5c6360ce09beda35c3d1d55e9 to your computer and use it in GitHub Desktop.
Save brshallo/f923b9b5c6360ce09beda35c3d1d55e9 to your computer and use it in GitHub Desktop.
Given lm object, extract the high level variable, the levels, and the reference level, as well as the estimate of the impact for this variable (assuming using effects coding).
library(tidyverse)
library(broom)
library(car)

param <- getOption("contrasts")
go_deviance <- param
# traditional `contr.sum` does not name levels, so use function from `car` package
go_deviance["unordered"] <- "contr.Sum"

options(contrasts = go_deviance)

tidy_with_ref <- function(mod_obj, sep_regex = "\\[S\\."){
  
  tidy_coefs <- tidy(mod_obj) %>% 
    separate(term, c("variable", "level"), sep_regex, remove = FALSE) %>% 
    mutate(level = str_sub(level, 1, -2))
  
  xlevels <- mod_obj$xlevels  
  
  missing_levels <- xlevels %>% 
    enframe() %>% 
    unnest() %>% 
    set_names(c("variable", "level"))
  
  missing_value <- tidy_coefs %>% 
    filter(!is.na(level)) %>% 
    group_by(variable) %>% 
    summarise(estimate = sum(-estimate))
  
  missing_levels_values <- left_join(missing_levels, missing_value) %>% 
    anti_join(tidy_coefs, by = c("variable", "level"))
  
  bind_rows(tidy_coefs, missing_levels_values) %>% 
    arrange(variable, level)
}

mod_fit <- glm(cty ~ manufacturer + hwy + class, data = mpg)

tidy_with_ref(mod_fit)
#> Warning: Expected 2 pieces. Missing pieces filled with `NA` in 2 rows [1,
#> 16].
#> Joining, by = "variable"
#> # A tibble: 24 x 7
#>    term           variable   level   estimate std.error statistic   p.value
#>    <chr>          <chr>      <chr>      <dbl>     <dbl>     <dbl>     <dbl>
#>  1 (Intercept)    (Intercep~ <NA>    -0.813      0.535   -1.52     1.30e- 1
#>  2 class[S.2seat~ class      2seater -2.08       0.435   -4.77     3.47e- 6
#>  3 class[S.compa~ class      compact  0.00201    0.214    0.00939  9.93e- 1
#>  4 class[S.midsi~ class      midsize -0.681      0.199   -3.42     7.47e- 4
#>  5 class[S.miniv~ class      minivan  0.475      0.327    1.45     1.47e- 1
#>  6 class[S.picku~ class      pickup   1.47       0.236    6.21     2.81e- 9
#>  7 class[S.subco~ class      subcom~ -0.0172     0.224   -0.0767   9.39e- 1
#>  8 <NA>           class      suv      0.830     NA       NA       NA       
#>  9 hwy            hwy        <NA>     0.738      0.0217  34.0      1.00e-87
#> 10 manufacturer[~ manufactu~ audi    -0.972      0.272   -3.57     4.46e- 4
#> # ... with 14 more rows

Created on 2019-02-28 by the reprex package (v0.2.1)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment