Week 4: Exploratory Data Analysis and Visualization

FiveThirtyEight bad-drivers data

Week 4

Exploratory Data Analysis and Visualization

Tidyverse and useful EDA libraries

Agenda

• Intro: CUNY Data Science first Data Science seminar
• Due this week
• Due next week
• Data science in context
• EDA with the tidyverse and skimr
• Visualizing patterns, unusual values, and relationships
• Missing values and next steps

Due this week

Day	Item
Tuesday	Meetup, 8:00 p.m. ET
Thursday	Project 1 Approach due
Sunday	Scenario Design Quiz due
Sunday	Scenario Design Discussion due
Sunday	Project 1 Code Base due
Monday	Project 1 Video Explainer due

Due next week

Day	Item
Thursday	5A Airline Delays: Approach
Thursday	5B ELO Calculations: Approach
Sunday	Discussion 5A: Untidy Dataset
Sunday	Discussion/Quiz 5B: AI Nightmare
Sunday	5A Airline Delays: Code Base
Sunday	5B ELO Calculations: Code Base
Monday	5A Airline Delays: Video Explainer
Monday	5B ELO Calculations: Video Explainer
Next Week	PROJECT 2 DUE

CUNY Data Science seminar

Tuesday, 8:00 p.m. ET

Today’s data

FiveThirtyEight’s bad-drivers.csv contains one row for each state and the District of Columbia.

It includes fatal-collision rates, shares involving speeding and alcohol impairment, plus insurance premiums and losses.

Source: https://github.com/fivethirtyeight/data/tree/master/bad-drivers

The EDA mindset

Exploratory data analysis helps us learn what is in a data set before we model it.

• What does each row represent?
• Which variables need cleaning?
• Are values missing, unusual, or implausible?
• Which comparisons deserve a closer look?

Cleaning and renaming columns

clean_names() first converts the original labels to lowercase snake case. rename() then gives the measures short, readable names we can use throughout the lesson.

drivers <- read_csv("https://raw.githubusercontent.com/fivethirtyeight/data/refs/heads/master/bad-drivers/bad-drivers.csv", show_col_types = FALSE) |>
  clean_names() |>
  rename(
    fatal_collisions = number_of_drivers_involved_in_fatal_collisions_per_billion_miles,
    speeding = percentage_of_drivers_involved_in_fatal_collisions_who_were_speeding,
    alcohol_impaired = percentage_of_drivers_involved_in_fatal_collisions_who_were_alcohol_impaired,
    not_distracted = percentage_of_drivers_involved_in_fatal_collisions_who_were_not_distracted,
    no_previous_accidents = percentage_of_drivers_involved_in_fatal_collisions_who_had_not_been_involved_in_any_previous_accidents,
    premium = car_insurance_premiums,
    losses = losses_incurred_by_insurance_companies_for_collisions_per_insured_driver
  )

Result: the cleaned data

glimpse() shows the result of the full import, cleaning, and renaming pipeline: 51 rows and eight readable variables.

drivers |> glimpse()

Rows: 51
Columns: 8
$ state                 <chr> "Alabama", "Alaska", "Arizona", "Arkansas", "Cal…
$ fatal_collisions      <dbl> 18.8, 18.1, 18.6, 22.4, 12.0, 13.6, 10.8, 16.2, …
$ speeding              <dbl> 39, 41, 35, 18, 35, 37, 46, 38, 34, 21, 19, 54, …
$ alcohol_impaired      <dbl> 30, 25, 28, 26, 28, 28, 36, 30, 27, 29, 25, 41, …
$ not_distracted        <dbl> 96, 90, 84, 94, 91, 79, 87, 87, 100, 92, 95, 82,…
$ no_previous_accidents <dbl> 80, 94, 96, 95, 89, 95, 82, 99, 100, 94, 93, 87,…
$ premium               <dbl> 784.55, 1053.48, 899.47, 827.34, 878.41, 835.50,…
$ losses                <dbl> 145.08, 133.93, 110.35, 142.39, 165.63, 139.91, …

A readable data dictionary

Variable	Meaning
fatal_collisions	Drivers involved in fatal collisions per billion miles
speeding	Percent of those drivers who were speeding
alcohol_impaired	Percent of those drivers who were alcohol-impaired
premium	Average car-insurance premium in dollars
losses	Insurer losses per insured driver in dollars

A compact EDA report with skimr

