Memo: Visualize Regression Models

R for Data Science

Dependencies

library(tidyverse)
library(sjPlot) 

The package sjPlot provides the function plot_model. It plots model objects. Learn more here: https://strengejacke.github.io/sjPlot/index.html

1 Linear regression

Basic model with summary() view:

iris %>% 
  lm(formula = Sepal.Length ~ Sepal.Width + Petal.Length + Species, data = .) %>% 
  summary()

Call:
lm(formula = Sepal.Length ~ Sepal.Width + Petal.Length + Species, 
    data = .)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.82156 -0.20530  0.00638  0.22645  0.74999 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        2.39039    0.26227   9.114 5.94e-16 ***
Sepal.Width        0.43222    0.08139   5.310 4.03e-07 ***
Petal.Length       0.77563    0.06425  12.073  < 2e-16 ***
Speciesversicolor -0.95581    0.21520  -4.442 1.76e-05 ***
Speciesvirginica  -1.39410    0.28566  -4.880 2.76e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3103 on 145 degrees of freedom
Multiple R-squared:  0.8633,    Adjusted R-squared:  0.8595 
F-statistic: 228.9 on 4 and 145 DF,  p-value: < 2.2e-16

Plot the coefficients instead:

iris %>% 
  lm(formula = Sepal.Length ~ Sepal.Width + Petal.Length + Species, data = .) %>% 
  plot_model(model = .)

Spruce it up:

iris %>% 
  lm(formula = Sepal.Length ~ Sepal.Width + Petal.Length + Species, data = .) %>% 
  plot_model(model = ., color = "black", sort.est = TRUE, show.values = TRUE) +
  geom_hline(yintercept = 0, color = "gray") + 
  labs(title = "Look at this model", subtitle = "Wow such insight") + 
  theme_sjplot2() + 
  theme_light()

2 Logitistic Regression

HDMA data (you need to install the package AER):

data("HMDA", package = "AER")

Basic model with summary() view:

HMDA %>% 
  glm(formula = deny ~ pirat + mhist + chist + hirat + afam, data = ., family = binomial(link = "logit")) %>% 
  summary()

Call:
glm(formula = deny ~ pirat + mhist + chist + hirat + afam, family = binomial(link = "logit"), 
    data = .)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.2199  -0.4964  -0.3530  -0.2685   2.9264  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -5.0159     0.3258 -15.397  < 2e-16 ***
pirat         5.6869     0.9885   5.753 8.75e-09 ***
mhist2        0.5426     0.1750   3.101 0.001929 ** 
mhist3        0.5873     0.4589   1.280 0.200671    
mhist4        0.7653     0.5683   1.347 0.178067    
chist2        0.7350     0.1925   3.818 0.000135 ***
chist3        1.0396     0.2760   3.767 0.000165 ***
chist4        1.2909     0.3156   4.090 4.31e-05 ***
chist5        1.3157     0.2191   6.005 1.91e-09 ***
chist6        1.7262     0.2013   8.576  < 2e-16 ***
hirat        -0.7388     1.1525  -0.641 0.521476    
afamyes       0.8709     0.1582   5.504 3.70e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1744.2  on 2379  degrees of freedom
Residual deviance: 1480.9  on 2368  degrees of freedom
AIC: 1504.9

Number of Fisher Scoring iterations: 5

Plot the odds-ratios instead:

HMDA %>% 
  glm(formula = deny ~ pirat + mhist + chist + hirat + afam, data = ., family = binomial(link = "logit")) %>% 
  plot_model(model = .) 

Fine-tune the plot:

HMDA %>% 
  glm(formula = deny ~ pirat + mhist + chist + hirat + afam, data = ., family = binomial(link = "logit")) %>% 
  plot_model(model = ., color = "black", sort.est = TRUE) +
  geom_hline(yintercept = 1) + 
  labs(title = "Probability a mortgage application is denied", subtitle = "Boston Housing Data, 1990") + 
  theme_minimal()

Plot predicted probabilities of a continuous variable:

HMDA %>% 
  glm(formula = deny ~ pirat + mhist + chist + hirat + afam, data = ., family = binomial(link = "logit")) %>% 
  plot_model(model = ., type = "pred", terms="pirat [all]")