Package {speccurvieR}

speccurvieR: Easy, Fast, and Pretty Specification Curve Analysis

Description

Making specification curve analysis easy, fast, and pretty. It improves upon existing offerings with additional features and 'tidyverse' integration. Users can easily visualize and evaluate how their models behave under different specifications with a high degree of customization. For a description and applications of specification curve analysis see Simonsohn, Simmons, and Nelson (2020) doi:10.1038/s41562-020-0912-z.

Author(s)

Maintainer: Zayne Sember zayne@mit.edu [copyright holder]

See Also

Useful links:

• https://github.com/zaynesember/speccurvieR

• Report bugs at https://github.com/zaynesember/speccurvieR/issues

CalCOFI Bottle Data

Description

A subset of data from the California Cooperative Oceanic Fisheries Investigations. Each observation describes a sample of ocean water collected.

Usage

bottles

Format

## 'bottles' A data frame with 500 rows and 62 columns:

Cst_Cnt

Cast count

Btl_Cnt

Bottle Count

Sta_ID

Line and Station

Depth_ID

Depth ID

Depthm

Bottle depth in meters

T_degC

Water temperature in degrees Celsius

Salnty

Salinity (Practical Salinity Scale 1978)

O2ml_L

Milliliters of oxygen per liter of seawater

STheta

Potential Density (Sigma Theta), Kg/M³

O2Sat

Oxygen percent saturation

Oxy_µmol/Kg

Oxygen micromoles per kilogram seawater

BtlNum

Niskin bottle sample was collected from

RecInd

Record Indicator

T_prec

Temperature Precision

T_qual

Quality Code

S_prec

Salinity Precision

S_qual

Quality Code

P_qual

Quality Code

O_qual

Quality Code

SThtaq

Quality Code

O2Satq

Quality Code

ChlorA

Micrograms Chlorophyll-a per liter seawater

Chlqua

Quality Code

Phaeop

Micrograms Phaeopigment per liter seawater

Phaqua

Quality Code

PO4uM

Micromoles Phosphate per liter of seawater

PO4q

Quality Code

SiO3uM

Micromoles Silicate per liter of seawater

SiO3qu

Quality Code

NO2uM

Micromoles Nitrite per liter of seawater

NO2q

Quality Code

NO3uM

Micromoles Nitrate per liter of seawater

NO3q

Quality Code

NH3uM

Micromoles Ammonia per liter of seawater

NH3q

Quality Code

C14As1

14C Assimilation of Replicate 1

C14A1p

Precision of 14C Assimilation of Replicate 1

C14A1q

Quality Code

C14As2

14C Assimilation of Replicate 2

C14A2p

Precision of 14C Assimilation of Replicate 2

C14A2q

Quality Code

DarkAs

14C Assimilation of Dark/Control Bottle

DarkAp

Precision of 14C Assimilationof Dark/Control Bottle

Darkaq

Quality Code

MeanAs

Mean 14C Assimilation of Replicates 1 and 2

MeanAp

Precision of Mean 14C Assimilation of Replicates 1 and 2

MeanAq

Quality Code

IncTim

Elapsed incubation time of primary productivity experiment

LightP

Light intensities of the incubation tubes

R_Depth

Reported Depth (from pressure) in meters

R_Temp

Reported (Potential) Temperature in degrees Celsius

R_Sal

Reported Salinity (from Specific Volume Anomoly, M³/Kg)

R_DYNHT

Reported Dynamic Height in units of dynamic meters

R_Nuts

Reported Ammonium concentration

R_Oxy_µmol/Kg

Reported Oxygen micromoles/kilogram

DIC1

Dissolved Inorganic Carbon micromoles per kilogram solution

DIC2

Dissolved Inorganic Carbon on a replicate sample

TA1

Total Alkalinity micromoles per kilogram solution

TA2

Total Alkalinity on a replicate sample

pH1

pH (the degree of acidity/alkalinity of a solution)

pH2

pH on a replicate sample

DIC Quality Comment

Quality Comment

Source

<https://calcofi.org/data/oceanographic-data/bottle-database/>

Extracts the control variable names and coefficients from an lm model summary.

Description

Extracts the control variable names and coefficients from a model summary.

Usage

control_extractor(model, x, feols_model = FALSE)

Arguments

Value

A dataframe with two columns, 'term' contains the name of the control and 'coef' contains the coefficient estimate.

Examples

m <- summary(lm(Salnty ~ STheta + T_degC, bottles))
control_extractor(model = m, x = "STheta");

m <- summary(lm(Salnty ~ STheta*T_degC + O2Sat, bottles))
control_extractor(model = m, x = "STheta");

Removes duplicate control variables

Description

Removes duplicate control variables from user input.

Usage

duplicate_remover(controls, x)

Arguments

Value

A vector of strings containing control variable names

Examples

duplicate_remover(controls = c("control1", "control2*control3"),
                  x = "independentVariable");

Builds models formulae with every combination of control variables possible.

Description

Builds models formulae with every combination of control variables possible.

Usage

formula_builder(y, x, controls, fixed_effects = NA, ...)

Arguments

Value

A vector of formula objects using every possible combination of controls.

Examples

formula_builder("dependentVariable", "independentVariable",
                c("control1", "control2"));
formula_builder("dependentVariable", "independentVariable",
                c("control1*control2"), fixed_effects="month");

Paste together controls and independent variable

Description

'paste_factory()' constructs the right hand side of the regression as a a string i.e. "x + control1 + control2".

Usage

paste_factory(controls, x)