skim() reports missing values, ranges, and distribution summaries. It is often the fastest first EDA step.

drivers |> skim()

Data summary

Name	drivers
Number of rows	51
Number of columns	8
_______________________
Column type frequency:
character	1
numeric	7
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
state	0	1	4	20	0	51	0

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
fatal_collisions	0	1	15.79	4.12	5.90	12.75	15.60	18.50	23.90	▁▇▇▇▃
speeding	0	1	31.73	9.63	13.00	23.00	34.00	38.00	54.00	▅▃▇▅▁
alcohol_impaired	0	1	30.69	5.13	16.00	28.00	30.00	33.00	44.00	▁▃▇▃▁
not_distracted	0	1	85.92	15.16	10.00	83.00	88.00	95.00	100.00	▁▁▁▂▇
no_previous_accidents	0	1	88.73	6.96	76.00	83.50	88.00	95.00	100.00	▃▆▇▅▆
premium	0	1	886.96	178.30	641.96	768.43	858.97	1007.94	1301.52	▆▇▂▃▂
losses	0	1	134.49	24.84	82.75	114.64	136.05	151.87	194.78	▂▆▇▆▂

Numeric summaries

summarise() collapses many rows into a few summary values.

drivers |>
  summarise(
    states = n(),
    mean_collision_rate = mean(fatal_collisions),
    median_collision_rate = median(fatal_collisions),
    min_collision_rate = min(fatal_collisions),
    max_collision_rate = max(fatal_collisions)
  )

# A tibble: 1 × 5
  states mean_collision_rate median_collision_rate min_collision_rate
   <int>               <dbl>                 <dbl>              <dbl>
1     51                15.8                  15.6                5.9
# ℹ 1 more variable: max_collision_rate <dbl>

Ordering states

arrange(desc(...)) sorts from largest to smallest. Use it to inspect extremes, not to draw a causal conclusion.

drivers |>
  select(state, fatal_collisions) |>
  arrange(desc(fatal_collisions)) |>
  slice_head(n = 8)

# A tibble: 8 × 2
  state          fatal_collisions
  <chr>                     <dbl>
1 North Dakota               23.9
2 South Carolina             23.9
3 West Virginia              23.8
4 Arkansas                   22.4
5 Kentucky                   21.4
6 Montana                    21.4
7 Louisiana                  20.5
8 Oklahoma                   19.9

Collision-rate distribution

The histogram counts states in rate ranges. The density plot shows the same distribution as a smooth curve. patchwork places the two plots side by side.

plot_histogram <- drivers |>
  ggplot(aes(x = fatal_collisions)) +
  geom_histogram(binwidth = 2, boundary = 0, fill = "navy", color = "white") +
  labs(
    title = "Histogram",
    x = "Fatal-collision rate\nper billion miles",
    y = "Number of states"
  ) +
  theme_minimal(base_size = 16)

plot_density <- drivers |>
  ggplot(aes(x = fatal_collisions)) +
  geom_density(fill = "navy", color = "navy", alpha = 0.7) +
  labs(
    title = "Density plot",
    x = "Fatal-collision rate\nper billion miles",
    y = "Density"
  ) +
  theme_minimal(base_size = 16)

plot_collision_distribution <- plot_histogram + plot_density +
  plot_layout(ncol = 2)

Collision-rate distribution: rendered plot

Comparing all states

plot_state_comparison <- drivers |>
  ggplot(aes(x = reorder(state, fatal_collisions), y = fatal_collisions)) +
  geom_col(fill = "navy") +
  coord_flip() +
  labs(
    title = "Fatal-collision rates by state",
    x = NULL,
    y = "Drivers involved in fatal collisions per billion miles"
  ) +
  theme_minimal(base_size = 12)

Comparing all states: rendered plot

A relationship to explore

Does the share of drivers who were speeding rise with the fatal-collision rate?

plot_speeding_collision <- drivers |>
  ggplot(aes(x = speeding, y = fatal_collisions)) +
  geom_point(size = 2.5, alpha = 0.75, color = "navy") +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  labs(
    title = "Speeding and fatal-collision rates",
    x = "Drivers involved in fatal collisions who were speeding (%)",
    y = "Fatal collisions per billion miles"
  ) +
  theme_minimal(base_size = 16)

Speeding and fatal-collision rates: rendered plot

Label unusual observations

Labels can make a plot easier to discuss. Here, geom_label_repel() labels only the three states with the largest collision rates and prevents overlapping labels.

