EPI 553 | Principles of Statistical Inference II
2026-04-23
Figure 1: Anscombe’s Quartet: four datasets with identical summary statistics but radically different patterns.
Every analysis you do should begin and end with looking at the data.
ggplot(penguins, aes(x = bill_length_mm, y = bill_depth_mm, color = species)) +
geom_point(alpha = 0.7, size = 2.5) +
geom_smooth(method = "lm", se = FALSE, linewidth = 1) +
facet_wrap(~ island) +
scale_color_manual(values = epi_colors) +
labs(title = "Bill Dimensions of Palmer Penguins",
subtitle = "By species and island, with linear trend lines",
x = "Bill length (mm)", y = "Bill depth (mm)",
color = "Species")Figure 2
Cleveland & McGill (1984) ranked visual encodings by accuracy:
Figure 4
Prefer dot plots and bar charts to pie charts. Prefer faceting to stacking. Use color for categories, not for quantities (unless you use a perceptually uniform scale like viridis).
dat <- tibble(group = LETTERS[1:6],
value = c(22, 18, 17, 16, 14, 13))
p1 <- ggplot(dat, aes(x = "", y = value, fill = group)) +
geom_col(width = 1) +
coord_polar("y") +
scale_fill_brewer(palette = "Set2") +
labs(title = "Pie chart") +
theme_void(base_size = 14) +
theme(plot.title = element_text(face = "bold", hjust = 0.5))
p2 <- ggplot(dat, aes(x = value, y = fct_reorder(group, value))) +
geom_col(fill = "steelblue", width = 0.7) +
geom_text(aes(label = value), hjust = -0.3, fontface = "bold") +
scale_x_continuous(expand = expansion(mult = c(0, 0.1))) +
labs(title = "Sorted bar chart", x = NULL, y = NULL) +
theme_epi553()
p1 + p2 +
plot_annotation(title = "Pie charts encode angles. Bar charts encode position.",
subtitle = "Which one lets you instantly rank the groups?",
theme = theme(plot.title = element_text(face = "bold", size = 16),
plot.subtitle = element_text(color = "grey40")))Figure 5
library(ggrepel)
penguin_means <- penguins |>
drop_na() |>
summarise(mass = mean(body_mass_g), .by = c(species, year))
ggplot(penguin_means, aes(x = year, y = mass, color = species)) +
geom_line(linewidth = 1.4) +
geom_point(size = 3.5) +
geom_text_repel(
data = filter(penguin_means, year == max(year)),
aes(label = species),
hjust = 0, nudge_x = 0.15, direction = "y",
segment.color = NA, fontface = "bold", size = 5
) +
scale_color_manual(values = epi_colors) +
scale_x_continuous(breaks = 2007:2009,
expand = expansion(mult = c(0.05, 0.25))) +
labs(title = "Average Body Mass of Palmer Penguins, 2007-2009",
subtitle = "Direct labels replace the legend; gridlines are softened",
x = NULL, y = "Body mass (g)",
caption = "Source: palmerpenguins R package") +
theme(legend.position = "none",
panel.grid.major.x = element_blank())Figure 6
ggplot(penguins, aes(x = species, y = body_mass_g)) +
geom_boxplot() +
theme_grey(base_size = 14)Figure 7
Honest. Also forgettable.
ggplot(penguins, aes(x = species, y = body_mass_g, fill = species)) +
geom_boxplot(alpha = 0.5, outlier.shape = NA) +
geom_jitter(width = 0.2, alpha = 0.4, size = 1.5) +
scale_fill_manual(values = epi_colors) +
theme_epi553() +
theme(legend.position = "none")Figure 8
Now we see how many penguins are in each species and where the outliers actually live.
ggplot(penguins, aes(x = species, y = body_mass_g, fill = species)) +
geom_boxplot(alpha = 0.5, outlier.shape = NA, color = "grey20") +
geom_jitter(width = 0.2, alpha = 0.4, color = "grey30", size = 1.5) +
scale_fill_manual(values = epi_colors) +
labs(title = "Body Mass of Palmer Penguins by Species",
subtitle = "Boxplots with individual observations overlaid",
x = NULL, y = "Body mass (g)") +
theme(legend.position = "none")Figure 9
library(ggdist)
ggplot(drop_na(penguins, body_mass_g),
aes(x = species, y = body_mass_g, fill = species)) +
stat_halfeye(adjust = 0.5, width = 0.6, .width = 0,
justification = -0.3, point_color = NA) +
geom_boxplot(width = 0.15, outlier.shape = NA, alpha = 0.7) +
geom_jitter(width = 0.05, alpha = 0.3, size = 1.2) +
scale_fill_manual(values = epi_colors) +
coord_cartesian(xlim = c(1.2, NA), clip = "off") +
labs(title = "Gentoo Penguins Are Substantially Heavier",
subtitle = "Distribution, median, and individual observations of body mass by species",
x = NULL, y = "Body mass (g)",
caption = "Source: palmerpenguins R package") +
theme(legend.position = "none",
panel.grid.major.x = element_blank())Figure 10
The final chart shows the distribution (raincloud), the summary (boxplot), the raw data (jitter), and a conclusion baked into the title. That is the difference between a chart and a finding.
Figure 11
Figure 12
Set2, Dark2, Okabe-Itoviridis, mako, BluesRdBu, BrBG (only when there is a meaningful zero)About 8% of men and 0.5% of women have some form of color vision deficiency.
Figure 14
Rule of thumb: use viridis-family palettes by default. Check with colorBlindness::cvdPlot() or coblis.com.
Many groups, but only a few matter to your story. Highlighting fades out the context.
library(gghighlight)
gapminder_like <- penguins |>
drop_na() |>
summarise(
mean_mass = mean(body_mass_g),
.by = c(species, year)
)
ggplot(gapminder_like,
aes(x = year, y = mean_mass,
color = species)) +
geom_line(linewidth = 1.4) +
geom_point(size = 3.5) +
gghighlight(
species == "Gentoo",
unhighlighted_params = list(
color = "grey80",
linewidth = 0.8)) +
scale_color_manual(
values = c("Gentoo" = "#2E86AB"))The eye is drawn instantly to the highlighted line, but the reader still has full context.
Boxplots hide bimodality. Histograms hide group differences. Raincloud plots show everything.
library(ggdist)
ggplot(drop_na(penguins, body_mass_g),
aes(x = body_mass_g, y = species, fill = species)) +
stat_halfeye(adjust = 0.6, .width = 0, justification = -0.2,
point_color = NA, alpha = 0.7) +
geom_boxplot(width = 0.15, outlier.shape = NA, alpha = 0.5) +
geom_jitter(height = 0.07, alpha = 0.3, size = 1.2) +
scale_fill_manual(values = epi_colors) +
labs(title = "Raincloud Plot of Penguin Body Mass",
subtitle = "Distribution + boxplot + raw observations in one chart",
x = "Body mass (g)", y = NULL) +
theme(legend.position = "none")Figure 18
When you want to compare distributions across many groups in a small space.
library(ggridges)
ggplot(diamonds, aes(x = price, y = cut, fill = cut)) +
geom_density_ridges(alpha = 0.8, scale = 1.1) +
scale_x_log10(labels = label_dollar()) +
scale_fill_viridis_d(option = "mako") +
labs(title = "Distribution of Diamond Prices by Cut",
subtitle = "Ridge plots make many distributions easy to compare",
x = "Price (log scale)", y = NULL) +
theme(legend.position = "none")Figure 19
Confidence intervals as skinny error bars often get missed. Show the whole distribution.
library(ggdist)
library(distributional)
estimates <- tibble(
predictor = c("Smoking", "Exercise", "Income", "Sleep", "Age"),
estimate = c(0.65, -0.45, -0.20, -0.30, 0.05),
se = c(0.10, 0.08, 0.07, 0.09, 0.04)
)
ggplot(estimates, aes(y = fct_reorder(predictor, estimate),
xdist = dist_normal(estimate, se))) +
stat_halfeye(.width = c(0.5, 0.95), fill = "#2E86AB",
slab_alpha = 0.6, point_size = 3.5) +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey40") +
labs(title = "Coefficient Plot with Full Sampling Distributions",
subtitle = "Thick band = 50% interval, thin line = 95% interval",
x = "Log odds ratio", y = NULL)Figure 20
This is far more honest than a forest plot of point estimates with whiskers, because the reader can see the shape of the uncertainty.
p_scatter <- ggplot(penguins, aes(bill_length_mm, body_mass_g, color = species)) +
geom_point(alpha = 0.7, size = 2) +
scale_color_manual(values = epi_colors) +
theme(legend.position = "none") +
labs(x = "Bill length (mm)", y = "Body mass (g)")
p_top <- ggplot(penguins, aes(bill_length_mm, fill = species)) +
geom_density(alpha = 0.6) +
scale_fill_manual(values = epi_colors) +
theme_void() + theme(legend.position = "none")
p_right <- ggplot(penguins, aes(body_mass_g, fill = species)) +
geom_density(alpha = 0.6) +
scale_fill_manual(values = epi_colors) +
coord_flip() +
theme_void() + theme(legend.position = "none")
(p_top + plot_spacer() + p_scatter + p_right +
plot_layout(ncol = 2, widths = c(4, 1), heights = c(1, 4))) +
plot_annotation(
title = "Marginal Density + Scatterplot",
subtitle = "Composed with patchwork: p_top + plot_spacer() + p_scatter + p_right",
theme = theme(plot.title = element_text(face = "bold", size = 16),
plot.subtitle = element_text(color = "grey40", family = "Fira Code",
size = 11))
)Figure 21
The fastest way to make a chart “editorial”: annotate the finding directly on the chart.
library(ggtext)
ggplot(penguins, aes(x = flipper_length_mm, y = body_mass_g, color = species)) +
geom_point(alpha = 0.7, size = 2.2) +
annotate("curve", x = 175, y = 5500, xend = 210, yend = 5100,
arrow = arrow(length = unit(0.25, "cm")),
curvature = -0.3, color = "grey30", linewidth = 0.6) +
annotate("label", x = 173, y = 5650,
label = "Gentoos cluster in the\nheavy / long-flipper corner",
hjust = 0, size = 4.2, color = "grey20", lineheight = 0.95,
fill = "#fff8e1", label.size = 0, fontface = "italic") +
scale_color_manual(values = epi_colors) +
labs(title = "Body mass scales with flipper length, but <span style='color:#46166B;'>**species matter**</span>",
subtitle = "Annotated with a callout arrow and label box",
x = "Flipper length (mm)", y = "Body mass (g)",
color = NULL) +
theme(plot.title = element_markdown(face = "bold", size = 16),
legend.position = "top")Figure 22
Public health data are almost always spatial. The sf package gives ggplot2 first-class support.
library(sf)
library(tigris)
options(tigris_use_cache = TRUE)
ny_counties <- counties(state = "NY", cb = TRUE, class = "sf")
ny_counties$fmd_prev <- runif(nrow(ny_counties), 8, 18)
ggplot(ny_counties) +
geom_sf(aes(fill = fmd_prev), color = "white", linewidth = 0.2) +
scale_fill_viridis_c(option = "rocket", direction = -1,
name = "FMD %") +
labs(title = "Frequent Mental Distress Prevalence by NY County",
subtitle = "Hypothetical data for demonstration",
caption = "Source: U.S. Census Bureau (TIGER/Line)") +
theme_void(base_size = 14) +
theme(plot.title = element_text(face = "bold"),
plot.subtitle = element_text(color = "grey40"))Figure 23
Heatmaps turn a matrix of numbers into a pattern you can see. Essential for correlation tables, confusion matrices, and time-by-group summaries.
# Correlation matrix of numeric penguin variables
cor_data <- penguins |>
drop_na() |>
select(where(is.numeric)) |>
cor() |>
as.data.frame() |>
rownames_to_column("var1") |>
pivot_longer(-var1, names_to = "var2", values_to = "cor")
ggplot(cor_data, aes(x = var1, y = var2, fill = cor)) +
geom_tile(color = "white", linewidth = 0.8) +
geom_text(aes(label = round(cor, 2),
color = abs(cor) > 0.6),
size = 4.5, fontface = "bold") +
scale_fill_gradient2(low = "#D55E00", mid = "white", high = "#0072B2",
midpoint = 0, limits = c(-1, 1),
name = "Correlation") +
scale_color_manual(values = c("TRUE" = "white", "FALSE" = "grey20"),
guide = "none") +
scale_x_discrete(labels = \(x) str_replace_all(x, "_", "\n")) +
scale_y_discrete(labels = \(x) str_replace_all(x, "_", "\n")) +
labs(title = "Correlation Heatmap of Penguin Measurements",
subtitle = "Color intensity encodes strength; text encodes exact value",
x = NULL, y = NULL) +
coord_fixed() +
theme(axis.text.x = element_text(angle = 0, hjust = 0.5, size = 10),
axis.text.y = element_text(size = 10),
panel.grid = element_blank())Figure 24
When you need to show the difference between two time points or conditions, dumbbell charts are more effective than grouped bars.
# Simulated epi example: disease rates before/after intervention
intervention_data <- tibble(
county = c("Albany", "Saratoga", "Rensselaer", "Schenectady",
"Columbia", "Greene", "Warren", "Washington"),
before = c(15.2, 12.8, 18.1, 16.5, 11.3, 14.7, 9.8, 13.2),
after = c(11.1, 10.2, 12.4, 13.8, 9.5, 11.0, 8.1, 10.9)
) |>
mutate(change = after - before,
county = fct_reorder(county, change))
ggplot(intervention_data) +
geom_segment(aes(x = before, xend = after,
y = county, yend = county),
color = "grey60", linewidth = 1.2) +
geom_point(aes(x = before, y = county), color = "#D55E00",
size = 4) +
geom_point(aes(x = after, y = county), color = "#0072B2",
size = 4) +
annotate("text", x = 19, y = 8.3, label = "Before",
color = "#D55E00", fontface = "bold", size = 4.5) +
annotate("text", x = 19, y = 7.7, label = "After",
color = "#0072B2", fontface = "bold", size = 4.5) +
labs(title = "Every County Improved After the Intervention",
subtitle = "Rate per 1,000 population, before vs. after community health program",
x = "Rate per 1,000", y = NULL) +
scale_x_continuous(limits = c(7, 20))Figure 25
Dumbbell charts encode direction, magnitude, and rank simultaneously. Far more effective than grouped bar charts for before/after comparisons.
A wall of numbers:
| Term | Estimate | SE | p |
|---|---|---|---|
| (Intercept) | -1.23 | 0.41 | 0.003 |
| Smoking | 0.65 | 0.10 | <0.001 |
| Exercise | -0.45 | 0.08 | <0.001 |
| Income | -0.20 | 0.07 | 0.004 |
| Sleep | -0.30 | 0.09 | 0.001 |
| Age | 0.05 | 0.04 | 0.211 |
A coefficient plot communicates direction, magnitude, uncertainty, and significance in one glance. A table requires row-by-row mental math.
Fit a logistic regression on the penguins data and plot the odds ratios directly.
library(broom)
# Fit a logistic regression: predict heavy penguin (above median mass)
model_data <- penguins |>
drop_na() |>
mutate(heavy = as.integer(body_mass_g > median(body_mass_g)))
fit <- glm(heavy ~ bill_length_mm + bill_depth_mm + flipper_length_mm +
species + sex,
data = model_data, family = binomial)
# Tidy the model output (use Wald CIs for stability)
model_tidy <- tidy(fit, exponentiate = TRUE) |>
filter(term != "(Intercept)") |>
mutate(
conf.low = exp(log(estimate) - 1.96 * std.error),
conf.high = exp(log(estimate) + 1.96 * std.error),
term = case_match(term,
"bill_length_mm" ~ "Bill length (mm)",
"bill_depth_mm" ~ "Bill depth (mm)",
"flipper_length_mm" ~ "Flipper length (mm)",
"speciesChinstrap" ~ "Chinstrap vs. Adelie",
"speciesGentoo" ~ "Gentoo vs. Adelie",
"sexmale" ~ "Sex (male vs. female)"
),
significant = p.value < 0.05
)
ggplot(model_tidy, aes(x = estimate, y = fct_reorder(term, estimate),
color = significant)) +
geom_vline(xintercept = 1, linetype = "dashed", color = "grey50") +
geom_pointrange(aes(xmin = conf.low, xmax = conf.high),
size = 0.8, linewidth = 1.1) +
geom_text(aes(label = sprintf("OR = %.2f", estimate)),
vjust = -1, size = 3.8, show.legend = FALSE) +
scale_color_manual(values = c("TRUE" = "#e64173", "FALSE" = "grey60"),
guide = "none") +
scale_x_log10() +
labs(title = "Predictors of Above-Median Body Mass (Logistic Regression)",
subtitle = "Odds ratios with 95% CI on log scale; red = p < 0.05",
x = "Odds Ratio (log scale)", y = NULL,
caption = "Model: glm(heavy ~ bill + flipper + species + sex, family = binomial)")Figure 27
Show what your model predicts across the range of a key variable, holding others at their means.
library(broom)
# Generate predictions across flipper length range
pred_grid <- tibble(
flipper_length_mm = seq(170, 235, length.out = 200),
bill_length_mm = mean(model_data$bill_length_mm),
bill_depth_mm = mean(model_data$bill_depth_mm),
species = "Adelie",
sex = "female"
)
preds <- augment(fit, newdata = pred_grid, type.predict = "response",
se_fit = TRUE) |>
mutate(lower = pmax(.fitted - 1.96 * .se.fit, 0),
upper = pmin(.fitted + 1.96 * .se.fit, 1))
ggplot(preds, aes(x = flipper_length_mm, y = .fitted)) +
geom_ribbon(aes(ymin = lower, ymax = upper),
fill = "#2E86AB", alpha = 0.2) +
geom_line(color = "#2E86AB", linewidth = 1.3) +
geom_rug(data = model_data,
aes(x = flipper_length_mm, y = heavy),
sides = "tb", alpha = 0.15, color = "grey40") +
scale_y_continuous(labels = label_percent()) +
labs(title = "Predicted Probability of Above-Median Mass by Flipper Length",
subtitle = "Logistic regression (Adelie, female); rug marks show observed data",
x = "Flipper length (mm)", y = "Predicted probability",
caption = "Shaded band = approximate 95% confidence interval")Figure 28
Predicted probability curves are one of the most effective ways to communicate logistic regression results to non-statisticians.
The ggeffects package automates predicted value plots for any model class.
library(ggeffects)
# One line to get predicted values
preds <- ggpredict(fit,
terms = c("flipper_length_mm",
"sex"))
# Built-in plot method
plot(preds) +
labs(
title = "Marginal Effect of Flipper Length by Sex"
)
# Or extract data for custom ggplot
as.data.frame(preds) |>
ggplot(aes(x, predicted,
color = group)) +
geom_ribbon(...) +
geom_line(...)
performance + see PackagesThe performance package (from easystats) provides model diagnostics; see visualizes them with ggplot2.
# Fit a linear model for diagnostics demo
lm_fit <- lm(body_mass_g ~ bill_length_mm + bill_depth_mm + flipper_length_mm + species,
data = drop_na(penguins))
library(performance)
library(see)
check_model(lm_fit, check = c("linearity", "normality", "qq", "homogeneity"))Figure 30
check_model() replaces the base R plot(model) with a modern, multi-panel diagnostic dashboard. One function call, publication-ready output.
When comparing nested or competing models, visualize the fit statistics side by side.
# Fit competing models
m1 <- glm(heavy ~ flipper_length_mm, data = model_data, family = binomial)
m2 <- glm(heavy ~ flipper_length_mm + species, data = model_data, family = binomial)
m3 <- glm(heavy ~ flipper_length_mm + species + sex, data = model_data, family = binomial)
m4 <- fit # full model from earlier
# Helper to get Wald CIs
tidy_wald <- function(mod, label) {
tidy(mod, exponentiate = TRUE) |>
mutate(conf.low = exp(log(estimate) - 1.96 * std.error),
conf.high = exp(log(estimate) + 1.96 * std.error),
model = label)
}
# Compare coefficients across models
models_tidy <- bind_rows(
tidy_wald(m1, "Model 1:\nFlipper only"),
tidy_wald(m2, "Model 2:\n+ Species"),
tidy_wald(m3, "Model 3:\n+ Sex"),
tidy_wald(m4, "Model 4:\n+ Bill measures")
) |>
filter(term == "flipper_length_mm")
ggplot(models_tidy, aes(x = estimate, y = model)) +
geom_vline(xintercept = 1, linetype = "dashed", color = "grey50") +
geom_pointrange(aes(xmin = conf.low, xmax = conf.high),
color = "#2E86AB", size = 1, linewidth = 1.1) +
geom_text(aes(label = sprintf("OR = %.2f (%.2f, %.2f)",
estimate, conf.low, conf.high)),
vjust = -1.2, size = 3.8, color = "grey30") +
scale_x_log10() +
labs(title = "How Stable Is the Flipper Length Effect Across Models?",
subtitle = "Odds ratio for flipper_length_mm as covariates are added",
x = "Odds Ratio (log scale)", y = NULL,
caption = "Stable estimates across nested models suggest robust association")Figure 31
A waffle chart is a part-to-whole alternative to pie charts, popularized by the #30DayChartChallenge. Each square = 1 unit.
library(waffle)
penguin_counts <- penguins |>
drop_na() |>
count(species) |>
mutate(n_scaled = round(n / 5)) # each square = 5 penguins
waffle(
c("Adelie" = penguin_counts$n_scaled[1],
"Chinstrap" = penguin_counts$n_scaled[2],
"Gentoo" = penguin_counts$n_scaled[3]),
rows = 5,
size = 1,
colors = c("#FF6B35", "#A23B72", "#2E86AB"),
title = "Palmer Penguins by Species",
xlab = "1 square = 5 penguins"
) +
theme(plot.title = element_text(face = "bold", size = 16),
legend.position = "bottom")Figure 32
Challenge for you: pick a prompt from the #30DayChartChallenge and create a chart using a dataset from this semester. Share it!
