Coefficient plots for linear models

Coefficient plots are often more useful than tables¹ but plotting raw coefficients can be misleading when the predictors are on different scales.

The packages arm and modelsummary, and later, ggstats are used to illustrate these plots.

In the process, I discover some other problems with naive use of coefficient plots. I compare plotting:

• raw coefficients
• standardized coefficients
• meaningfully scaled coefficients

Load packages

library(dplyr)
library(ggplot2)
library(arm)           # coefplot()
library(modelsummary)  # modelplot()
library(ggstats)       # ggcoef_model()

Occupational prestige data

How does occupational prestige depend on education (years), income ($), %women, …?

I use a classic example from the carData package. The observations are averages for occupational categories rather than individuals. The data come from the 1971 Canadian census.

data(Prestige, package="carData")
head(Prestige, 5)

##                     education income women prestige census type
## gov.administrators      13.11  12351 11.16     68.8   1113 prof
## general.managers        12.26  25879  4.02     69.1   1130 prof
## accountants             12.77   9271 15.70     63.4   1171 prof
## purchasing.officers     11.42   8865  9.11     56.8   1175 prof
## chemists                14.62   8403 11.68     73.5   2111 prof

Fit a simple model

mod0 <- lm(prestige ~ education + income + women,
           data=Prestige)
summary(mod0)

## 
## Call:
## lm(formula = prestige ~ education + income + women, data = Prestige)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -19.825  -5.333  -0.136   5.159  17.504 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6.794334   3.239089   -2.10    0.039 *  
## education    4.186637   0.388701   10.77  < 2e-16 ***
## income       0.001314   0.000278    4.73  7.6e-06 ***
## women       -0.008905   0.030407   -0.29    0.770    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.85 on 98 degrees of freedom
## Multiple R-squared:  0.798,  Adjusted R-squared:  0.792 
## F-statistic:  129 on 3 and 98 DF,  p-value: <2e-16

Displaying coefficients

I’m interested in the coefficients in this model, so here is just a reminder that you can extract the numerical values using coef(), or display them in a table similar to the “Coefficients” portion above using broom::tidy() or arm::display():

coef(mod0)

## (Intercept)   education      income       women 
##   -6.794334    4.186637    0.001314   -0.008905

broom::tidy(mod0)

## # A tibble: 4 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept) -6.79     3.24        -2.10  3.85e- 2
## 2 education    4.19     0.389       10.8   2.59e-18
## 3 income       0.00131  0.000278     4.73  7.58e- 6
## 4 women       -0.00891  0.0304      -0.293 7.70e- 1

arm::display(mod0, detail = TRUE)

## lm(formula = prestige ~ education + income + women, data = Prestige)
##             coef.est coef.se t value Pr(>|t|)
## (Intercept) -6.79     3.24   -2.10    0.04   
## education    4.19     0.39   10.77    0.00   
## income       0.00     0.00    4.73    0.00   
## women       -0.01     0.03   -0.29    0.77   
## ---
## n = 102, k = 4
## residual sd = 7.85, R-Squared = 0.80

Visualize coefficients

The plots below show default coefficient plots for this model using modelsummary::modelplot() and arm::coefplot(). At the end, I show some examples using the ggstats package.

They are disappointing and misleading because these show the raw coefficients. It looks like only education has a non-zero effect, but the effect of income is also highly significant.

These plots are misleading because the predictors are on different scales, so it makes little sense to compare the change in prestige for a 1 year increase in education with the change for a $1 increase in income.

theme_set(theme_minimal(base_size = 14))
modelplot(mod0, 
          coef_omit="Intercept", 
          color="blue", size=1) +
  labs(title="Raw coefficients for mod0")

arm::coefplot(mod0, col.pts="red", cex.pts=1.5)

Other graphical features to note here:

• coefplot() includes a dashed vertical reference line at x = 0 by default, but modelplot() does not.
• The \(Y\) axis is ordered in reverse to what appears in tabular output.

Fit a more complex model

