--- title: "Confidence intervals and p-values" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Confidence intervals and p-values} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) library(smartcor) ``` ## 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 payload} r = smart_cor(mtcars$mpg, mtcars$wt) ``` 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. ```{r payload-tidy} 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]) ``` 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. ```{r source-map} 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. ```{r bootstrap-crosscheck} 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")] ``` 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. ```{r conf-level} 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")] ) ``` ## 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. ```{r nulls} 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. ```{r compare} compare_methods(mtcars$gear, mtcars$carb) ``` 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)`. ```{r linearity-p} lin = test_linearity(mtcars$wt, mtcars$mpg) unlist(lin[c("p_y_on_x", "p_x_on_y", "p.value")]) ``` 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.