Arguments

Value

A string concatenating independent and control variables separated by '+'.

Examples

paste_factory(controls = c("control1", "control2"),
              x = "independentVariable");

Plots the AIC across model specifications.

Description

plot_aic() plots the Akaike information criterion across model specifications. Only available for nonlinear regression models.

Usage

plot_aic(sca_data, title = "", show_index = TRUE, plot_vars = TRUE)

Arguments

Value

If 'plot_vars = TRUE' a 'patchwork' object combining the plot and the variable panel; if 'plot_vars = FALSE' a ggplot object.

Examples

plot_aic(sca_data = sca(y = "Salnty", x = "T_degC",
                       controls = c("ChlorA", "O2Sat"),
                       data = bottles, progress_bar = TRUE, parallel = FALSE),
                     title = "AIC");
plot_aic(sca_data = sca(y = "Salnty", x = "T_degC",
                       controls = c("ChlorA*O2Sat"),
                       data = bottles, progress_bar = FALSE,
                       parallel = FALSE),
                       show_index = FALSE, plot_vars = FALSE);

plot_aic(sca_data = sca(y = "Salnty", x = "T_degC",
                         controls = c("ChlorA*NO3uM", "O2Sat*NO3uM"),
                         data = bottles,
                         progress_bar = TRUE, parallel = TRUE, workers = 2));

Plots the independent variable's coefficient against model fit.

Description

plot_coef_fit() plots the independent variable's coefficient against a measure of model fit across specifications, revealing whether better-fitting models tend to produce systematically different estimates (i.e. whether your best-fitting specifications are outliers).

Usage

plot_coef_fit(sca_data, metric = NULL, title = "")

Arguments

Value

A ggplot object.

Examples

plot_coef_fit(sca(y = "Salnty", x = "T_degC",
                controls = c("ChlorA", "O2Sat", "NO2uM"),
                data = bottles, progress_bar = FALSE));

Plots control variable distributions.

Description

plot_control_distributions() plots the distribution of coefficients for each control variable included in the model specifications.

Usage

plot_control_distributions(
  sca_data,
  title = "",
  type = "density",
  zero_line = TRUE
)

Arguments

Value

A ggplot object.

Examples

plot_control_distributions(sca_data = sca(y="Salnty", x="T_degC",
                                    controls = c("ChlorA", "O2Sat"),
                                    data = bottles,
                                    progress_bar = TRUE, parallel = FALSE),
                         title = "Control Variable Distributions")
plot_control_distributions(sca_data = sca(y = "Salnty", x="T_degC",
                                    controls = c("ChlorA*O2Sat"),
                                    data = bottles,
                                    progress_bar = FALSE, parallel = FALSE),
                         type = "hist")

plot_control_distributions(sca_data = sca(y = "Salnty", x = "T_degC",
                                    controls = c("ChlorA*NO3uM",
                                                 "O2Sat*NO3uM"),
                                    data = bottles, progress_bar = TRUE,
                                    parallel = TRUE, workers = 2),
                         type = "density")

Plots a specification curve.

Description

plot_curve() takes the data frame output of sca() and produces a ggplot of the independent variable's coefficient (as indicated in the call to sca()) across model specifications. By default a panel is added showing which control variables are present in each model. The combined plot is returned as a 'patchwork' object, so it can be further customised with ggplot2 and patchwork operators (e.g. '& theme_sca(base_size = 14)').

Usage

plot_curve(
  sca_data,
  title = "",
  show_index = TRUE,
  plot_vars = TRUE,
  ylab = "Coefficient",
  plot_se = "bar",
  median_line = FALSE,
  point_size = NULL
)

Arguments

Value

If 'plot_vars = TRUE' a 'patchwork' object combining the curve and the variable panel; if 'plot_vars = FALSE' a ggplot object. Both can be further customised with ggplot2 (and patchwork) operators.

Examples

plot_curve(sca_data = sca(y="Salnty", x="T_degC", c("ChlorA", "O2Sat"),
                         data=bottles, progress_bar=TRUE, parallel=FALSE),
                     title = "Salinity and Temperature Models",
                     show_index = TRUE, plot_vars = TRUE,
                     ylab = "Coefficient value", plot_se = "ribbon");
plot_curve(sca_data = sca(y="Salnty", x="T_degC",
                         c("ChlorA*O2Sat", "ChlorA", "O2Sat"),
                         data=bottles, progress_bar=FALSE, parallel=FALSE),
                     show_index = TRUE, plot_vars = TRUE,
                     plot_se = "ribbon");

plot_curve(sca_data = sca(y="Salnty", x="T_degC",
                         c("ChlorA*NO3uM", "O2Sat", "ChlorA", "NO3uM"),
                         data=bottles,
                         progress_bar = TRUE, parallel = TRUE, workers=2),
          plot_se="");

Plots the deviance of residuals across model specifications.

Description

plot_deviance() plots the deviance of residuals across model specifications. Only available for linear regression models.

Usage

plot_deviance(sca_data, title = "", show_index = TRUE, plot_vars = TRUE)

Arguments

Value

If 'plot_vars = TRUE' a 'patchwork' object combining the plot and the variable panel; if 'plot_vars = FALSE' a ggplot object.

Examples

plot_deviance(sca_data = sca(y = "Salnty", x = "T_degC",
                            controls = c("ChlorA", "O2Sat"),
                            data = bottles, progress_bar = TRUE,
                            parallel = FALSE),
                     title = "Model Deviance");
