---
layout: post
title: "Plotting a regression line on a scatter plot of smoking and drinking data with ggplot2 (CC355)"
blurb: "Measuring the relationship between smoking and excessive drinking"
author: "PD Schloss"
date: 2025-04-17 8:00
comments: false
youtube: oWiKv_p3ktM
start_hash: 
end_hash: 
---

Pat uses R to recreate a plot from Philip Bump showing the correlation between smoking and drinking across US counties using functions from the tidyverse. To pull this off he uses the dplyr, ggplot2, and showtext packages. The functions he uses from these packages include aes, annotate, coord_cartesian, drop_na, element_blank, element_line, element_text, filter, font_add_google, geom_abline, geom_hline, geom_label, geom_point, geom_smooth, geom_vline, ggplot, ggsave, inner_join, labs, library, lm, median, mutate, read_csv, round, select, showtext_auto, showtext_opts, starts_with, summarize, theme, tribble, and unit. The newsletter describing this visualization at a 30,000 ft view can be found [here](https://shop.riffomonas.org/posts/how-do-you-even-pronounce-bivariate-choropleth-here-s-an-alternative). You can find the original article [here](https://s2.washingtonpost.com/camp-rw/?trackId=undefined&s=67dec26d37beb82d8ae1cdb7). The UWPHI data can be found [here](https://www.countyhealthrankings.org/health-data/methodology-and-sources/data-documentation). If you have a figure that you would like to see me discuss in a future newsletter and episode of Code Club, email me at pat@riffomonas.org!

<iframe style="margin: 0 auto;display:block;" width="560" height="315" src="https://www.youtube.com/embed/{{ page.youtube }}" frameborder="0" allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>


## Code

```R
library(tidyverse)
library(showtext)

font_add_google("Libre Franklin", "franklin")
showtext_opts(dpi = 300)
showtext_auto()

codes_of_interest <- tribble(
  ~fipscode, ~hjust, ~vjust,
  "02180",  -0.15,  0.5,
  "38079",  0.5,  -0.5, 
  "02070",  -0.15,  0.5,
  "19115",  -0.15,  0.5,
  "19019",  -0.15,  0.5,
  "55109",  0.5,  1.5,
  "06041",  0.5,  1.5,
  "51610",  0.5,  1.5,
  "49043",  0.5,  1.5,
  "49051",  0.5,  1.5,
  "49049",  1.15,  0.5,
  "09005",  0,  0,
  "24033",  1.15,  0.5,
  "49031",  1.15,  0.5,
  "28021",  1.15,  0.5,
  "28063",  1.15,  0.5,
  "28053",  1.15,  0.5,
  "21237",  1.15,  0.5,
  "46121",  1.15,  0.5,
  "38085",  -0.15,  0.5
)

smoking_drinking <- read_csv("https://www.countyhealthrankings.org/sites/default/files/media/document/analytic_data2025.csv",
    skip = 1) %>%
  filter(countycode != "000") %>%
  select(fipscode, state, county,
         smoking = v009_rawvalue, drinking = v049_rawvalue) %>%
  drop_na() %>%
  mutate(smoking = round(smoking, digits = 2),
         drinking = round(drinking, digits = 2))

median_values <- smoking_drinking %>%
  summarize(med_smoking = median(smoking),
            med_drinking = median(drinking))

interest <- inner_join(smoking_drinking, codes_of_interest, by = "fipscode")

model <- lm(smoking~drinking, data = smoking_drinking)
intercept <- model$coefficients[1]
slope <- model$coefficients[2]

smoking_drinking %>%
  ggplot(aes(x = drinking, y = smoking)) +
  geom_vline(xintercept = median_values$med_drinking, color = "gray80",
             linewidth = 0.2) +
  geom_hline(yintercept = median_values$med_smoking, color = "gray80",
             linewidth = 0.2) +
  geom_point(size = 1, alpha = 0.2, color = "#479126") +
  geom_point(data = interest, 
             aes(x = drinking, y = smoking), fill = "#499427", shape = 21) +
  geom_abline(slope = slope, intercept = intercept,
              linewidth = 0.4) +
  # geom_smooth(formula = y~x, method = "lm") +
  geom_label(data = interest, 
             aes(x = drinking, y = smoking, label = fipscode),
             hjust = interest$hjust, vjust = interest$vjust,
             family = "franklin", size = 8.5, size.unit = "pt",
             label.padding = unit(1, "pt"), label.size = unit(0, "pt")) +
  annotate(geom = "text",
           x = c(-0.03, 0.317),
           y = c(0.416, -0.035),
           hjust = c(-0.05, 1),
           vjust = c(1.2, -0.5),
           label = c("Adult smoking", "Excessive drinking"),
           family = "franklin", size = 8.5, size.unit = "pt") +
  annotate(geom = "text",
           x = c(0, 0.317),
           y = c(0.36, 0.01),
           hjust = c(0, 1),
           label = c("Less drinking,\nmore smokers",
                     "More drinking,\nfewer smokers"),
           family = "franklin", size = 9, size.unit = "pt",
           fontface = "bold", lineheight = 1) +
  coord_cartesian(clip = "off", expand = FALSE,
                  xlim = c(-0.03, 0.317), ylim = c(-0.035, 0.416)) +
  labs(
    title = "How smoking compares with excessive drinking",
    subtitle = "Median county values are shown with horizontal and vertical lines.",
    caption = "Source: University of Wisconsin Population Health Institute"
  ) +
  theme(
    text = element_text(family = "franklin"),
    panel.grid = element_blank(),
    panel.background = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    axis.title = element_blank(),
    axis.line = element_line(linewidth = 0.3),
    plot.title = element_text(face = "bold", margin = margin(t = 3, b = 12)),
    plot.subtitle = element_text(margin = margin(b = 8)),
    plot.caption = element_text(hjust = 0, margin = margin(t = 7))
  )

ggsave("uwphi_biplot.png", width = 6, height = 6.39)
```