median_speeding <- median(drivers$speeding)

plot_labeled_states <- drivers |>
  ggplot(aes(x = speeding, y = fatal_collisions)) +
  geom_point(size = 2.5, color = "navy") +
  geom_vline(
    xintercept = median_speeding,
    linetype = "dotted",
    color = "red"
  ) +
  annotate(
    "label",
    x = median_speeding + 1,
    y = 7,
    label = paste("Median speeding =", median_speeding, "%"),
    fill = "white",
    color = "red",
    hjust = 0
  ) +
  geom_label_repel(
    data = drivers |> filter(fatal_collisions >= 23),
    aes(label = state),
    fill = "white",
    seed = 607
  ) +
  labs(
    title = "States with the highest fatal-collision rates",
    x = "Speeding (%)",
    y = "Fatal collisions per billion miles"
  ) +
  theme_minimal(base_size = 16)

Labeled observations: rendered plot

Alcohol impairment and collision rates

plot_alcohol_collision <- drivers |>
  ggplot(aes(x = alcohol_impaired, y = fatal_collisions)) +
  geom_point(size = 2.5, alpha = 0.75, color = "navy") +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  labs(
    title = "Alcohol impairment and fatal-collision rates",
    x = "Drivers involved in fatal collisions who were alcohol-impaired (%)",
    y = "Fatal collisions per billion miles"
  ) +
  theme_minimal(base_size = 16)

Alcohol impairment and collision rates: rendered plot

Insurance premiums and losses

plot_premium_losses <- drivers |>
  ggplot(aes(x = premium, y = losses)) +
  geom_point(size = 2.5, alpha = 0.75, color = "navy") +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  scale_x_continuous(labels = scales::dollar) +
  scale_y_continuous(labels = scales::dollar) +
  labs(
    title = "Insurance premiums and insurer losses",
    x = "Average car-insurance premium",
    y = "Insurer losses per insured driver"
  ) +
  theme_minimal(base_size = 16)

Insurance premiums and losses: rendered plot

Stitching relationship plots with GGally

GGally::ggpairs() stitches scatterplots, distributions, and correlations into one grid. Use it to screen several relationships, then make one focused chart for your audience.

plot_relationship_grid <- drivers |>
  select(fatal_collisions, speeding, alcohol_impaired, premium, losses) |>
  ggpairs(
    lower = list(continuous = wrap("points", color = "navy")),
    diag = list(continuous = wrap("densityDiag", fill = "navy", color = "navy"))
  )

GGally relationship grid: rendered plot

Ridge plots with ggridges

Ridge plots compare the distribution of one numeric measure across groups. Here, each ridge shows fatal-collision rates for a broad U.S. region.

region_lookup <- tibble(
  state = c(state.name, "District of Columbia"),
  region = c(as.character(state.region), "South")
)

plot_ridges <- drivers |>
  left_join(region_lookup, by = "state") |>
  ggplot(aes(x = fatal_collisions, y = region)) +
  geom_density_ridges(fill = "navy", color = "navy", alpha = 0.7) +
  labs(
    title = "Fatal-collision rates by broad U.S. region",
    x = "Drivers involved in fatal collisions per billion miles",
    y = NULL
  ) +
  theme_minimal(base_size = 16)

Ridge plots: rendered plot

Correlation is a summary, not an explanation

drivers |>
  select(-state) |>
  cor() |>
  round(2)

                      fatal_collisions speeding alcohol_impaired not_distracted
fatal_collisions                  1.00    -0.03             0.20           0.01
speeding                         -0.03     1.00             0.29           0.13
alcohol_impaired                  0.20     0.29             1.00           0.04
not_distracted                    0.01     0.13             0.04           1.00
no_previous_accidents            -0.02     0.01            -0.25          -0.20
premium                          -0.20     0.04            -0.02           0.02
losses                           -0.04    -0.06            -0.08          -0.06
                      no_previous_accidents premium losses
fatal_collisions                      -0.02   -0.20  -0.04
speeding                               0.01    0.04  -0.06
alcohol_impaired                      -0.25   -0.02  -0.08
not_distracted                        -0.20    0.02  -0.06
no_previous_accidents                  1.00    0.08   0.04
premium                                0.08    1.00   0.62
losses                                 0.04    0.62   1.00

A correlation describes how two numeric variables move together in this data set. It does not establish why they move together.

Correlation heat map