plot_deviance(sca_data = sca(y = "Salnty", x = "T_degC",
                            controls = c("ChlorA*O2Sat"),
                            data = bottles, progress_bar = FALSE,
                            parallel = FALSE),
                     show_index = FALSE, plot_vars = FALSE);

plot_deviance(sca_data = sca(y = "Salnty", x="T_degC",
                         controls = c("ChlorA*NO3uM", "O2Sat*NO3uM"),
                         data = bottles, progress_bar = TRUE, parallel = TRUE,
                         workers = 2));

Plots how each control influences the independent variable's coefficient.

Description

plot_influence() shows, for every control variable, the distribution of the independent variable's coefficient across the specifications that include versus exclude that control. It makes clear which modelling choices move the estimate, and by how much.

Usage

plot_influence(sca_data, title = "")

Arguments

Value

A ggplot object with one facet per control comparing the coefficient when that control is excluded versus included.

Examples

plot_influence(sca(y = "Salnty", x = "T_degC",
                  controls = c("ChlorA", "O2Sat", "NO2uM"),
                  data = bottles, progress_bar = FALSE));

Plots a specification curve under multiple standard error types.

Description

plot_multi_se() estimates every specification (as 'sca()' does) and, for each, computes the independent variable's standard error under several types via 'se_compare()'. It then plots the specification curve faceted by standard error type: the coefficient estimates are identical across facets, but the confidence intervals – and hence which specifications are "significant" – change with the choice of standard error, showing how sensitive your conclusions are to that choice.

Usage

plot_multi_se(
  y,
  x,
  controls,
  data,
  types = c("iid", "HC3"),
  cluster = NULL,
  fixed_effects = NULL,
  level = 0.95,
  title = ""
)

Arguments

Value

A ggplot object: the specification curve faceted by standard error type, points coloured by significance under each type.

Examples

plot_multi_se(y = "Salnty", x = "T_degC", controls = c("ChlorA", "O2Sat"),
            data = bottles, types = c("iid", "HC1", "HC3"));

Plots the adj. R-squared across model specifications.

Description

plot_r2_adj() plots the adjusted R-squared across model specifications. Only available for linear regression models. Note when fixed effects are are specified the within adjusted R-squared is used (i.e. 'fixest::r2()' with 'type="war2"').

Usage

plot_r2_adj(sca_data, title = "", show_index = TRUE, plot_vars = TRUE)

Arguments

Value

If 'plot_vars = TRUE' a 'patchwork' object combining the plot and the variable panel; if 'plot_vars = FALSE' a ggplot object.

Examples

plot_r2_adj(sca_data = sca(y = "Salnty", x = "T_degC",
                         controls = c("ChlorA", "O2Sat"),
                         data = bottles, progress_bar = TRUE,
                         parallel = FALSE),
                     title = "Adjusted R^2");
plot_r2_adj(sca_data = sca(y="Salnty", x="T_degC",
                         controls = c("ChlorA*O2Sat"),
                         data = bottles, progress_bar = FALSE,
                         parallel = FALSE),
                     show_index = FALSE, plot_vars = FALSE);

plot_r2_adj(sca_data = sca(y = "Salnty", x = "T_degC",
                         controls = c("ChlorA*NO3uM", "O2Sat*NO3uM"),
                         data = bottles,
                         progress_bar = TRUE, parallel = TRUE, workers = 2));

Plots RMSE across model specifications.

Description

plot_rmse() plots the root mean square error across model specifications. Only available for linear regression models.

Usage

plot_rmse(sca_data, title = "", show_index = TRUE, plot_vars = TRUE)

Arguments

Value

If 'plot_vars = TRUE' a 'patchwork' object combining the plot and the variable panel; if 'plot_vars = FALSE' a ggplot object.

Examples

plot_rmse(sca_data = sca(y="Salnty", x="T_degC", c("ChlorA", "O2Sat"),
                         data=bottles, progress_bar=TRUE, parallel=FALSE),
                     title = "RMSE");
plot_rmse(sca_data = sca(y="Salnty", x="T_degC", c("ChlorA*O2Sat"),
                         data=bottles, progress_bar=FALSE, parallel=FALSE),
                     show_index = FALSE, plot_vars = FALSE);

plot_rmse(sca_data = sca(y="Salnty", x="T_degC",
                         c("ChlorA*NO3uM", "O2Sat*NO3uM"), data=bottles,
                         progress_bar = TRUE, parallel=TRUE, workers=2));

Plots the number of observations across model specifications.

Description

plot_samplesizes() plots 'n_obs', the number of observations each specification was fit on, against the specification index. It makes visible whether specifications were fit on different samples – with default listwise deletion they often are – which is a reason to consider 'sca(common_sample = TRUE)'.

Usage

plot_samplesizes(sca_data, title = "")

Arguments

Value

A ggplot object.

See Also

[sca()] and its 'common_sample' argument.

Examples

plot_samplesizes(sca(y = "Salnty", x = "T_degC",
                     controls = c("ChlorA", "O2Sat", "NO2uM"),
                     data = bottles, progress_bar = FALSE, parallel = FALSE))

Plot the null distribution of a joint-inference test

Description

'plot_sca_test()' visualises the output of [sca_test()]. For each test statistic it draws the permutation null distribution with the observed value marked and the permutation p-value annotated, making it easy to see how extreme the observed specification curve is relative to the sharp null.

Usage

plot_sca_test(test_result, type = "histogram", title = "")

Arguments

Value

A ggplot object with one facet per test statistic.

Examples