For this example, graphical analysis suggested that the relation between prestige and income was non-linear, which could be corrected using log(income). As well, other plots suggested an interaction with type of occupation.

mod1 <- lm(prestige ~ education + women +
             log(income)*type, data=Prestige)
summary(mod1)

## 
## Call:
## lm(formula = prestige ~ education + women + log(income) * type, 
##     data = Prestige)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.814  -4.685   0.656   3.966  16.569 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -152.2059    23.2499   -6.55  3.5e-09 ***
## education               2.9282     0.5883    4.98  3.1e-06 ***
## women                   0.0883     0.0323    2.73   0.0076 ** 
## log(income)            18.9819     2.8285    6.71  1.7e-09 ***
## typeprof               85.2642    30.4582    2.80   0.0063 ** 
## typewc                 29.4133    36.5075    0.81   0.4226    
## log(income):typeprof   -9.0124     3.4102   -2.64   0.0097 ** 
## log(income):typewc     -3.8334     4.2603   -0.90   0.3706    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.27 on 90 degrees of freedom
##   (4 observations deleted due to missingness)
## Multiple R-squared:  0.875,  Adjusted R-squared:  0.865 
## F-statistic: 90.1 on 7 and 90 DF,  p-value: <2e-16

Compare models

arm::coefplot() allows you to compare several models in the same plot, using add=TRUE for the second and subsequent ones. The intercept is ignored by default.

arm::coefplot(mod0, col.pts="red", cex.pts=1.5)
arm::coefplot(mod1, add=TRUE, col.pts="blue", cex.pts=1.5)

But WAIT: income was entered as log(income) in mod1. The plot above is silently wrong, because it doesn’t include log(income) or the interaction terms. It is was probably stupid to try to compare the coefficients, and this should have raised a warning or error.

All the terms will appear if the most complex model is plotted first.

arm::coefplot(mod1, col.pts="red", cex.pts=1.5)
arm::coefplot(mod0, add=TRUE, col.pts="blue", cex.pts=1.5)

Standardize the variables, giving \(\beta\) coefficients

An alternative is to present the standardized coefficients. These are interpreted as the standardized change in prestige for a one standard deviation change in the predictors.

One way to do this is to transform all variables to standardized (\(z\)) scores. The syntax ultimately uses scale to transform all the numeric variables.

Prestige_std <- Prestige |>
  as_tibble() |>
  mutate(across(where(is.numeric), scale))

Fit the same models to the standardized variables

This accords better with the fitted (standardized) coefficients. income and education are now both large and significant compared with the effect of women.

But a nagging doubt remains: How can I interpret the numerical size of coefficients? The direct answer is that they give the expected change in standard deviations of prestige for a one standard deviation change in each of the predictors.

mod0_std <- lm(prestige ~ education + income + women,
               data=Prestige_std)

mod1_std <- lm(prestige ~ education + women +
               log(income)*type, 
               data=Prestige_std)

Plot the standardized coefficients

These plots look more like what I was expecting to show the relative magnitudes of the coefficients and CIs so one could see which differ from zero.

modelplot(mod1_std, 
          coef_omit="Intercept", color="blue", size=1) +
  labs(title="Standardized coefficients for mod0") +
  geom_vline(xintercept = 0, linetype = 2) 

Using standardize = "refit"

It turns out there is an easier way to get plots of standardized coefficients. modelsummary() extracts coefficients from model objects using the parameters package, and that package offers several options for standardization: See model parameters documentation. We can pass the standardize="refit" (or other) argument directly to modelsummary() or modelplot(), and that argument will be forwarded to parameters. To compare models, pass them to modelplot as a list.

modelplot(list(mod0 = mod0, mod1 = mod1),
          standardize = "refit",
          coef_omit="Intercept", size=1) +
  labs(title="Standardized coefficients for mod0 and mod1") +
  geom_vline(xintercept = 0, linetype = 2) 

A small annoyance: This gives a warning because standardize is handled via ... and that is pushed to both parameters and ggplot2. ## More meaningful units A better substantative comparison of the coefficients could use understandable scales for the predictors, e.g., months of education, $50,000 income or 10% of women’s participation.