A heat map uses color to show the direction and strength of correlations. GGally::ggcorr() draws the plot directly from a correlation matrix. We keep only one triangular half, since every correlation otherwise appears twice.

correlation_matrix <- drivers |>
  select(-state) |>
  cor()

correlation_matrix[upper.tri(correlation_matrix, diag = TRUE)] <- NA

plot_correlation_heatmap <- GGally::ggcorr(
  data = drivers |> select(-state),
  cor_matrix = correlation_matrix,
  low = "red",
  mid = "white",
  high = "navy",
  label = TRUE,
  label_round = 2,
  label_size = 3
  )

Correlation heat map: rendered plot

Outliers in data science

An outlier is an observation that falls unusually far from the other observations in a variable.

• An outlier is not automatically an error.
• It can reflect a data-entry issue, an unusual case, or a real pattern worth investigating.
• Before removing or changing one, check the data source and consider how it affects the question you are studying.

Faceted boxplots for possible outliers

pivot_longer() places selected numeric variables into one column. facet_wrap() creates a separate small panel, called a facet, for each variable.

• Each facet contains one boxplot.
• The dark points beyond the whiskers are possible outliers to investigate.
• scales = "free_y" lets each panel use a readable y-axis even though the variables have different units.

plot_numeric_outliers <- drivers |>
  select(-state, -not_distracted) |>
  pivot_longer(
    cols = everything(),
    names_to = "measure",
    values_to = "value"
  ) |>
  mutate(measure = str_replace_all(measure, "_", " ")) |>
  ggplot(aes(y = value)) +
  geom_boxplot(fill = "navy") +
  facet_wrap(~ measure, , ncol = 4) +
  labs(
    title = "Boxplots of selected numeric variables",
    x = NULL,
    y = NULL
  ) +
  theme_minimal(base_size = 13) +
  theme(
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank()
  )

Numeric variables: rendered boxplots

EDA workflow

1. Identify the unit of observation and the question.
2. Inspect names, types, and missingness.
3. Summarize variables and inspect extremes.
4. Visualize distributions and relationships.
5. Record data-quality questions before modeling.

Practice

Using drivers, create a plot that answers one question you care about.

• Choose two numeric variables.
• Make a scatterplot with clear axis labels.
• Add a trend line only if it helps your question.
• Write one sentence that distinguishes an observed pattern from a causal claim.

Missing values in EDA

Missing values can change summaries, plots, and models. Explore them before choosing an imputation method.

Load data with missing values

airquality is built into R. It records daily New York air-quality measures and already contains missing values in Ozone and Solar.R.

airquality_data <- airquality |>
  as_tibble()

Count missing values

airquality_data |>
  miss_var_summary()

# A tibble: 6 × 3
  variable n_miss pct_miss
  <chr>     <int>    <num>
1 Ozone        37    24.2 
2 Solar.R       7     4.58
3 Wind          0     0   
4 Temp          0     0   
5 Month         0     0   
6 Day           0     0   

Visualize missing values

naniar::vis_miss() displays observed and missing cells. Dark cells mark missing values.

plot_missingness <- airquality_data |>
  select(Ozone, Solar.R, Wind, Temp) |>
  vis_miss()

Missing values: rendered plot

Common imputation techniques

Technique	R function or package	Useful when	Main caution
Remove incomplete rows	drop_na()	Only a small amount is missing	Can discard useful data or introduce bias
Mean or median replacement	replace_na()	A simple numeric baseline is needed	Makes the data look less variable than it is
Mode or "Unknown" category	replace_na()	A categorical value is missing	Can distort category proportions
Forward or backward fill	tidyr::fill()	Rows have a meaningful order, such as time	Do not use for unrelated observations
k-nearest-neighbors imputation	VIM::kNN()	Similar observations can help estimate a value	Results depend on scaling and the distance rule
Multiple imputation	mice() and complete()	An analysis needs uncertainty from missingness reflected	Requires careful assumptions and diagnostics

EDA toolkit recap

• Tidyverse: import, clean, reshape, and summarize data with readr, dplyr, and tidyr.
• Data quality: use skimr for quick profiles and naniar to explore missing values.
• ggplot2: make rich visualizations - distributions, bar charts, scatterplots, boxplots, facets, and trend lines.ggrepel makes labels readable.
• Visualization extensions: patchwork combines plots; GGally explores pairs and correlations; ggridges compares group densities;