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 of oceanographic measurements.
data(bottles)
names(bottles)
#> [1] "Cst_Cnt" "Btl_Cnt" "Sta_ID"
#> [4] "Depth_ID" "Depthm" "T_degC"
#> [7] "Salnty" "O2ml_L" "STheta"
#> [10] "O2Sat" "Oxy_µmol.Kg" "BtlNum"
#> [13] "RecInd" "T_prec" "T_qual"
#> [16] "S_prec" "S_qual" "P_qual"
#> [19] "O_qual" "SThtaq" "O2Satq"
#> [22] "ChlorA" "Chlqua" "Phaeop"
#> [25] "Phaqua" "PO4uM" "PO4q"
#> [28] "SiO3uM" "SiO3qu" "NO2uM"
#> [31] "NO2q" "NO3uM" "NO3q"
#> [34] "NH3uM" "NH3q" "C14As1"
#> [37] "C14A1p" "C14A1q" "C14As2"
#> [40] "C14A2p" "C14A2q" "DarkAs"
#> [43] "DarkAp" "DarkAq" "MeanAs"
#> [46] "MeanAp" "MeanAq" "IncTim"
#> [49] "LightP" "R_Depth" "R_TEMP"
#> [52] "R_Sal" "R_DYNHT" "R_Nuts"
#> [55] "R_Oxy_µmol.Kg" "DIC1" "DIC2"
#> [58] "TA1" "TA2" "pH1"
#> [61] "pH2" "DIC.Quality.Comment"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.
s <- sca(y = "Salnty", x = "T_degC",
controls = c("O2Sat", "ChlorA", "NO2uM", "SiO3uM"),
data = bottles, progress_bar = FALSE)
dim(s)
#> [1] 15 15
s[1:3, c("coef", "se", "p", "sig.level")]
#> coef se p sig.level
#> Salnty ~ T_degC + NO2uM -0.06985492 0.005696327 4.988171e-26 p < .005
#> Salnty ~ T_degC + ChlorA + NO2uM -0.06210476 0.007568722 3.912279e-13 p < .005
#> Salnty ~ T_degC + ChlorA -0.05997833 0.007428642 6.504071e-13 p < .005With four controls there are 15 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:
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.
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.
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:
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:
plot_coef_fit() plots the focal coefficient against
model fit, so you can see whether your best-fitting specifications give
systematically different estimates:
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:
se_compare(formula = "Salnty ~ T_degC + ChlorA", data = bottles,
types = c("iid", "HC0", "HC3"))
#> estimate iid HC0 HC3
#> (Intercept) 34.2940251811 0.097594017 0.107876468 0.111322464
#> T_degC -0.0599783335 0.007428642 0.008370367 0.008669478
#> ChlorA 0.0006514447 0.012449618 0.005327664 0.011852569se_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.
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.
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
#> Specification curve joint-inference test (Simonsohn, Simmons & Nelson 2020)
#>
#> Focal variable: T_degC
#> Specifications: 7
#> Permutations: 199 used (0 failed) | blocked within FE: no
#> Null: shuffle x
#> Direction: two.sided alpha = 0.05
#>
#> Statistic Observed p-value
#> Median estimate 0.0543 0.0050
#> Share significant 1.0000 0.0050
#> Stouffer Z 0.9213 0.6450
#>
#> p-values are permutation-based; resolution floor = 0.0050.
#> Interpretation (SSN): conclude a robust effect when the median test AND
#> at least one of {share significant, Stouffer} are significant.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:
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.
direction = "positive" or "negative" only when
you have an a-priori predicted direction; the package never chooses the
direction from the data.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.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.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:
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")
#> spec observed p_raw p_adj significant_adj
#> 1 ChlorA + NO2uM -0.06210476 0.018 0.050 FALSE
#> 2 NO2uM -0.06199380 0.020 0.060 FALSE
#> 3 ChlorA -0.06167098 0.016 0.046 TRUE
#> 4 NO2uM + O2Sat 0.03066608 0.002 0.002 TRUE
#> 5 O2Sat 0.03141611 0.002 0.002 TRUE
#> 6 ChlorA + NO2uM + O2Sat 0.05577464 0.002 0.002 TRUE
#> 7 ChlorA + O2Sat 0.05610969 0.002 0.002 TRUEp_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:
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:
result_fwer <- sca_minp(result_fwer, method = "step_down")
result_fwer
#> Specification curve joint-inference test (Simonsohn, Simmons & Nelson 2020)
#>
#> Focal variable: T_degC
#> Specifications: 7
#> Permutations: 499 used (0 failed) | blocked within FE: no
#> Null: Freedman-Lane (control superset, common sample)
#> Direction: two.sided alpha = 0.05
#>
#> Statistic Observed p-value
#> Median estimate 0.0307 0.0800
#> Share significant 1.0000 0.0100
#> Stouffer Z -1.8034 1.0000
#>
#> p-values are permutation-based; resolution floor = 0.0020.
#>
#> After correcting for searching 7 specifications, 7 remain statistically significant
#> (smallest corrected p = 0.0020).
#> Interpretation (SSN): conclude a robust effect when the median test AND
#> at least one of {share significant, Stouffer} are significant.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.
speccurvieR objects come with tidy() and
glance() methods (the same generics broom and
modelsummary use), so results drop straight into a
reporting pipeline:
tidy(result)
#> term estimate p.value direction null_type
#> 1 median 0.05428491 0.005 two.sided shuffle_x
#> 2 share_significant 1.00000000 0.005 two.sided shuffle_x
#> 3 stouffer 0.92130513 0.645 two.sided shuffle_xsca_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:
sca_table(result)
#> Focal variable T_degC
#> Specifications 7
#> Observations per specification
#> Permutations 199
#> Null hypothesis shuffle x
#> Direction two.sided
#> alpha 0.05
#> Median estimate 0.0543 (p = 0.0050)
#> Share significant 100.0% (p = 0.0050)
#> Stouffer Z 0.921 (p = 0.6450)
#>
#> p-values are permutation-based; resolution floor = 0.0050.And sca_report() writes a one-paragraph,
manuscript-ready summary:
sca_report(result)
#> [1] "A specification curve analysis estimated the effect of T_degC across 7 specifications. The median estimate was 0.0543 (100.0% of specifications were statistically significant at the 0.05 level). Joint inference via shuffle-x permutation tests (199 permutations) rejected the null of no effect: the median estimate (p = 0.0050), the share of significant specifications (p = 0.0050), and Stouffer's combined test (Z = 0.921, p = 0.6450)."Both sca_table() and sca_report() also
accept an sca() curve, in which case they give a
descriptive (non-inferential) summary.
Simonsohn, U., Simmons, J. P., & Nelson, L. D. (2020). Specification curve analysis. Nature Human Behaviour, 4, 1208–1214. https://doi.org/10.1038/s41562-020-0912-z