--- title: "Specification Curve Analysis with speccurvieR" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Specification Curve Analysis with speccurvieR} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5, message = FALSE, warning = FALSE ) library(speccurvieR) ``` ## What is specification curve analysis? When you fit a regression you make many choices: which controls to include, which fixed effects, which standard errors. Each choice is defensible, and each can move your estimate. **Specification curve analysis** (Simonsohn, Simmons, and Nelson 2020) makes that garden of forking paths explicit: it estimates the model under *every* reasonable combination of choices, plots the resulting curve of estimates, and lets you ask whether your result is robust to the choices you *didn't* make. `speccurvieR` makes this easy, fast, and pretty. This vignette walks through the whole workflow on the bundled `bottles` data — a sample of the [CalCOFI bottle database](https://calcofi.org/data/oceanographic-data/bottle-database/) of oceanographic measurements. ```{r} data(bottles) names(bottles) ``` ## Estimating the curve `sca()` is the workhorse. Give it a dependent variable `y`, a focal independent variable `x`, and a set of `controls`; it estimates the model with every combination of those controls and returns a tidy data frame — one row per specification. ```{r} s <- sca(y = "Salnty", x = "T_degC", controls = c("O2Sat", "ChlorA", "NO2uM", "SiO3uM"), data = bottles, progress_bar = FALSE) dim(s) s[1:3, c("coef", "se", "p", "sig.level")] ``` With four controls there are `r nrow(s)` specifications. Each row carries the focal coefficient, its standard error and p-value, model-fit measures, and indicator columns recording which controls were included. If you prefer, you can specify the whole model as a formula. The first right-hand-side term is the focal variable, the rest are controls, and anything after a `|` is treated as fixed effects: ```{r} s_formula <- sca(Salnty ~ T_degC + O2Sat + ChlorA + NO2uM + SiO3uM, data = bottles, progress_bar = FALSE) ``` `sca()` also supports weights, `glm()` families (via `family`/`link`), fixed effects through `fixest::feols()`, and parallel estimation (`parallel = TRUE`) for very large curves. ## Visualizing the curve `plot_curve()` draws the specification curve: each specification's estimate, ranked, coloured by significance, with a panel below showing which controls are active in each model. ```{r} plot_curve(s, title = "Effect of temperature on salinity") ``` You can request just the curve (`plot_vars = FALSE`) and, because every plot is a `ggplot2` object, customize it freely. Model fit is easy to inspect across specifications, as are the distributions of the control coefficients: ```{r} plot_rmse(s) plot_control_distributions(s) ``` ## Diagnostic plots `speccurvieR` adds a few diagnostics that other packages do not. `plot_influence()` shows, for each control, how *including* versus *excluding* it shifts the focal coefficient — making clear which modelling choices actually move your estimate: ```{r} plot_influence(s) ``` `plot_coef_fit()` plots the focal coefficient against model fit, so you can see whether your best-fitting specifications give systematically different estimates: ```{r} plot_coef_fit(s) ``` ## Comparing standard errors A specification curve compares *point estimates*; `se_compare()` complements it by comparing *standard errors*. Pass a formula and the standard-error types you want — heteroskedasticity-consistent, clustered, bootstrapped, or all of them: ```{r} se_compare(formula = "Salnty ~ T_degC + ChlorA", data = bottles, types = c("iid", "HC0", "HC3")) ``` `se_compare()` mirrors `sca()`'s model interface, so it also accepts `glm()` families. `plot_se()` visualizes the result, showing each coefficient's estimate with a confidence interval under every standard-error type. ## Joint-inference test Looking at a curve tells you whether results are robust *descriptively*. But how do you know the curve as a whole is more extreme than you would expect by chance? `sca_test()` answers this with the permutation-based joint-inference test of Simonsohn, Simmons, and Nelson (2020). It tests the sharp null that the focal variable has **no effect in any specification**. To do so it repeatedly shuffles the focal variable — breaking its association with the outcome — re-estimates the *entire* curve each time, and compares the observed curve to the resulting null distribution. (When fixed effects are present the shuffle is blocked within the first fixed effect, because the focal variable is only exchangeable within those fixed-effect groups.) Three statistics summarize the curve, each with its own permutation p-value: the median estimate, the share of specifications that are statistically significant (restricted to the predicted direction when you set `direction`), and a Stouffer combination of the per-specification p-values. ```{r} result <- sca_test(y = "Salnty", x = "T_degC", controls = c("O2Sat", "ChlorA", "NO2uM"), data = bottles, n_permutations = 199, seed = 1, keep_curves = TRUE, progress_bar = FALSE) result ``` `plot_sca_test()` shows each statistic's null distribution with the observed value marked, so you can see how far into the tail the real curve sits: ```{r, fig.height = 3} plot_sca_test(result) ``` And `plot_sca_test_specs()` — requested via `keep_curves = TRUE` above — gives the specification-curve view: each specification's observed estimate against its *own* null band, with specifications outside their band highlighted. ```{r} plot_sca_test_specs(result) ``` ### A few caveats - The test is two-sided by default. Use `direction = "positive"` or `"negative"` only when you have an a-priori predicted direction; the package never chooses the direction from the data. - Permutation p-values have a resolution floor of `1 / (n_permutations + 1)` (more precisely, of the permutations that estimate successfully). Use a few hundred permutations for exploration and `n_permutations >= 1000` to resolve small p-values; `parallel = TRUE` spreads the work across workers. - The default null shuffles the focal variable, which also breaks its correlation with the controls and is miscalibrated under collinearity. For observational data, `null_type = "freedman_lane"` (permute the residuals of a reduced `y ~ controls` model) or `null_type = "residual_bootstrap"` preserve that correlation; both are for linear models and fit every specification on a common sample. ### Which specifications are real? The joint test asks whether the curve *as a whole* is more extreme than chance. A natural next question is *which individual* specifications are real. You cannot just read the per-specification p-values off the curve and call the significant ones robust: with dozens of correlated specifications, some are bound to look significant by chance — that is the multiple-comparisons problem all over again, now *across* specifications. `speccurvieR` corrects for it with the min-P / max-statistic permutation method of Westfall and Young (1993). It reuses the permutations the joint test already ran (so there is nothing extra to compute) and, for each one, records the most extreme specification anywhere in the curve — the null distribution of "the best result a search could turn up by chance" — then compares each specification to it. When you run `sca_test()` with `keep_curves = TRUE` **and** a confound-preserving null, these family-wise-error-rate-adjusted p-values are attached automatically: ```{r} result_fwer <- sca_test(y = "Salnty", x = "T_degC", controls = c("O2Sat", "ChlorA", "NO2uM"), data = bottles, null_type = "freedman_lane", n_permutations = 499, seed = 1, keep_curves = TRUE, progress_bar = FALSE) as.data.frame(result_fwer, what = "specs") ``` `p_raw` is each specification's uncorrected p-value against its *own* null; `p_adj` is the corrected one; `significant_adj` flags those that survive. The specification curve plot now highlights the survivors in a third colour: ```{r} plot_sca_test_specs(result_fwer) ``` The default correction (`"single_step"`) controls the family-wise error rate in the weak sense. `sca_minp()` recomputes it — without re-running the permutations — with Westfall and Young's more powerful step-down refinement, or a different threshold: ```{r} result_fwer <- sca_minp(result_fwer, method = "step_down") result_fwer ``` A caution on interpretation: each specification is compared to its *own* null, which under a confound-preserving null is centred on the value that specification produces when the focal variable has no partial effect. A flagged *under-controlled* specification therefore signals an association beyond that conditional null — it is **not** evidence of a causal effect of the size it reports, because omitting a confounder re-routes part of the focal coefficient. ## Reporting and export speccurvieR objects come with `tidy()` and `glance()` methods (the same generics `broom` and `modelsummary` use), so results drop straight into a reporting pipeline: ```{r} tidy(result) ``` `sca_table()` produces a publication-ready results block. By default it returns a small, dependency-free data frame, but it can also render to Markdown, LaTeX, or (if the package is installed) `gt`, `kableExtra`, or `flextable`: ```{r} sca_table(result) ``` And `sca_report()` writes a one-paragraph, manuscript-ready summary: ```{r} sca_report(result) ``` Both `sca_table()` and `sca_report()` also accept an `sca()` curve, in which case they give a descriptive (non-inferential) summary. ## References Simonsohn, U., Simmons, J. P., & Nelson, L. D. (2020). Specification curve analysis. *Nature Human Behaviour*, 4, 1208–1214.