---
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.
<https://doi.org/10.1038/s41562-020-0912-z>
