How to plot a specification curve

Updates

Updated 2025-11-24 to account for breaking package changes.

Like many researchers, I often want to plot a range of coefficient estimates to figure out if the results I’m finding are robust to other sensible specification and functional form choices. This kind of estimate is called a specification curve (Simonsohn, Simmons, and Nelson 2015), and I am far from the first to do it. In fact, there are even a couple packages available: Joachim Gassen’s rdfanalysis and Philipp Masur’s specr (I haven’t used either, yet).

I wanted to roll my own, though. Since my code is fairly simple (less than 60 lines including comments and generous spacing) and uses the tidyverse, it may be helpful to other people too.

pacman::p_load(tidyverse, cowplot, fastDummies)
theme_set(theme_cowplot())
set.seed(42)

First, we create some fake estimates. This includes the point estimates, standard errors, and any variables that describe the specifications we “estimated”.

fe <- c("Unit", "Time", "Unit + Time")
controls <- c("Basic", "Some", "Full")

estimates <- as_tibble(expand.grid(`Fixed Effects` = fe, `Controls` = controls))
estimates <- estimates %>% mutate(est = rnorm(n()), se = runif(n(), 0, 0.1))
estimates <- estimates %>% mutate(ci_l = est - 1.96 * se, ci_h = est + 1.96 * se)

Next, we construct a typical coefficient plot that shows how the estimates and confidence intervals change across different categories.

spec_cols <- c("Fixed Effects", "Controls") # Could be set by user or figured out

estimates <- estimates %>%
  arrange(est) %>%
  mutate(h_order = 1:n()) # Sort on point estimates for horizontal ordering

coef_plot <- ggplot(estimates, aes(x = h_order, y = est)) +
  geom_linerange(aes(ymin = ci_l, ymax = ci_h), linewidth = 1, alpha = 0.5) +
  geom_point(fill = "white", shape = 21) +
  labs(y = "Coefficient") +
  theme(axis.title.x = element_blank(), axis.ticks.x = element_blank(), axis.line.x = element_blank(), axis.text.x = element_blank())

Then we visualize the specifications themselves as a series of dot plots below the coefficient plot.

make_spec_plot <- function(category) {
  # Creates a dot plot for one specification category
  specs <- dummy_cols(estimates, select_columns = category, remove_selected_columns = T) %>%
    select(h_order, starts_with(category)) %>%
    pivot_longer(starts_with(category), names_prefix = paste0(category, "_")) %>%
    mutate(name = factor(name, levels = rev(unique(name))))

  spec_plot <- ggplot(specs, aes(x = h_order, y = name, alpha = value)) +
    geom_point() +
    scale_alpha_continuous(guide = "none") +
    theme(axis.title.x = element_blank(), axis.ticks.x = element_blank(), axis.line.x = element_blank(), axis.text.x = element_blank()) +
    theme(axis.text.y = element_text(size = 6), axis.title.y = element_blank(), axis.ticks.y = element_blank(), axis.line.y = element_blank())
  spec_plot
}
spec_plots <- lapply(spec_cols, make_spec_plot)

Finally, we combine them for the final result:

combined_plot <- plot_grid(
  plotlist = c(list(coef_plot), spec_plots),
  labels = c("", spec_cols), label_size = 8, label_fontface = "italic", vjust = 0.5, hjust = -0.1,
  rel_heights = c(1, 0.2, 0.2), ncol = 1, align = "v"
)

combined_plot

This approach is flexible and can be easily adapted to include more specification categories, different aesthetics, or other customizations. Happy plotting!