Confidence intervals and p-values

What smartcor reports

A correlation coefficient alone does not tell you much. You also need to know how precise it is and what it is being tested against. So every smart_cor() result includes a confidence interval, a p-value, the null hypothesis, and the name of the method used to compute each of them. The interval and the test depend on the method. For example, a Kendall interval is not computed the same way as a Pearson interval.

An example

Take two continuous variables.

r = smart_cor(mtcars$mpg, mtcars$wt)
#> 
#> ── Variable types ──
#> 
#> ℹ Detected as continuous (numeric with 25 unique values).
#> mtcars$mpg: "continuous"
#> ℹ Detected as continuous (numeric with 29 unique values).
#> mtcars$wt: "continuous"
#> 
#> ── Method selection ──
#> 
#> ✔ Selected: Pearson Correlation

The printed output shows the interval, the p-value, the test, the null hypothesis, and where each number comes from. tidy() returns the same information as columns.

library(tibble)  # for printing
infcols = c("estimate", "statistic", "p.value", "p_method", "null_hypothesis",
             "ci_lower", "ci_upper", "conf_level", "ci_method", "ci_source")
as.data.frame(tidy(r)[infcols])
#>     estimate statistic      p.value               p_method null_hypothesis
#> 1 -0.8676594 -9.559044 1.293959e-10 t-test on r (cor.test)         rho = 0
#>     ci_lower   ci_upper conf_level ci_method            ci_source
#> 1 -0.9249087 -0.7906939       0.95 bootstrap Percentile bootstrap

Two of these columns describe the interval. ci_method is "analytic" for a closed-form or asymptotic interval and "bootstrap" for a resampled one. ci_source names the formula or procedure behind it. With these two columns you can check how an interval was computed, not only its endpoints.

Analytic intervals by default

Fifteen of the eighteen methods have a closed-form or asymptotic confidence interval, and smartcor uses it by default. The other three, Theil’s U, the Kruskal-Wallis \(\eta^2_H\), and Freeman’s theta, have no standard analytic interval, so smartcor bootstraps them. The table below is built by the package itself. For each method we force that estimator on a suitable pair and read back how its interval was computed, so the table always matches the code.

set.seed(1)
colour = factor(sample(c("red", "blue", "green"), 90, replace = TRUE))
shape  = factor(sample(c("circle", "square", "triangle"), 90, replace = TRUE))
cyl    = factor(mtcars$cyl)

one = function(method, x, y) {
  tidy(smart_cor(x, y, method = method, seed = 1, verbose = FALSE))[
    c("method_label", "ci_method", "ci_source", "p_method", "null_hypothesis")]
}

spec = list(
  one("pearson",        mtcars$mpg,  mtcars$wt),
  one("spearman",       mtcars$mpg,  mtcars$wt),
  one("kendall",        mtcars$mpg,  mtcars$wt),
  one("point_biserial", mtcars$mpg,  mtcars$vs),
  one("biserial",       mtcars$mpg,  mtcars$vs),
  one("phi",            mtcars$vs,   mtcars$am),
  one("tetrachoric",    mtcars$vs,   mtcars$am),
  one("yules_q",        mtcars$vs,   mtcars$am),
  one("polychoric",     mtcars$gear, mtcars$carb),
  one("polyserial",     mtcars$mpg,  mtcars$gear),
  one("gamma",          mtcars$gear, mtcars$carb),
  one("rank_biserial",  mtcars$vs,   mtcars$gear),
  one("eta_squared",    mtcars$mpg,  cyl),
  one("kw_eta_squared", mtcars$mpg,  cyl),
  one("freemans_theta", cyl,         mtcars$gear),
  one("cramers_v",      colour,      shape),
  one("theils_u",       colour,      shape),
  one("tschuprows_t",   colour,      shape)
)
map = do.call(rbind, spec)
knitr::kable(map[c("method_label", "ci_method", "ci_source")],
             row.names = FALSE, caption = "How each method's confidence interval is computed.")

The analytic intervals come from the usual sources:

• Fisher’s \(z\) transform for Pearson and point-biserial.
• Bonett and Wright (2000) for Spearman.
• Fieller, Hartley, and Pearson (1957) for Kendall.
• The asymptotic standard errors of Olsson (1979) and Olsson, Drasgow, and Dorans (1982), on the Fisher-\(z\) scale, for polychoric, tetrachoric, polyserial, and biserial.
• Brown and Benedetti (1977) for Goodman-Kruskal’s gamma and Yule’s Q.
• Cliff (1996) for rank-biserial.
• A noncentral chi-square interval for phi, Cramer’s V, and Tschuprow’s T.
• A noncentral \(F\) interval for the correlation ratio \(\eta^2\).

The noncentral intervals invert the chi-square or \(F\) test for the noncentrality parameter. They are asymmetric and never go below zero, which matters when the association is weak. The point-biserial interval uses the Fisher-\(z\) variance \(1/(n-3)\), which assumes bivariate normality. A binary variable cannot be normal, so this interval is only approximate.

The four estimators added in 1.2.0

• Biserial. The estimate is the polyserial correlation with the binary variable treated as an ordinal variable with two categories. By default it uses polycor::polyserial(..., std.err = TRUE); we recommend polyserial_ml = TRUE, which fits the full maximum-likelihood estimate. The p-value is a Wald test on the Fisher-\(z\) scale using the standard error from the fit, and the interval is the matching Wald interval transformed back. If the fit fails, smartcor reports the point-biserial coefficient instead and says so in method and rationale.
• Correlation ratio \(\eta^2\). The p-value comes from the one-way ANOVA \(F\) test. The interval inverts the noncentral \(F\) distribution for the noncentrality parameter \(\lambda\) and maps each bound to \(\lambda / (\lambda + n - k)\), as effectsize::eta_squared() does for a two-sided interval.
• Kruskal-Wallis \(\eta^2_H\). The p-value comes from the Kruskal-Wallis test. The interval is a percentile bootstrap.
• Freeman’s theta. The p-value is a one-sided permutation test. It shuffles the ordinal variable n_perm times (500 by default) and recomputes theta, as we do for Theil’s U. The interval is a percentile bootstrap.

With bootstrap = FALSE, the last two have no interval and return NA. Both use seed when resampling.

The bootstrap argument

bootstrap takes three values.

• "auto" (the default) uses the analytic interval when one exists and bootstraps otherwise. There is one exception: when the linearity rule recommends robust inference for Pearson, it uses a paired percentile-bootstrap interval.
• TRUE always uses a percentile-bootstrap interval.
• FALSE uses only the analytic interval and returns NA when there is none.

The bootstrap only changes the confidence interval. The p-value stays the same. For Pearson, the p-value is the standard Pearson test even when the rationale recommends robust inference.

n_boot sets the number of resamples (500 by default). Below we bootstrap the same Pearson pair and compare the interval with the Fisher-\(z\) interval above.

rb = smart_cor(mtcars$mpg, mtcars$wt, bootstrap = TRUE, n_boot = 1000, verbose = FALSE)
tidy(rb)[c("estimate", "ci_lower", "ci_upper", "ci_method", "ci_source")]
#> # A tibble: 1 × 5
#>   estimate ci_lower ci_upper ci_method ci_source           
#>      <dbl>    <dbl>    <dbl> <chr>     <chr>               
#> 1   -0.868   -0.926   -0.792 bootstrap Percentile bootstrap

The two intervals are close, so the Fisher-\(z\) approximation works well for this pair. Set seed if you want the bootstrap interval to be reproducible. A call with seed restores your .Random.seed when it finishes. conf_level sets the confidence level, and a lower level gives a narrower interval.

rbind(
  `90%` = tidy(smart_cor(mtcars$mpg, mtcars$wt, conf_level = 0.90, verbose = FALSE))[c("ci_lower", "ci_upper")],
  `95%` = tidy(smart_cor(mtcars$mpg, mtcars$wt, conf_level = 0.95, verbose = FALSE))[c("ci_lower", "ci_upper")]
)
#> # A tibble: 2 × 2
#>   ci_lower ci_upper
#> *    <dbl>    <dbl>
#> 1   -0.918   -0.803
#> 2   -0.926   -0.793