result <- sca_test(y = "Salnty", x = "T_degC", controls = c("ChlorA", "O2Sat"),
                   data = bottles, n_permutations = 100, progress_bar = FALSE)
plot_sca_test(result)

Plot the specification curve against its per-specification null band

Description

'plot_sca_test_specs()' draws the observed specification curve (each specification's focal coefficient, ranked) together with a shaded band giving that specification's own null distribution under the permutation test. This is the specification-curve view of the joint-inference test of Simonsohn, Simmons, and Nelson (2020): specifications whose observed estimate falls outside their null band are highlighted, so it is easy to see which parts of the curve are more extreme than chance.

When the result carries per-specification family-wise-error-rate-adjusted p-values (i.e. it was computed with 'keep_curves = TRUE' and a confound-preserving null, so [sca_minp()] has run – automatically or explicitly), the points are coloured in three tiers – within the chance band, beyond it but not significant after correction, and significant after multiple-comparison correction (enlarged) – so the specifications that survive correction stand out. Otherwise the usual two-tier within/outside-band colouring is used.

It requires an [sca_test()] result computed with 'keep_curves = TRUE'.

Usage

plot_sca_test_specs(test_result, level = 0.95, title = "")

Arguments

Value

A ggplot object.

References

Simonsohn, U., Simmons, J. P., & Nelson, L. D. (2020). Specification curve analysis. Nature Human Behaviour, 4, 1208-1214. doi:10.1038/s41562-020-0912-z

Examples

result <- sca_test(y = "Salnty", x = "T_degC", controls = c("ChlorA", "O2Sat"),
                   data = bottles, n_permutations = 100, keep_curves = TRUE,
                   progress_bar = FALSE)
plot_sca_test_specs(result)

Plots standard error estimates across types.

Description

plot_se() takes the data frame output of 'se_compare()' and plots, for each coefficient, the estimate together with a confidence interval derived from every available type of standard error. This makes it easy to see how inference about a coefficient changes with the choice of standard error.

Usage

plot_se(se_data, level = 0.95, intercept = FALSE, title = "")

Arguments

Value

A ggplot object with one facet per coefficient; within each facet the estimate is plotted against each standard error type with a confidence interval, and the colour indicates whether that interval excludes zero.

Examples

plot_se(se_compare(formula = "Salnty ~ T_degC + ChlorA", data = bottles,
                  types = c("iid", "HC0", "HC3")));
plot_se(se_compare(formula = "Salnty ~ T_degC + ChlorA", data = bottles,
                  types = "HC1", cluster = c("Sta_ID", "Depth_ID")),
       level = 0.9);

Plot the variance decomposition of a specification curve

Description

'plot_variance()' visualises [sca_variance()] as a bar chart of the share of the focal coefficient's variance attributable to each control choice (plus a '"Residual"' bar), making it easy to see which modelling choices drive the spread of estimates across the curve.

Usage

plot_variance(
  sca_data,
  estimate = "coef",
  method = c("lmg", "type2"),
  residual = TRUE,
  title = ""
)

Arguments

Value

A ggplot object.

See Also

[sca_variance()] for the underlying table.

Examples

s <- sca(y = "Salnty", x = "T_degC",
         controls = c("ChlorA", "O2Sat", "NO2uM"),
         data = bottles, progress_bar = FALSE, parallel = FALSE)
plot_variance(s)

Plots the variables in each model.

Description

plot_vars() plots the variables included in each model specification in order of model index. Returns a ggplot object that can then be combined with the output of other functions like plot_rmse() if further customization of each plot is desired.

Usage

plot_vars(sca_data, title = "", color_controls = FALSE)

Arguments

Value

A ggplot object.

Examples

plot_vars(sca_data = sca(y = "Salnty", x = "T_degC",
                        controls = c("ChlorA", "O2Sat"),
                        data = bottles, progress_bar = TRUE,
                        parallel = FALSE),
                     title = "Model Variable Specifications");
plot_vars(sca_data = sca(y = "Salnty", x = "T_degC",
                        controls = c("ChlorA*O2Sat"),
                        data = bottles, progress_bar = FALSE,
                        parallel = FALSE),
                     color_controls = TRUE);

plot_vars(sca_data = sca(y = "Salnty", x = "T_degC",
                        controls = c("ChlorA*NO3uM", "O2Sat*NO3uM"),
                        data = bottles,
                        progress_bar = TRUE, parallel = TRUE, workers = 2));

Print a specification curve joint-inference test

Description

Print a specification curve joint-inference test

Usage

## S3 method for class 'sca_test'
print(x, ...)

Arguments

Value

'x', invisibly.

Objects exported from other packages

Description

These objects are imported from other packages. Follow the links below to see their documentation.

generics

glance, tidy

Perform specification curve analysis

Description

sca() is the workhorse function of the package–this estimates models with every possible combination of the controls supplied and returns a data frame where each row contains the pertinent information and parameters for a given model by default. This data frame can then be input to plot_curve() or any other plotting function in the package. Alternatively, if 'return_formulae = TRUE', it returns a list of formula objects with every possible combination of controls.

Usage

sca(
  y,
  x,
  controls,
  data,
  weights = NULL,
  family = "linear",
  link = NULL,
  fixed_effects = NULL,
  common_sample = FALSE,
  return_formulae = FALSE,
  progress_bar = TRUE,
  parallel = FALSE,
  workers = 2,
  ...
)

Arguments

Value

When 'return_formulae' is 'FALSE', a dataframe where each row contains the independent variable coefficient estimate, standard error, test statistic, p-value, model specification, measures of model fit, and 'n_obs', the number of observations the specification was fit on.

Examples

sca(y = "Salnty", x = "T_degC", controls = c("ChlorA", "O2Sat"),
    data = bottles, progress_bar = TRUE, parallel = FALSE);
# Equivalent call using the formula interface:
sca(Salnty ~ T_degC + ChlorA + O2Sat, data = bottles, progress_bar = FALSE);
# Formula interface with an interaction control and fixed effects:
sca(Salnty ~ T_degC + ChlorA + ChlorA*O2Sat | Sta_ID, data = bottles,
    progress_bar = FALSE);

sca(y = "Salnty", x = "T_degC", controls = c("ChlorA*NO3uM", "O2Sat*NO3uM"),
    data = bottles, progress_bar = TRUE, parallel = TRUE, workers = 2);

sca(y = "Salnty", x = "T_degC", controls = c("ChlorA", "O2Sat*NO3uM"),
    data = bottles, progress_bar = TRUE, parallel = FALSE,
    return_formulae = TRUE);

Family-wise-error-rate-adjusted p-values for each specification

Description

'sca_minp()' adds, to an existing [sca_test()] result, a multiple-comparison- corrected p-value for *every* specification in the curve. It answers the question the per-specification null band only shows descriptively: *which* individual specifications are more extreme than chance once you account for having searched all of them?

It reuses the permutation null [sca_test()] already computed (so it never re-runs the permutations) and requires that result to have been produced with 'keep_curves = TRUE' and a confound-preserving null ('null_type = "freedman_lane"' or '"residual_bootstrap"'). When those prerequisites are met, [sca_test()] already attaches the single-step adjustment automatically; call 'sca_minp()' to recompute it with a different method or significance threshold.

Usage

sca_minp(result, method = c("single_step", "step_down"), alpha = NULL)

Arguments

Details

Reporting a specification as "significant" because its ordinary p-value is below 0.05 is misleading when the curve contains many correlated specifications: searching dozens of them all but guarantees some will look significant by chance. 'sca_minp()' corrects for this with the min-P / max-statistic permutation method of Westfall and Young (1993). For each permutation it records the most extreme specification anywhere in the curve, building the null distribution of "the best result a search could turn up by chance", and compares each observed specification to that distribution.

Each specification is compared to its *own* permutation null rather than to zero. This matters under the confound-preserving nulls: an under-controlled specification's null distribution is centred on its confounded value, not on zero, so a specification is flagged only when its estimate is more extreme than that specification could itself produce under the null. Consequently a flagged under-controlled specification indicates an association beyond the superset-conditional sharp null – it is **not** a causal effect of size (estimate minus null centre), because omitting a confounder re-routes part of the focal variable's coefficient.

'method = "single_step"' (the default, and what [sca_test()] attaches automatically) controls the family-wise error rate in the weak sense, under the global null of no effect in any specification, with no further assumptions. ‘method = "step_down"' applies Westfall and Young’s free step-down refinement, which is more powerful and controls the family-wise error rate in the strong sense under subset pivotality (exact for '"freedman_lane"' via the common control-superset conditioning; approximate for '"residual_bootstrap"' and under partial alternatives). Adjusted p-values have a resolution floor of '1 / (n_eff + 1)', where 'n_eff' is the number of permutations that estimated at least one specification (equal to 'n_used' unless some permutations were wholly non-estimable).

Value

The input 'result', augmented with a 'fwer' element under 'null_curves': a list with 'specs' (a data frame with one row per specification giving its term-set key 'spec', 'observed' coefficient, uncorrected own-null p-value 'p_raw', FWER-adjusted p-value 'p_adj', and logical 'significant_adj') and 'summary' (the method, alpha, whether weak or strong FWER is controlled, the number of distinct specifications tested, the number significant after correction, and the smallest adjusted p-value).

References

Westfall, P. H., & Young, S. S. (1993). Resampling-based multiple testing: Examples and methods for p-value adjustment. Wiley.

See Also

[sca_test()], which attaches the single-step adjustment automatically; [plot_sca_test_specs()] to see which specifications survive correction.

Examples

result <- sca_test(y = "Salnty", x = "T_degC",
                   controls = c("ChlorA", "O2Sat", "NO2uM"), data = bottles,
                   null_type = "freedman_lane", n_permutations = 199,
                   keep_curves = TRUE, seed = 1, progress_bar = FALSE)
# The single-step adjustment is already attached; recompute it step-down.
result <- sca_minp(result, method = "step_down")
as.data.frame(result, what = "specs")

Write a results paragraph for a specification curve analysis

Description

'sca_report()' returns a one-paragraph, manuscript-ready description of a specification curve. For an [sca_test()] object it states the median estimate, the share of significant specifications, the joint-inference statistics with their permutation p-values, and (when the full set of statistics is available) whether the null of no effect is rejected. For an [sca()] curve it gives a purely descriptive summary.

Usage

sca_report(x, digits = 3, ...)

Arguments

Value

A length-one character string.

See Also

[sca_table()] for a tabular summary.

Examples

s <- sca(y = "Salnty", x = "T_degC", controls = c("ChlorA", "O2Sat"),
         data = bottles, progress_bar = FALSE, parallel = FALSE)
sca_report(s)

result <- sca_test(y = "Salnty", x = "T_degC", controls = c("ChlorA", "O2Sat"),
                   data = bottles, n_permutations = 50, seed = 1,
                   progress_bar = FALSE)
sca_report(result)

Render a specification curve result as a publication table

Description

'sca_table()' turns the output of [sca()] or [sca_test()] into a compact results block suitable for a manuscript: the number of specifications, the median estimate, the share of significant specifications, and (for [sca_test()]) the joint-inference statistics and their permutation p-values.

The default 'format = "data.frame"' returns a small two-column ('label' / 'value') data frame with a 'print' method and adds no dependencies. The other formats render that frame with an optional package: '"markdown"' / '"latex"' use 'knitr', and '"gt"', '"kableExtra"', and '"flextable"' use the package of the same name (install it first). A [tidy()]/[glance()]-aware tool such as 'modelsummary' can consume the object directly without 'sca_table()'.

Usage

sca_table(
  x,
  format = c("data.frame", "markdown", "latex", "gt", "kableExtra", "flextable"),
  digits = 3,
  ...
)

Arguments

Value

For 'format = "data.frame"', a data frame of class '"sca_table"'; otherwise an object of the corresponding rendering package.

See Also

[sca_report()] for a prose summary; [tidy()]/[glance()] for tidy data frames.

Examples

s <- sca(y = "Salnty", x = "T_degC", controls = c("ChlorA", "O2Sat"),
         data = bottles, progress_bar = FALSE, parallel = FALSE)
sca_table(s)

result <- sca_test(y = "Salnty", x = "T_degC", controls = c("ChlorA", "O2Sat"),
                   data = bottles, n_permutations = 50, seed = 1,
                   progress_bar = FALSE)
sca_table(result)

Joint-inference test for a specification curve

Description

'sca_test()' performs the permutation-based joint-inference test of Simonsohn, Simmons, and Nelson (2020). It tests the sharp null hypothesis that the focal independent variable 'x' has no effect on 'y' in any specification. The focal variable is repeatedly permuted (shuffled) and the entire specification curve is re-estimated on each permuted data set, building a null distribution for curve-level summary statistics. The observed statistics are then compared to that null distribution to obtain permutation p-values.

The permutation is *blocked within the first fixed effect* when 'fixed_effects' are supplied, because exchangeability of 'x' only holds within fixed-effect groups; otherwise it is a global permutation.

Usage

sca_test(
  y,
  x,
  controls,
  data,
  weights = NULL,
  family = "linear",
  link = NULL,
  fixed_effects = NULL,
  n_permutations = 500,
  test_stats = c("median", "share_significant", "stouffer"),
  direction = "two.sided",
  alpha = 0.05,
  sca_data = NULL,
  keep_curves = FALSE,
  null_type = c("shuffle_x", "freedman_lane", "residual_bootstrap"),
  common_sample = FALSE,
  parallel = FALSE,
  workers = 2,
  seed = NULL,
  progress_bar = TRUE
)

Arguments

Value

An object of class '"sca_test"': a list with elements 'observed' (a named vector of observed statistics), 'null_distribution' (a data frame with one row per usable permutation and one column per statistic, plus 'n_valid'), 'p_values' (a named vector of permutation p-values), 'params' (a list of run metadata), and, when 'keep_curves = TRUE', 'null_curves' (a list with 'spec', a data frame of each specification's key and observed coefficient, and 'null_coef', a specifications-by-permutations matrix of the permuted focal coefficients). When 'keep_curves = TRUE' *and* the null is confound-preserving ('null_type = "freedman_lane"' or '"residual_bootstrap"'), 'null_curves' also gains 'fwer', the per-specification family-wise-error-rate-adjusted p-values attached automatically by [sca_minp()] (see there for its structure).

References

Simonsohn, U., Simmons, J. P., & Nelson, L. D. (2020). Specification curve analysis. Nature Human Behaviour, 4, 1208-1214. doi:10.1038/s41562-020-0912-z

Freedman, D., & Lane, D. (1983). A nonstochastic interpretation of reported significance levels. Journal of Business & Economic Statistics, 1(4), 292-298.

Winkler, A. M., Ridgway, G. R., Webster, M. A., Smith, S. M., & Nichols, T. E. (2014). Permutation inference for the general linear model. NeuroImage, 92, 381-397.

See Also

[plot_sca_test()] to visualise the null distributions of the test statistics, and [plot_sca_test_specs()] for the per-specification null-band plot (requires 'keep_curves = TRUE').

Examples

# Test whether temperature robustly predicts salinity across specifications.
result <- sca_test(y = "Salnty", x = "T_degC",
                   controls = c("ChlorA", "O2Sat"),
                   data = bottles, n_permutations = 100, progress_bar = FALSE)
result

Tidy and glance methods for speccurvieR objects

Description

[generics::tidy()] and [generics::glance()] methods for the data frame returned by [sca()] (class '"sca"') and the object returned by [sca_test()]. These are the same generics that 'broom::tidy()' / 'broom::glance()' and 'modelsummary' dispatch on, so the results can feed those tools.

Note that the tidied column names follow the broom convention ('estimate' / 'std.error' / 'p.value'), which differs from the raw column names of [sca()] output ('coef' / 'se' / 'p'). For an 'sca' curve every row shares the same focal 'term', so a tool that keys on 'term' (such as 'modelsummary') needs each row distinguished (e.g. by 'spec_id'); the per-statistic rows of 'tidy.sca_test()' are already distinct. The 'controls' column lists the model terms included (so an interaction control 'a*b' appears as its expanded terms), which can differ from the user-facing control count in 'glance()$n_controls'.

Usage

## S3 method for class 'sca'
tidy(x, ...)

## S3 method for class 'sca'
glance(x, ...)

## S3 method for class 'sca_test'
tidy(x, ...)

## S3 method for class 'sca_test'
glance(x, ...)

## S3 method for class 'sca_test'
as.data.frame(
  x,
  row.names = NULL,
  optional = FALSE,
  ...,
  what = c("summary", "null", "params", "specs")
)

Arguments

Value

A data frame. For 'tidy()', one row per specification ('sca') or per test statistic ('sca_test'); for 'glance()', a one-row summary.

See Also

[sca_table()] and [sca_report()] for formatted reporting.

Examples

s <- sca(y = "Salnty", x = "T_degC", controls = c("ChlorA", "O2Sat"),
         data = bottles, progress_bar = FALSE, parallel = FALSE)
tidy(s)
glance(s)

Decompose the variance of the specification curve by analytic choice

Description

'sca_variance()' quantifies how much of the variation in the focal coefficient across a specification curve is attributable to each analytic choice – here, the inclusion or exclusion of each control variable. It is the descriptive complement to [sca_test()]: where the joint-inference test asks whether the curve as a whole is more extreme than chance, the variance decomposition asks which modelling choices *drive* the spread of estimates.

The focal coefficient (the 'coef' column of [sca()] output) is regressed on the 0/1 control-inclusion indicators, and the variance explained is partitioned across controls. By default this uses the LMG (Shapley-value) decomposition of R-squared, whose shares are order-invariant and sum exactly to the model R-squared; the remaining '1 - R^2' is reported as a '"Residual"' row capturing interactions among choices and unexplained variation.

Usage

sca_variance(
  sca_data,
  estimate = "coef",
  method = c("lmg", "type2"),
  residual = TRUE
)

Arguments

Details

**Method.** With ‘method = "lmg"' (default) each control’s share is the Shapley value of R-squared: its average marginal contribution to R-squared over all orderings of the controls. These shares are non-negative, order-invariant, and sum to the full-model R-squared, so the reported percentages (controls plus '"Residual"') sum to 100. 'method = "type2"' instead reports each control's "drop-one" partial sum of squares, normalised to sum to 100; because partial sums of squares do not add up under a non-orthogonal design, those percentages are unique contributions rather than an exact partition. (The regression is purely additive – each control, including an interaction control, enters as a single 0/1 indicator with no product terms – so Type II and Type III partial sums of squares coincide; the '"type2"' label is kept only for familiarity.) '"shapley"' is accepted as an alias for '"lmg"'. The LMG cost grows as 2 to the power of the number of varying controls (like ‘sca()'’s own combinatorial growth) and is capped at 20 controls; for larger sets use '"type2"', which is linear.

**Choice unit.** The decomposition is over *controls* (one binary in/out choice each), read from the indicator columns of [sca()] output. An interaction control such as '"a*b"' is a single choice (its column '"a:b"'). Within one [sca()] call the model family, fixed effects, and focal variable are constant, so only control choices (plus '"Residual"') appear.

**Relationship to specr.** This is analogous to 'specr::icc_specs()' / 'specr::plot_variance()' and returns a comparable choice-by-percent table and bar chart. It deliberately differs in method: 'icc_specs()' fits a random-effects model and reads intraclass-correlation-based variance components, which treats each analytic choice as one grouping factor with many levels. speccurvieR's only enumerated choice is each control's binary status, for which two-level random-effect variance components are unreliable, so an LMG R-squared decomposition is used instead and no 'lme4' dependency is added. The output is intentionally *not* labelled "ICC".

**Caveats.** The shares are descriptive, not inferential: they summarise an already-estimated set of point estimates and carry no sampling-uncertainty interpretation. A share measures how much a control's in/out status *moves* the focal coefficient, not the control's predictive importance and not the direction of any shift. For glm or fixed-effects curves the coefficient is on the model's native (e.g. log-odds) scale, so shares are not comparable across families. At least two specifications with a varying control are required; with exactly two controls the three specifications are saturated by the two main effects, so the residual is zero. If an interaction control is supplied *together with* its component main-effect controls, the controls are not cleanly separable and the shares should be interpreted with care.

Value

A data frame with one row per control (and a final '"Residual"' row when 'residual = TRUE'), sorted by descending share, with columns: 'choice' (the control name, or '"Residual"'), 'variance' (the portion of the focal coefficient's variance attributed to that choice), and 'percent' (its share as a percentage; the column sums to 100).

See Also

[plot_variance()] to visualise the decomposition; [sca_test()] for the joint-inference test.

Examples

s <- sca(y = "Salnty", x = "T_degC",
         controls = c("ChlorA", "O2Sat", "NO2uM"),
         data = bottles, progress_bar = FALSE, parallel = FALSE)
sca_variance(s)
sca_variance(s, method = "type2", residual = FALSE)

Prepares the output of 'sca()' for plotting.

Description

Takes in the data frame output by 'sca()' and returns a list with the data frame and labels to make a plot to visualize the controls included in each spec curve model.

Usage

scp(sca_data)

Arguments

Value

A list containing a data frame, control coefficients, and control names.

Examples

scp(sca(y = "Salnty", x = "T_degC", controls = c("ChlorA", "O2Sat"),
        data = bottles, progress_bar=TRUE, parallel=FALSE));

Estimates bootstrapped standard errors for regression models

Description

Takes in a data frame, regression formula, and bootstrapping parameters and estimates bootstrapped standard errors for models with and without fixed effects. The model is refit on 'n_samples' resamples of 'sample_size' observations drawn with replacement, and the standard deviation of the resampled coefficients is rescaled by 'sqrt(sample_size / nrow(data))' to estimate the standard error of the full-sample estimator (an m-out-of-n bootstrap; when 'sample_size' equals 'nrow(data)' this is the ordinary nonparametric bootstrap).

Usage

se_boot(
  data,
  formula,
  n_x,
  n_samples,
  sample_size,
  weights = NULL,
  fam_obj = NULL
)

Arguments

Value

A named list containing bootstrapped standard errors for each coefficient.

Examples

se_boot(data = bottles, formula = "Salnty ~ T_degC + ChlorA + O2Sat",
        n_x = 3, n_samples = 4, sample_size = 300)


se_boot(data = data.frame(x1 = rnorm(50000, mean=4, sd=10),
                          x2 = rnorm(50000, sd=50),
                          ID = rep(1:100, 500),
                          area = rep(1:50, 1000),
                          y = rnorm(50000)),
        formula = "y ~ x1 + x2 | ID",
        n_x = 2, n_samples = 10, sample_size = 1000)

Compare different kinds of standard errors

Description

se_compare() takes in a regression formula (with or without fixed effects), data, and the types of standard errors desired, including clustered, heteroskedasticity-consistent, and bootstrapped. It then returns a data frame with coefficient and standard error estimates for easy comparison and plotting.

Usage

se_compare(
  formula,
  data,
  weights = NULL,
  family = "linear",
  link = NULL,
  types = "all",
  cluster = NULL,
  clustered_only = FALSE,
  fixed_effects_only = FALSE,
  boot_samples = NULL,
  boot_sample_size = NULL,
  ...
)

Arguments

Value

A data frame where row represents an independent variable in the model and each column a type of standard error. Coefficient estimates for each variable are also included (column '"estimate"' for non-fixed effects model and column '"estimate_FE"' for fixed effects models). Columns are automatically named to specify the standard error type.

Some examples:

"iid" = normal standard errors, i.e. assuming homoskedasticity

"CL_FE" = standard errors clustered by the first fixed effect

"bootstrap_k8n300_FE" = bootstrapped standard errors for a fixed effects model where 'boot_samples = 8' and 'boot_sample_size = 300'

"CL_Depth_ID_FE" = standard errors clustered by the variable "Depth_ID" for a model with fixed effects

"HC0_Sta_ID" = HC0 standard errors clustered by the variable "Sta_ID"

"HC0_Depth_ID_BY_Sta_ID" = HC0 standard errors two-way clustered by "Depth_ID" and "Sta_ID"

Note: for fixed effects models the "(Intercept)" row will be all 'NA' because the intercept is not reported by 'feols()' when fixed effects are present.

Examples

se_compare(formula = "Salnty ~ T_degC + ChlorA + O2Sat | Sta_ID",
           data = bottles, types = "all", cluster = c("Depth_ID", "Sta_ID"),
           fixed_effects_only = FALSE, boot_samples=c(4, 8, 10),
           boot_sample_size=c(300, 500))

se_compare(formula = "Salnty ~ T_degC + ChlorA + O2Sat", data = bottles,
           types = "bootstrapped", boot_samples = c(8, 10),
           boot_sample_size = c(300, 500))

se_compare(formula = "Salnty ~ T_degC + ChlorA", data = bottles,
           types = c("HC0", "HC1", "HC3"))

# Two-way (and multiway) clustering: pass a list, where each element names
# the dimensions to cluster on jointly. Here: one-way by Sta_ID, one-way by
# Depth_ID, and two-way by both.
se_compare(formula = "Salnty ~ T_degC + ChlorA", data = bottles,
           types = "HC1",
           cluster = list("Sta_ID", "Depth_ID", c("Sta_ID", "Depth_ID")))

# Logistic regression: compare standard error types for a binary outcome.
bottles$saline <- as.integer(bottles$Salnty >
                               stats::median(bottles$Salnty, na.rm = TRUE))
se_compare(formula = "saline ~ T_degC + ChlorA", data = bottles,
           family = "binomial", types = c("iid", "HC0", "HC3"))

Deprecated functions in speccurvieR

Description

As of version 1.0.0 the package uses snake_case names throughout. These camelCase aliases are retained for backward compatibility: they emit a deprecation warning and forward to their snake_case replacements (translating any deprecated camelCase argument names automatically). They are scheduled for removal in version 2.0.0.

Usage

plotCurve(...)

plotVars(...)

plotRMSE(...)

plotR2Adj(...)

plotAIC(...)

plotDeviance(...)

plotControlDistributions(...)

controlExtractor(...)

unAsIs(...)

Arguments

Value

The value returned by the corresponding replacement function.

A clean, consistent ggplot2 theme for speccurvieR plots.

Description

‘theme_sca()' is the shared theme applied by the package’s plotting functions. It is exported so the same look can be reused or tweaked when customising the plots they return.

Usage

theme_sca(
  base_size = 11,
  base_family = getOption("speccurvieR.base_family", "")
)

Arguments

Value

A ggplot2 theme object.

Examples

library(ggplot2)
ggplot(bottles, aes(T_degC, Salnty)) + geom_point() + theme_sca();

Removes the 'AsIs' class attribute from the input.

Description

Removes the 'AsIs' class attribute from the input. Taken from: <https://stackoverflow.com/a/12866609>

Usage

un_as_is(x)

Arguments

Value

An object without the 'AsIs' class attribute.

Examples

un_as_is(x = I(c(1:4)));