Note that the effect of this is just to multiply the coefficients and their standard errors by a factor. The statistical conclusions of significance are unchanged.

Prestige_scaled <- Prestige |>
  mutate(education = 12 * education,
         income = income / 50,
         women = women / 10)

mod0_scaled <- lm(prestige ~ education + income + women,
               data=Prestige_scaled)

arm::display(mod0_scaled, detail = TRUE)

## lm(formula = prestige ~ education + income + women, data = Prestige_scaled)
##             coef.est coef.se t value Pr(>|t|)
## (Intercept) -6.79     3.24   -2.10    0.04   
## education    0.35     0.03   10.77    0.00   
## income       0.07     0.01    4.73    0.00   
## women       -0.09     0.30   -0.29    0.77   
## ---
## n = 102, k = 4
## residual sd = 7.85, R-Squared = 0.80

arm::coefplot(mod0_scaled, col.pts="red", cex.pts=1.5,
              main = "Coefficients for prestige with scaled predictors",
              varnames=c("intercept", 
                         "education\n(/month)",
                         "income\n(/$50K)",
                         "women\n(/10%)")
               )

Re-ordering terms

There is a mismatch between the tabular displays of coefficients, where terms are listed in their order in the model, and the plots shown above, where terms appear on the vertical axis in the reverse order, because the \(Y\) axis starts at 0. With modelplot(), we can reorder the terms using scale_y_discrete().

modelplot(list(mod0 = mod0, mod1 = mod1),
          standardize = "refit",
          coef_omit="Intercept", size=1) +
  labs(title="Standardized coefficients for mod0 and mod1") +
  scale_y_discrete(limits = rev) +
  geom_vline(xintercept = 0, linetype = 2) +
  theme(legend.position = c(0.9, 0.2))

Factor levels

There are several features of the plots shown above that are handled better in other packages. A main one is how terms involving factor variables are handled: In the Prestige data, occupation type is a factor, with levels "bc", "prof", "wc" meaning blue-collar, professional and white-collar jobs.

By default, lm() treats the first level ("bc") as the baseline category, so the coefficients for type are labeled typeprof and typewc. This is ugly and non-informative for presentation purposes, where it would be better to group the coefficients for the levels under the type factor.

To illustrate, I ignore using log(income), and specify some models using type and its interactions.

mod0 <- lm(prestige ~ education + income + women,
          data=Prestige_std)
mod2 <- update(mod0, ~ . + type)
mod3 <- update(mod0, ~ . + type + type * income)

The ggstats package provides some nicer versions of coefficient plots that handle factors in a more reasonable way, as levels nested within the factor. By default it shows \(p\)-values and significance stars in the labels, and (redundantly) uses different symbols for terms significant at \(p \le 0.05\) or not.

library(ggstats)
ggcoef_model(mod2) +
  labs(x = "Standarized Coefficient")

I find the \(p\)-values and significance stars distracting in the plot, but these can be suppressed using options show_p_values = FALSE and signif_stars = FALSE.

Comparing models

To compare several models, use ggcoef_compare(), supplied with a list of models. This does a much nicer job of organizing the terms for main effects and interactions so they may be more readily compared.

models <- list(
  "Base model"      = mod0,
  "Add type"        = mod2,
  "Add interaction" = mod3)

ggcoef_compare(models) +
  labs(x = "Standarized Coefficient")

Wrapup

At first glance, printing and plotting coefficients from statistical models seems like an easy task for which there should be straightforward solutions, with reasonable defaults and the ability to customize displays easily.

modelplot() is nice in that you can ask it to plot standardized coefficients, but it doesn’t handle factors well. I think that, by default it should reverse the order of coefficients on the vertical axis and plot the vertical reference line at 0.

---

1. Kastellec, J. P., & Leoni, E. L. (2007). Using Graphs Instead of Tables in Political Science. Perspectives on Politics, 5, 755–771.↩︎