Each method tests a different null

A p-value means little unless you know the null hypothesis, and the null is not the same for every method. Pearson, Spearman, Kendall, and the latent-variable methods test whether the (latent) correlation is zero. The chi-square-based measures test whether the two variables are independent. Theil’s U tests whether one variable tells you nothing about the other. The correlation ratio tests whether the group means are equal, the Kruskal-Wallis \(\eta^2_H\) tests whether the groups have the same distribution, and Freeman’s theta tests whether the ordinal variable is independent of the groups. The null_hypothesis and p_method columns state this for every result.

knitr::kable(unique(map[c("method_label", "p_method", "null_hypothesis")]),
             row.names = FALSE, caption = "The null hypothesis and test behind each p-value.")

Comparing interval widths

The width of an interval is as informative as the estimate. compare_methods() lists every method that applies to a pair, each with its own interval.

compare_methods(mtcars$gear, mtcars$carb)
#> 
#> ── Comparing methods ──
#> 
#> mtcars$gear ("ordinal") x mtcars$carb ("ordinal")
#> N: 32
#> Methods: "polychoric", "kendall", "spearman", and "gamma"
#> 
#> "polychoric": 0.2453 <- recommended
#> "kendall": 0.0980
#> "spearman": 0.1149
#> "gamma": 0.1405
#> 
#> ── Method Comparison ───────────────────────────────────────────────────────────
#> Variables: mtcars$gear ("ordinal") x mtcars$carb ("ordinal")
#> N: 32
#> 
#> Polychoric Correlation <- recommended
#> r = 0.2453 p = 0.2487 95% CI [-0.1733, 0.5888]
#> Kendall's Tau-b
#> r = 0.0980 p = 0.5302 95% CI [-0.1455, 0.3303]
#> Spearman Rank Correlation
#> r = 0.1149 p = 0.5312 95% CI [-0.2447, 0.4467]
#> Goodman-Kruskal's Gamma
#> r = 0.1405 p = 0.5433 95% CI [-0.3046, 0.5353]

Here the Goodman-Kruskal gamma interval is much wider than the rank-based ones. Gamma drops all tied pairs, so it uses less of the data and its interval is less precise. Two methods can give similar coefficients and still differ a lot in precision.

P-values of the diagnostics

The diagnostics used to choose the method have their own p-values. tidy() keeps them as normality_p, latent_p, and linearity_p. They describe the tests behind the choice of method, not the correlation itself.

test_bivariate_normality() can combine several tests. For a continuous-ordinal pair, for example, it runs Jarque-Bera and the \(G^2\) test. It rejects if any of them rejects, and its p.value is the smallest of their p-values.

test_linearity() runs in both directions, y on x and x on y. A direction rejects linearity only when a residual-structure test (Box-Pierce, Ljung-Box, or BDS) and RESET both reject, so its p-value is \[p_{dir} = \max\left(\min(p_{BP}, p_{LB}, p_{BDS}),\ p_{RESET}\right).\] The pair is called nonlinear only when both directions reject, so the combined p-value is \(p = \max(p_{y|x}, p_{x|y})\). These are intersection-union p-values, so \(p < \alpha\) exactly when linearity is rejected. The BDS p-value is the smallest over the distance thresholds and embedding dimensions 2 and 3 of fNonlinear::bdsTest(residuals, m = 3).

lin = test_linearity(mtcars$wt, mtcars$mpg)
unlist(lin[c("p_y_on_x", "p_x_on_y", "p.value")])
#>   p_y_on_x   p_x_on_y    p.value 
#> 0.01262732 0.05617456 0.05617456

If Pearson is selected but Mardia’s test rejects bivariate normality, smartcor still reports the Pearson interval and p-value. The rationale then notes that they are approximate because the data are not bivariate normal.

Summary

Every result has a confidence interval, a p-value, a null hypothesis, and the name of the procedure behind each. bootstrap decides whether the interval is resampled; it never changes the p-value. Use ci_method, ci_source, and p_method to see how the interval and the p-value were computed.