Correlations
Aaron R. Caldwell
2026-10-06
correlations.RmdThe TOSTER package provides several functions for calculating and
analyzing correlations. These functions extend beyond traditional
correlation tests by offering equivalence testing capabilities and
robust correlation methods. The included functions are based on research
by Goertzen and Cribbie (2010) (z_cor_test &
compare_cor), and Wilcox (2011)
(boot_cor_test)1.
Simple Correlation Test
Basic tests of association can be performed with the
z_cor_test function. This function is styled after R’s
built-in cor.test function but uses Fisher’s z
transformation as the basis for all significance tests (p-values).
Despite this difference in methodology, the confidence intervals are
typically very similar to those produced by cor.test.
##
## Pearson's product-moment correlation
##
## data: mtcars$mpg and mtcars$qsec
## t = 2.5252, df = 30, p-value = 0.01708
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.08195487 0.66961864
## sample estimates:
## cor
## 0.418684
# TOSTER's z-transformed correlation test
z_cor_test(mtcars$mpg, mtcars$qsec)##
## Pearson's product-moment correlation with approximate SE
##
## data: mtcars$mpg and mtcars$qsec
## z = 2.4023, N = 32, p-value = 0.01629
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.08195487 0.66961864
## sample estimates:
## r
## 0.418684Like cor.test, the z_cor_test function
supports Spearman and Kendall correlation coefficients:
# Spearman correlation
z_cor_test(mtcars$mpg,
mtcars$qsec,
method = "spear") # Short form accepted; "spearman" also works##
## Spearman's rank correlation rho with approximate SE
##
## data: mtcars$mpg and mtcars$qsec
## z = 2.5882, N = 32, p-value = 0.009647
## alternative hypothesis: true rho is not equal to 0
## 95 percent confidence interval:
## 0.1222486 0.7111100
## sample estimates:
## rho
## 0.4669358
# Kendall correlation
z_cor_test(mtcars$mpg,
mtcars$qsec,
method = "kendall")##
## Kendall's rank correlation tau with approximate SE
##
## data: mtcars$mpg and mtcars$qsec
## z = 2.6134, N = 32, p-value = 0.008964
## alternative hypothesis: true tau is not equal to 0
## 95 percent confidence interval:
## 0.08145572 0.51634821
## sample estimates:
## tau
## 0.3153652Advantages of z_cor_test
The main advantage of z_cor_test over the standard
cor.test is its ability to perform equivalence testing
(TOST) or any hypothesis test where the null hypothesis isn’t zero. This
makes it particularly useful for research questions focused on
demonstrating practical equivalence or testing against specific
correlation thresholds. The main disadvantage is that it relies on the
Fisher’s z transformation, which can be less accurate when assumptions
are violated or when outliers are present. In such cases, the
bootstrapped methods may provide more reliable results.
# Equivalence test with null boundary of 0.4
z_cor_test(mtcars$mpg,
mtcars$qsec,
alternative = "e", # e for equivalence
null = .4)##
## Pearson's product-moment correlation with approximate SE
##
## data: mtcars$mpg and mtcars$qsec
## z = 0.12088, N = 32, p-value = 0.5481
## alternative hypothesis: equivalence
## null values:
## correlation correlation
## 0.4 -0.4
## 90 percent confidence interval:
## 0.1397334 0.6360650
## sample estimates:
## r
## 0.418684In this example, we’re testing whether the correlation is equivalent to zero within the boundaries of ±0.4.
Using Summary Statistics
A key advantage of TOSTER is the ability to perform correlation tests
using only summary statistics, which is particularly useful when
reviewing published literature or working with limited data access. The
corsum_test function enables this functionality:
# Testing a correlation of 0.121 from a sample of 105 paired observations
corsum_test(r = .121,
n = 105,
alternative = "e",
null = .4)##
## Pearson's product-moment correlation with approximate SE
##
## data: x and y
## z = -3.0506, N = 105, p-value = 0.001142
## alternative hypothesis: equivalence
## null values:
## correlation correlation
## 0.4 -0.4
## 90 percent confidence interval:
## -0.0412456 0.2770284
## sample estimates:
## r
## 0.121This example tests whether a correlation of 0.121 from a sample of 105 paired observations is equivalent to zero within the boundaries of ±0.4.
Bootstrapped Correlation Test
For more robust analyses when raw data is available, TOSTER provides
the boot_cor_test function. This bootstrapping approach
generally produces more reliable results than Fisher’s z-based tests,
especially when outliers are present or distribution assumptions are
violated.
set.seed(993) # Setting seed for reproducibility
boot_cor_test(mtcars$mpg,
mtcars$qsec,
alternative = "e",
null = .4)##
## Bootstrapped Pearson's product-moment correlation (studentized)
##
## data: mtcars$mpg and mtcars$qsec
## N = 32, p-value = 0.5233
## alternative hypothesis: equivalence
## null values:
## correlation correlation
## 0.4 -0.4
## 90 percent confidence interval:
## 0.2232696 0.5563482
## sample estimates:
## r
## 0.418684
# Bootstrapped Spearman correlation
boot_cor_test(mtcars$mpg,
mtcars$qsec,
method = "spear",
alternative = "e",
null = .4)##
## Bootstrapped Spearman's rank correlation rho (BCa)
##
## data: mtcars$mpg and mtcars$qsec
## N = 32, p-value = 0.6765
## alternative hypothesis: equivalence
## null values:
## rho rho
## 0.4 -0.4
## 90 percent confidence interval:
## 0.1906086 0.6610353
## sample estimates:
## rho
## 0.4669358
# Bootstrapped Kendall correlation
boot_cor_test(mtcars$mpg,
mtcars$qsec,
method = "ken", # Short form accepted
alternative = "e",
null = .4)##
## Bootstrapped Kendall's rank correlation tau (BCa)
##
## data: mtcars$mpg and mtcars$qsec
## N = 32, p-value = 0.2033
## alternative hypothesis: equivalence
## null values:
## tau tau
## 0.4 -0.4
## 90 percent confidence interval:
## 0.1060556 0.4765220
## sample estimates:
## tau
## 0.3153652Bootstrap Confidence Intervals
The boot_ci argument selects the bootstrap confidence
interval, and the p-value is always obtained by inverting that same
interval, so the two always agree.
-
"auto"(default):"stud"for Pearson’s r and"bca"for all other correlations. In simulations, the studentized interval kept Pearson’s r near nominal coverage for skewed, heteroscedastic, and heavy-tailed data, where BCa under-covered; for the rank and robust correlations BCa did as well or better. -
"bca": bias-corrected and accelerated percentile interval. -
"perc": percentile interval. -
"basic": basic (reflected percentile) interval. -
"stud": studentized (bootstrap-t) interval, available for Pearson, Spearman, and Kendall correlations.
The method that was used is returned in the boot_ci
element of the result. The Pearson examples above therefore used the
studentized interval, and the Spearman and Kendall examples used BCa.
Here is a studentized interval for Spearman’s rho:
set.seed(993)
boot_cor_test(mtcars$mpg,
mtcars$qsec,
method = "spearman",
boot_ci = "stud",
alternative = "e",
null = .4)##
## Bootstrapped Spearman's rank correlation rho (studentized)
##
## data: mtcars$mpg and mtcars$qsec
## N = 32, p-value = 0.6838
## alternative hypothesis: equivalence
## null values:
## rho rho
## 0.4 -0.4
## 90 percent confidence interval:
## 0.1764676 0.6526568
## sample estimates:
## rho
## 0.4669358The percentile and BCa intervals give the same answer on any monotone
scale. The basic and studentized intervals do not, so the
boot_scale argument sets the scale on which they are
computed: the Fisher z scale (boot_scale = "z", the
default), or the correlation scale (boot_scale = "r"). The
z scale keeps the interval within \([-1,
1]\) and is closer to variance-stabilizing.
Standard Errors for the Studentized Bootstrap
The studentized bootstrap standardizes every bootstrap replicate by its own standard error, \(t^* = (\hat\theta^* - \hat\theta) / \widehat{SE}^*\), and uses the distribution of \(t^*\) in place of a normal or \(t\) reference distribution. This only improves on the basic interval if \(\widehat{SE}^*\) is estimated from the data in each resample. The usual normal-theory standard errors of Fisher’s z (e.g., \(1/\sqrt{n-3}\) for Pearson’s r) depend only on the sample size and so cannot be used for this.
boot_cor_test() therefore uses an influence-function
(sandwich) standard error, which is computed from the data and does not
assume bivariate normality:
\[ \widehat{SE}(\hat\rho) = \frac{1}{n}\sqrt{\sum_{i=1}^{n} \psi_i^2}, \]
where \(\psi_i\) is the empirical influence of observation \(i\) on the coefficient. On the z scale the delta method gives \(\widehat{SE}(\hat z) = \widehat{SE}(\hat\rho) / (1 - \hat\rho^2)\). The influence functions are:
Pearson’s r. The starting point is the asymptotic distribution-free (fourth-moment) influence function (Steiger and Hakstian 1982). With \(z_{x,i}\) and \(z_{y,i}\) the standardized observations (using denominator \(n\)),
\[ \psi^0_i = z_{x,i} z_{y,i} - \frac{r}{2}\left(z_{x,i}^2 + z_{y,i}^2\right). \]
Under bivariate normality this reduces to the familiar \(\text{Var}(r) \approx (1-\rho^2)^2/n\). On its own, though, it gives standard errors that are too small in small samples and with heavy-tailed data, because a sample under-represents the extreme points that dominate the variance of \(r\). Following the HC4 heteroscedasticity-consistent estimator for regression (Cribari-Neto 2004), which Wilcox (2011) also uses for testing correlations, each influence value is inflated according to its leverage:
\[ \psi_i = \psi^0_i\,(1 - h_i)^{-\delta_i/2}, \qquad h_i = \frac{1}{n} + \frac{D_i^2}{n-1}, \qquad \delta_i = \min\left(4, \frac{n h_i}{3}\right), \]
where \(D_i\) is the Mahalanobis distance of \((x_i, y_i)\) from the bivariate mean (so the \(h_i\) sum to 3), and \(h_i\) is capped at 0.99. This leverage correction of the ADF influence function is an adaptation made for TOSTER. It was chosen because in simulations (normal, heteroscedastic, \(t_5\), \(t_3\), and discrete data; \(n\) = 30 and 80) it kept the studentized interval’s coverage near nominal, including for \(t_3\) data, where the uncorrected standard error under-covered badly. It is slightly conservative for strongly heteroscedastic data. When the data have infinite fourth moments, no standard error of \(r\) is consistent, and a robust correlation (below) is still the better choice.
Spearman’s rho. Spearman’s rho is Pearson’s r computed on the mid-distribution transforms \(u_i = (R(x_i) - 1/2)/n\) and \(v_i = (R(y_i) - 1/2)/n\), where \(R\) denotes midranks. Its influence function adds terms for estimating these transforms to the Pearson influence function of \((u, v)\). With \(\omega^x_{ij} = \mathbf{1}(x_j > x_i) + \tfrac{1}{2}\mathbf{1}(x_j = x_i)\) (and \(\omega^y_{ij}\) likewise), the influence of observation \(i\) on each moment is
\[ \begin{aligned} d_i(\overline{uv}) &= u_i v_i + \tfrac{1}{n}\textstyle\sum_j \omega^x_{ij} v_j + \tfrac{1}{n}\sum_j \omega^y_{ij} u_j - 3\,\overline{uv}, \\ d_i(\bar u) &= u_i + \tfrac{1}{n}\textstyle\sum_j \omega^x_{ij} - 2\bar u, \\ d_i(\overline{u^2}) &= u_i^2 + \tfrac{2}{n}\textstyle\sum_j \omega^x_{ij} u_j - 3\,\overline{u^2}, \end{aligned} \]
(similarly for \(v\)), and the delta method for \(r_s = C/\sqrt{S_u S_v}\), with \(C = \overline{uv} - \bar u \bar v\) and \(S_u = \overline{u^2} - \bar u^2\), gives
\[ \psi_i = \frac{d_i(C)}{\sqrt{S_u S_v}} - \frac{r_s}{2}\left(\frac{d_i(S_u)}{S_u} + \frac{d_i(S_v)}{S_v}\right), \]
where \(d_i(C) = d_i(\overline{uv}) - \bar v\, d_i(\bar u) - \bar u\, d_i(\bar v)\) and \(d_i(S_u) = d_i(\overline{u^2}) - 2\bar u\, d_i(\bar u)\). Without ties this is the influence function of Spearman’s rho derived by Croux and Dehon (2010). The tie terms keep the standard error accurate for discrete data, and for bootstrap resamples, which always contain ties.
Kendall’s tau. Kendall’s tau-b is a ratio of U-statistics, \(\tau_b = \bar a / \sqrt{\bar b^x \bar b^y}\), where
\[ a_i = \frac{1}{n-1}\sum_{j \ne i} \text{sign}(x_i - x_j)\,\text{sign}(y_i - y_j), \qquad b^x_i = \frac{1}{n-1}\sum_{j \ne i} \mathbf{1}(x_i \ne x_j), \]
and \(b^y_i\) is defined likewise. Hoeffding’s projection (Hoeffding 1948) and the delta method give
\[ \psi_i = 2\left[\frac{a_i - \bar a}{\sqrt{\bar b^x \bar b^y}} - \frac{\tau_b}{2}\left(\frac{b^x_i - \bar b^x}{\bar b^x} + \frac{b^y_i - \bar b^y}{\bar b^y}\right)\right]. \]
Without ties this is the usual U-statistic variance, \(\frac{4}{n}\widehat{\text{Var}}(a_i)\). Computing it takes \(O(n^2)\) operations per resample, so the studentized Kendall interval is slower in large samples.
Under independence, all three standard errors reduce to the familiar
values, \(\text{Var}(r) \approx
\text{Var}(r_s) \approx 1/n\) and \(\text{Var}(\tau) \approx 4/(9n)\). Because
each pivot divides by an estimated standard error, studentized intervals
can be wide or erratic in small samples (roughly \(n < 30\)). This matters most for the
rank correlations. Kendall’s standard error, for example, is nearly
unbiased even at \(n = 20\), but it
varies a lot between samples (especially with ties) and is strongly
correlated with the estimate. In simulations at \(n = 30\), the studentized Spearman and
Kendall intervals under-covered by up to about 2.5 percentage points for
heavy-tailed or tied data, while the BCa interval stayed closer to
nominal; at \(n = 80\) both were
accurate. For Spearman’s rho and Kendall’s tau with fewer than about 50
observations, use boot_ci = "bca", which is what
"auto" selects for these methods. The studentized interval
is not available for the robust correlations below, which lack a
closed-form standard error.
Robust Correlation Methods
The boot_cor_test function also provides access to
robust correlation methods that are less sensitive to outliers and
violations of normality:
# Winsorized correlation with 10% trimming
boot_cor_test(mtcars$mpg,
mtcars$qsec,
method = "win",
alternative = "e",
null = .4,
tr = .1) # Set trim amount (default is 0.2)##
## Bootstrapped Winsorized correlation wincor (BCa)
##
## data: mtcars$mpg and mtcars$qsec
## N = 32, p-value = 0.6623
## alternative hypothesis: equivalence
## null values:
## wincor wincor
## 0.4 -0.4
## 90 percent confidence interval:
## 0.1867593 0.6465992
## sample estimates:
## wincor
## 0.464062
# Percentage bend correlation
boot_cor_test(mtcars$mpg,
mtcars$qsec,
method = "bend",
alternative = "e",
null = .4,
beta = .15) # Beta parameter controlling resistance to outliers##
## Bootstrapped percentage bend correlation pb (BCa)
##
## data: mtcars$mpg and mtcars$qsec
## N = 32, p-value = 0.6037
## alternative hypothesis: equivalence
## null values:
## pb pb
## 0.4 -0.4
## 90 percent confidence interval:
## 0.1720451 0.6243750
## sample estimates:
## pb
## 0.4484488The Winsorized correlation reduces the impact of outliers by replacing extreme values with less extreme values. The percentage bend correlation is another robust method that downweights the influence of outliers in the calculation.
Comparing Correlations
TOSTER provides tools for comparing correlations between independent groups or studies. This is useful for testing differences in relationships across populations or for evaluating replication studies.
Summary Statistics Approach
When only summary statistics are available, the
compare_cor function can be used:
# Comparing correlation r1=0.8 from n=40 with r2=0.2 from n=100
compare_cor(r1 = .8,
df1 = 38, # df = n-2
r2 = .2,
df2 = 98) # df = n-2##
## Difference between two independent correlations (Fisher's z transform)
##
## data: Summary Statistics
## z = 4.6364, p-value = 3.545e-06
## alternative hypothesis: true difference between correlations is not equal to 0
## sample estimates:
## difference between correlations
## 0.6The compare_cor function supports different methods for
comparing correlations:
# Testing equivalence using Fisher's method
compare_cor(r1 = .8,
df1 = 38,
r2 = .2,
df2 = 98,
null = .2,
method = "f", # Fisher (can also use "fisher")
alternative = "e") # Equivalence test##
## Difference between two independent correlations (Fisher's z transform)
##
## data: Summary Statistics
## z = 3.5872, p-value = 0.9998
## alternative hypothesis: equivalence
## null values:
## difference between correlations difference between correlations
## 0.2 -0.2
## sample estimates:
## difference between correlations
## 0.6Available methods include:
-
Fisher’s z transformation
(
method = "fisher"or"f"): Tests the difference between correlations on the z-transformed scale. This is generally recommended for most applications. -
Kraatz’s method (
method = "kraatz"or"k"): Directly measures the difference between correlation coefficients.
While both methods are appropriate for general significance testing, they may have limited statistical power in some scenarios (Counsell and Cribbie 2015).
Bootstrapped Comparison
When raw data is available for both correlations, the
boot_compare_cor function offers a more robust approach
through bootstrapping:
set.seed(8922) # Setting seed for reproducibility
# Generating example data
x1 = rnorm(40)
y1 = rnorm(40)
x2 = rnorm(100)
y2 = rnorm(100)
# Bootstrap comparison with winsorized correlation
boot_compare_cor(
x1 = x1,
x2 = x2,
y1 = y1,
y2 = y2,
null = .2,
alternative = "e", # Equivalence test
method = "win" # Winsorized correlation
)##
## Bootstrapped difference in Winsorized correlation wincor
##
## data: x1 and y1 vs. x2 and y2
## n1 = 40, n2 = 100, p-value = 0.7739
## alternative hypothesis: true differnce in wincor is 0.2
## 90 percent confidence interval:
## -0.2970547 0.3978333
## sample estimates:
## wincor
## 0.06383164This approach has several advantages:
- It does not rely on the Fisher’s z-transformation approximation
- It can incorporate robust correlation methods
- It can provide more accurate confidence intervals, especially when typical assumptions are violated
Practical Recommendations
When choosing which correlation method to use in TOSTER:
-
If raw data is available:
- For most cases, use
boot_cor_testwith Pearson, Spearman, or Kendall methods - When outliers or distribution assumptions are concerns, consider the robust methods (winsorized or percentage bend)
- For most cases, use
-
If only summary statistics are available:
- Use
corsum_testfor single correlation analysis - Use
compare_corwith the Fisher method for comparing correlations
- Use
-
For equivalence testing:
- Carefully select meaningful boundaries (null values) based on your research context
- Consider what effect size would be practically insignificant in your field