Bootstrap Confidence Intervals for 11 Robust Correlations in the Presence of Outliers and Leverage Observations

Approach 1: Replacement by Medians

As noted above, a first robust approach is to directly replace the linear summation of ${{\displaystyle \sum }}_{i=1}^n\left(x_i-\overline{x}\right)\cdot \left(y_i-\overline{y}\right)$ that is sensitive to outliers by a robust counterpart, i.e., ${{\displaystyle \sum }}_{i=n_l}^{n_u}\left[\left(x_i-\text{median}\left(x\right)\right)\left(y_i-\text{median}\left(y\right)\right)\right]$, where $n_l$ and $n_u$ are the lowest and highest ranked scores, respectively in a trimmed sample, e.g., 10% top and 10% bottom are trimmed as suggested in Gnanadesikan and Kettenring (1972). Next, replacing the mean deviation component, ${{\displaystyle \sum }}_{i=1}^n{\left(x_i-\overline{x}\right)}^2\cdot {{\displaystyle \sum }}_{i=1}^n{\left(y_i-\overline{y}\right)}^2$, by a robust counterpart ${{\displaystyle \sum }}_{i=n_l}^{n_u}{\left(x_i-\text{median}\left(x\right)\right)}^2{{\displaystyle \sum }}_{i=n_l}^{n_u}{\left(y_i-\text{median}\left(y\right)\right)}^2$. Consequently, a direct robust sequel to r, called a robust correlational median estimator (r_CME; Falk, 1998; Shevlyakov & Vilchevsky, 2002), can be expressed as

A second type of robust correlation makes use of a deviation estimate, called the median absolute deviation, i.e., $\text{MAD}\left(z\right)=\text{median}\left[\left|z-\text{median}\left(z\right)\right|\right]$, which is more robust than the conventional standard deviation (SD) estimate. With this in mind, a robust MAD-based correlation (r_MAD; Pasman & Shevlyakov, 1987) can be expressed as

where u and v are the robust principal variables, i.e., $u=\left[x-\text{median}\left(x\right)\right]/\sqrt{2}\text{MAD}\left(x\right)+\left[y-\text{median}\left(y\right)\right]/\sqrt{2}\text{MAD}\left(y\right)$ and $v=\left[x-\text{median}\left(x\right)\right]/\sqrt{2}\text{MAD}\left(x\right)-\left[y-\text{median}\left(y\right)\right]/\sqrt{2}\text{MAD}\left(y\right)$.

A third robust correlation is considered a derivative of $r_{MAD}$, and is known as the median correlation coefficient (r_MED) (Shevlyakov & Smirnov, 2011),

where u and v are defined as in Equation (4).

A fourth robust correlation is called the biweight midcorrelation ( $r_{bm}$; Wilcox, 2012), which first converts $x_i$ to $e_i$ by $e_i=\left[x_i-\text{median}\left(x\right)\right]/9\cdot \text{MAD}(x)$ and $y_i$ to $f_i$ by $f_i=\left[y_i-\text{median}\left(y\right)\right]/9\cdot \text{MAD}(y)$. Second, assign $a_i=1$ if $-1\le e_i\le 1$, else $a_i=0$. Similarly, assign $b_i=1$ if $-1\le f_i\le 1$, else $b_i=0$. Consequently,

where $s_e$ and $s_f$ are the SDs of $e_i$ and $f_i$, respectively.

Approach 2: Counting Ranks or Signs of X-Y Pairs

Another approach to robust correlation involves counting the signs or ranks of x-y pairs. This approach is intended to minimize the effect of outliers and/or leverage observations on estimation of bivariate linearity. Blomqvist (1950) developed the quadrant correlation coefficient ( $r_Q$),

where sign refers to the sign of deviations from median, and $\pi \approx 3.14156$. The component, ${{\displaystyle \sum }}_{i=1}^n\text{sign}\left[x_i-\text{median}\left(x\right)\right]\cdot \text{sign}\left[y_i-\text{median}\left(y\right)\right]$, is a count of the cases that when x is above (or below) the median of all xs, its paired y is also above (or below) the median of all ys. The average of the total counts ( $\gamma$) can be scaled back to a number that measures linear association identical to the $\theta$-metric (from -1 to +1) via Shepard’s Theorem, $\text{sin}\left(0.5\cdot \pi \cdot \gamma \right)$ (Kendall & Stuart, 1977).

A second robust correlation based on counting the x-y pairs is called Kendall’s $\tau$ correlation ( $r_{\tau }$; Kendall, 1938),

where $n_C$ refers to the number of concordant pairs for X and Y, and $n_D$ refers to the number of discordant pairs for variables X and Y. Specifically, one can consider any x-y pairs, (x_i, y_i) and (x_j, y_j), where $i\ne j$. If the ranks for both pairs agree (i.e., either x_i > x_j and y_i > y_j or x_i < x_j and y_i < y_j), then they are regarded as concordant pairs; otherwise, they are regarded as discordant pairs (i.e., either x_i > x_j and y_i < y_j or x_i < x_j and y_i > y_j). Note that Equation (8) also involves Shepard’s Theorem, $\text{sin}\left(0.5\cdot \pi \cdot \gamma \right)$, so that the value can be scaled back to the $\theta$-metric (from -1 to +1).

In addition, Spearman’s correlation (r_s; Spearman, 1904) transforms X and Y scores into rank scores, computes the distance between these rank scores, and measures bivariate linearity. That is, when there is no tie,

where $d_i^2={\left[\text{rank}\left(x_i\right)-\text{rank}\left(y_i\right)\right]}^2$ is the squared difference between the rank of a score in X and a score in Y.

Approach 3: Utilizing Weighted Least Squares Estimation

Rousseeuw and Leroy (1987) proposed the weighted least squares correlation (r_wLS). First, a standard linear regression equation is fitted by $y_i={\beta }_0+{\beta }_1{x}_i+e_i$. Second, $r_i=y_i-{{\displaystyle \sum }}_{i=1}^n{\beta }_1{x}_i$ is defined as the residual score for the ith respondent. Third, $s=1.4826\left[1+5/\left(n-1\right)\right]\cdot \sqrt{\text{median}\left(r_i^2\right)}$ is defined as a determining factor. Fourth, one can compute a weight for each $r_i$, i.e., $w_i=1$, if $\left|r_i/s\right|\le 2.5$, else $w_i=0$. Consequently,

where ${\overline{x}}_w={{\displaystyle \sum }}_{i=1}^n{w}_i{x}_i/{{\displaystyle \sum }}_{i=1}^n{w}_i$ and ${\overline{y}}_w={{\displaystyle \sum }}_{i=1}^n{w}_i{y}_i/{{\displaystyle \sum }}_{i=1}^n{w}_i$.

Approach 4: Deletion of Outliers

This approach focuses on discarding a certain percent of the top and bottom x-y pairs which are regarded as O-LO in a sample data-set. According to Huber (1981), a trimmed correlation can be presented as

where $n^{*}=n\cdot \left(1-\pi \right)$ is the number of trimmed observations with $\pi$ proportion, and p_i and q_i are the ith order or rank scores of variables X and Y, respectively. According to Algina et al. (2005), $\pi$ is often set at 20% in order to obtain a reliable robust estimate.

Another trim-based robust correlation is called the Winsorized correlation (r_W; Wilcox, 2012). This method replaces all the values above (or below) a cut-off with a list of maximum (or minimum) scores that are left in the dataset. Specifically, researchers can rank-order the X scores such that x(1^st) > x(2^nd) >…> x(nth). The upper and lower cut-off points become $\left(n\pi \right)$ and $\left(1-n\pi \right)$, where $\pi$ is often fixed at 20% according to Wilcox (2012). Consequently, when $x\left(i\right)\le x\left(n\pi \right)$, the Winsorized scores are converted as $w_X\left(i\right)=x\left(n\pi \right)$. When $x\left(n\pi \right)<x\left(i\right)<x\left(1-n\pi \right)$, $w_X\left(i\right)=x\left(i\right)$. When $x\left(i\right)\ge x\left(1-n\pi \right)$, $w_X\left(i\right)=x\left(1-n\pi \right)$. Repeating the same Winsorized procedure for Y, one can estimate $r_w$ by

A third trimmed correlation is known as the percentage bend correlation (r_pbc; Wilcox, 2012), which provides a linear estimate that is not overly sensitive to slight deviation from a parametric (normal) distribution. According to Wilcox, first compute the median deviation scores $a_i=\left|x_i-\text{median}\left(X\right)\right|$ and rank $a_i$ in ascending order. Second, find the criterion by ${\widehat{c}}_x=a_{\left(m\right)}$, where $m=\left[\left(1-\pi \right)n\right]$ and $\pi =20\%$ is defined as in (8). Third, set $n_1$ be the number of $x_i$ such that $\left[x_i-\text{median}\left(X\right)\right]/{\widehat{c}}_x<-1$, and $n_2$ be the number of $x_i$ such that $\left[x_i-\text{median}\left(X\right)\right]/{\widehat{c}}_x>1$. Fourth, convert $x_i$ to $u_i$ by $u_i=\left(x_i-{\widehat{\varphi }}_x\right)/{\widehat{c}}_x$, where ${\widehat{\varphi }}_x=\left[{\widehat{c}}_x\left(n_2-n_1\right)+b_x\right]/\left(n-n_2-n_1\right)$ and $b_x={{\displaystyle \sum }}_{i=n_1+1}^{n-n_2}{x}_i$. Repeat the same procedure for converting $y_i$ to $v_i$ through $v_i=\left(y_i-{\widehat{\varphi }}_y\right)/{\widehat{c}}_y$, where ${\widehat{\varphi }}_y=\left[{\widehat{c}}_y\left(n_2-n_1\right)+b_y\right]/\left(n-n_2-n_1\right)$ and $b_y={{\displaystyle \sum }}_{i=n_1+1}^{n-n_2}{y}_i$. Fifth, set $g_i=\psi \left(u_i\right)$ and $h_i=\psi \left(v_i\right)$ with $\psi \left(u_i\right)=\text{max}\left[-1,\text{min}\left(1,u_i\right)\right]$ and $\left(v_i\right)=\text{max}\left[-1,\text{min}\left(1,v_i\right)\right]$ Then compute the percentage bend correlation through

Factor 1: Population Bivariate Linearity ( $\theta$)

Four levels— .10, .30, .50, and .80—were evaluated. The first three levels (.10, .30, .50) generally refer to a small, medium, and large effect size, respectively, which are commonly found in behavioral and social science research (Cohen, 1988). In addition, an extremely large value (.80) was also included to examine the impact of a very large effect on the correlation estimates.

Factor 2: Sample Size (n)

Three levels—50, 100, and 200—were examined, which correspond to relatively small, moderate, and large sample sizes frequently used in simulation studies (Li et al., 2011).

Factor 3: Proportion and SD ( $\zeta ,\varrho$)

A total of six conditions that correspond to 3 different proportions of outliers and/or leverage points with 2 different SD units were examined. An ideal condition (i.e., no outliers/leverage points) was also examined, i.e., x and y are normally distributed, and are linearly related with a level of $\theta$. In equation,

where x ~ N(0, 1), and e is the error term such that $e~N\left(0,{\sigma }_e\right)$, where ${\sigma }_e=\sqrt{1-{\theta }^2}$.

The ideal condition in Equation (19) can be violated in behavioral and social sciences. Often, some of the x and/or y observations deviate from the general pattern of the remaining observations. According to previous simulation studies, (e.g., Algina et al., 2005), it is common that 10% of the x and/or y scores are these deviant O-LO with an extreme distance (i.e., 10 times of the SD in original scores) relative to the remaining scores. In this study, the 10% and 10 SD-unit is denoted as (10%, 10-SD). To provide a more complete picture about how these two factors influence the accuracy of r and the 11 robust correlation estimates, this study also includes (5%, 5-SD), (5%, 10-SD), (10%, 5-SD), (20%, 5-SD), and (20%, 10-SD), for a total of six different metric-based conditions.

The method used in the current simulation is based on the fact that a 5%, 10% or 20% random sample of the original x and/or y scores are multiplied by 5 or 10. Given that there is an equal chance that a portion (5%, 10%, or 20%) of both positive and negative x and/or y scores can be manipulated as outliers/leverage points, and x and y follow symmetric distributions, a good correlation estimate should be robust to these outliers and/or leverage points in a sample, and provide an accurate estimate of the true linear association between X and Y ( $\theta$). Given these metric-based conditions, there are four different cases of O-LO discussed below.

Factor 4: Outliers and Leverage Observations (Distribution)

Four cases are evaluated. For Case A (uniform y-outliers), the y-coordinate is such that the point is an outlier, and the point does not fit to the linear trend of most other points. For each level of Factor 3, $\zeta$% of the data points were randomly selected and their y scores were multiple by $\varrho$. For Case B where the point is a leverage observation, $\zeta$% of the data points were randomly selected and both x and y in the same x-y pair are multiplied by $\varrho$. For Case C in which the x-coordinate of the point is an outlier and it does not fit to the linear trend of most other points, $\zeta$% of the data points were randomly selected and their x scores were multiple by $\varrho$. For Case D where both the x-coordinate and y-coordinate of the point are outliers and the x and y scores of the outlier points are correlated in the opposite direction of $\theta$, $\zeta$% of the x-y points were generated from Equation (19) where $\theta$ was replaced with - $\theta$. For each of these data points both the x and y in the same x-y pair were multiplied by $\varrho$ to become O-LO points.

In sum, the four factors were combined to produce a design with $4\times 3\times 6\times 4=288$ conditions that contain outliers and/or leverage points (O-LO). In addition, a total of $4\times 3=12$ ideal conditions in which x and y follow a normal distribution without O-LO are also evaluated. For each of the 300 conditions, 1,000 datasets were generated by first generating X scores as a random sample from a normal distribution with mean 0 and SD = 1 and e scores as a random sample from a normal distribution with mean 0 and SD = $\sqrt{1-{\theta }^2}$. Corresponding Y scores were generated through Equation (19) such that the expected linear relationship is $\theta$ (Chan & Chan, 2004). Each of these 1,000 datasets was resampled B = 2,000 times to construct the 95% BSI, BPI, and BCaI surrounding each of the 12 correlation estimates. The code was written and executed in the R Project programming environment (R Core Team, 2017), and is available in the Supplementary Materials (Part A).

Normality

As predicted, r produced highly accurate results with normal data (see Table 1). The biases ranged from -.017 to .013 with a mean of .000 [range-of-biases = (-.017, .013); mean-of-biases = .000]. All of the 12 conditions without OL-O yielded a bias within $\pm$.10. MAPE was also appropriate (.007). The remaining 11 robust correlation estimates were reasonable: 10 of them ( $r_{CME}$, $r_{MAD}$, $r_{MED}$, $r_{bm}$, $r_Q$, $r_{\tau }$, $r_s$, $r_{wLS}$, $r_{TRIM}$, and, $r_{pbc}$) produced a mean bias within $\pm$.10 and a MAPE smaller than .10. The only exception is $r_w$ with a mean bias of -.108 and a MAPE of .108. Hence, the majority of the robust correlations can reasonably estimate the level of linearity $\theta$ when the assumption of normality and no O-LO are present in the data.

Case A (Y-Outliers)

The performance of r substantially decreased when y-outliers were present in the data, with range-of-biases = (-.495, -.124) and mean-of-biases = -.305, showing a substantial downward bias. None of the 72 conditions produced an acceptable bias. Accordingly, MAPE was also less than optimal (.305). The 11 robust correlations noticeably outperformed r. Ten of them ( $r_{CME}$, $r_{MAD}$, $r_{MED}$, $r_{bm}$, $r_Q$, $r_{\tau }$, $r_s$, $r_{wLS}$, $r_{TRIM}$, and, $r_{pbc}$) produced a mean bias within $\pm$.10 and a MAPE smaller than .10, whereas $r_w$ resulted in a slightly larger mean bias of -.122 and a MAPE of .122. Comparatively, $r_Q$ appears to be the most accurate with a mean bias of -.018 and a MAPE of .026 under Case A. However, $r_Q$ did not perform well under Cases C and D. Rather, $r_{MAD}$, $r_{MED}$, and $r_{TRIM}$ performed best in these cases, and hence, they are further discussed in the following sections.

Case B (LO)

As predicted, r is reasonably robust to LO. The biases ranged from -.159 to .033 with a mean of -.027. Of the 72 conditions, 69 (or 95.8%) produced a bias within $\pm$.10, indicating a reasonable fit. MAPE was also acceptable (.029). Of the 11 robust correlations, only 7 ( $r_{MAD}$, $r_{MED}$, $r_Q$, $r_{\tau }$, $r_s$, $r_{TRIM}$, and, $r_{pbc}$) produced a bias within $\pm$.10. As noted above, $r_{MAD}$, $r_{MED}$, and $r_{TRIM}$ also behave appropriately under Cases C and D, and hence, they are further discussed. For $r_{MAD}$, the range-of-biases (-.171, .104) and the mean-of-biases was 027. Of the 72 conditions, 69 (or 95.8%) produced a bias within $\pm$.10. MAPE was also good (.032). For $r_{MED}$, the range-of-biases was (-.212, .041) and mean-of-biases was -.040. Of the 72 conditions, 69 (or 95.8%) yielded a bias within $\pm$.10. MAPE is also appropriate (.042). For $r_{TRIM}$, the mean-of-biases was (-.202, .045) and mean-of-biases was -.049. Of the 72 conditions, 67 (93.1%) produced a bias within $\pm$.10. MAPE is also acceptable (.051).

Case C (X-Outliers)

As in Case A, r proved to be a poorly performing estimate of the bivariate linear relationships under Case C conditions. The range-of-biases was (-.987, -.389) and mean-of-biases was -.776, showing a severe-downward bias. Of the 72 conditions, 0 (or 0%) produced an acceptable bias. MAPE was also very poor (.776). Regarding the 11 robust correlations, only 3 ( $r_{MAD}$, $r_{MED}$, and $r_{TRIM}$) yielded an acceptable bias within $\pm$.10. For $r_{MAD}$, the range-of-biases was (-.187, .055) and mean-of-biases was -.047. Of the 72 conditions, 66 (or 91.7%) produced a bias within $\pm$.10. MAPE was also appropriate (.050). For $r_{MED}$, the range-of-biases was (-.213, .036) and mean-of-biases was -.063. Of the 72 conditions, 61 (or 84.7%) yielded a bias within $\pm$.10. MAPE was also reasonable (.064). For $r_{TRIM}$, the range-of-biases was (-.176, .025) and mean-of-biases was -.069. Of the 72 conditions, 57 (79.2%) produced a bias within $\pm$.10. MAPE was also reasonable (.070).

Case D (LO with Opposite Correlation)

Unsurprisingly, Case D was found to be the data condition most detrimental to the accuracy of correlation estimates. For r, the range-of-biases was (-1.984, -.789) and mean-of-biases was -1.536. None of the 72 conditions produced an acceptable bias. MAPE was also high (1.536) indicating a poor fit. Comparing the 11 robust correlations, $r_{MAD}$, $r_{MED}$, and $r_{TRIM}$ produced reasonable results. For $r_{MAD}$, the range-of-biases was (-.210, -.024) and mean-of-biases was -.071. Of the 72 conditions, 55 (or 76.4%) produced a bias within $\pm$.10. MAPE was also reasonable (.072). For $r_{MED}$, the range-of-biases was (-.225, .003) and mean-of-biases was -.091. Of the 72 conditions, 47 (or 65.3%) yielded a bias within $\pm$.10. MAPE was also reasonable (.091). For $r_{TRIM}$, the range-of-biases was (-.226, .014) and mean-of-biases was -.094. Of the 72 conditions, 44 (61.1%) produced a bias within $\pm$.10. MAPE was also reasonable (.094).

Normality

The 3 bootstrap CIs constructed for r performed acceptably, with a mean CP of .937 for BSI, .942 for BPI, and .944 for BCaI (Table 2). Of the 12 conditions, 9 (or 75%) yielded an acceptable CP for BSI, and 11 (or 91.7%) resulted in an acceptable CP for both BPI and BCaI. The remaining 11 robust correlations also produced good results, in general, with mean CPs ranging from .927 to .968. The bootstrap CIs constructed for the three key robust correlations, $r_{MAD}$, $r_{MED}$, and $r_{TRIM}$, were appropriate. For $r_{MAD}$, the BSI and BCaI were appropriate, with a mean CP of .968 and .943, respectively. For $r_{MED}$, the mean CPs were .965 and .932 for BSI and BCaI, respectively, which are appropriate. For $r_{TRIM}$, the mean coverage probabilities were .958, .967, and .939, respectively, for the BSI, BPI, and BCaI, which are also good.

Case A (Y-Outliers)

Given that r is not robust to Case A, the associated bootstrap CIs are also less than optimal. The mean CPs are only .583, .616, and .50 for BSI, BPI, and BCaI, respectively. Of the 72 conditions, only 1 (or 1.4%), 10 (or 13.9%), and 4 (or 5.6%) conditions based on BSI, BPI, and BCaI, respectively, yielded a CP within (.925, .975), indicating a poor fit. Comparatively, the 11 robust correlations yielded much more reasonable CPs, ranging from .886 to .967. Regarding the bootstrap CIs for $r_{MAD}$, one of the 3 robust correlation estimates as discussed above, the mean CPs were .967 and .931 for BSI, and BCaI respectively, showing acceptable fit. The BSI and BCaI for $r_{MED}$, was also good with mean CPs of .966 and .930, respectively. The CIs for the last of the acceptably performing robust correlations, $r_{TRIM}$, also performed well, with mean CPs of .960, .958, and .937 for BSI, BPI, and BCaI, respectively.

Case B (LO)

Although r was found to be reasonably robust to leverage observations, the associated BSI, BPI, and BCaI were less than optimal. The mean CPs were .856, .909, and .886 based on BSI, BPI, and BCaI, respectively. Of the 72 conditions, only 4.2% (BSI), 20.8% (BPI), and 19.4% (BCaI) of these conditions produced a CP within (.925, .975). Comparatively, the CPs yielded by BSI, BPI, and BCaI surrounding the 11 robust correlations (except r_wLS) were more reasonable. Specifically, the BSI, and BCaI constructed for $r_{MAD}$ resulted in good mean CPs of .966 and .941, respectively. Of the 72 conditions, 81.9% (BSI) and 95.8% (BCaI) yielded a CPs within (.925, .975), respectively. Regarding $r_{MAD}$, BSI and BCaI resulted in good mean CPs of .968 and .933, respectively. Of the 72 conditions, 76.4% and 95.5% of them, respectively, produced a CP within (.925, .975). For $r_{TRIM}$, the mean CPs are .964, .974, and .945 for BSI, BPI, and BCaI, respectively, indicating an excellent fit. Of the 72 conditions, 75%, 55.6%, and 98.6% of conditions yielded a CP within [.925, .975].

Case C (X-Outliers)

Similar to Case A, the bootstrap CIs surrounding r were not robust to x-outliers. The mean CPs were .588 for BSI, .641 for BPI, and .589 for BCaI, which are undesirable. Of the 72 conditions, 0%, 1.4%, and 1.4% resulted in a CP within (.925, .975) based on BSI, BPI, and BCaI, respectively. Comparatively, the bootstrap CIs constructed for the robust correlations (except r_wLS) were reasonable. Specifically, the BSI and BCaI constructed for $r_{MAD}$ (or $r_{MED}$) are good, with mean CPs of .967 and .938 (or .967 and .930) respectively. Of the 72 conditions, 77.8% and 97.2% (or 80.6% and 76.4%) of the BSIs and BCaIs, respectively, produced a CP within (.925, .975). Regarding $r_{TRIM}$, the mean CPs were .964, .968, and .943 for BSI, BPI, and BCaI, respectively. Of the 72 conditions, 76.4% (BSI), 69.4% (BPI), and 91.7% (BCaI) produced a CP within (.925, .975).

Case D (LO With Opposite Correlation)

r is not robust to Case D. The mean CPs were .359, .444, and .389 for BSI, BPI, and BCaI, respectively. Of the 72 conditions, 0% from BSI, BPI, or BCaI fell within the range of (.925, .975). Comparatively, the bootstrap CIs for $r_{MAD}$, $r_{MED}$, and $r_{TRIM}$ were appropriate. For $r_{MAD}$, the mean CPs were .966 and.932 for BSI and BCaI, respectively. Of the 72 conditions, 77.8% and 84.7% yielded a CP within (.925, .975). For $r_{MED}$, the mean CPs were .965 and .923 for BSI and BCaI, respectively. Of the 72 conditions, 76.4% and 63.9% yielded a CP within (.925, .975). The bootstrap CIs for $r_{TRIM}$ appear to have performed the best. The mean CPs were .959, .952, and .932 for BSI, BPI, and BCaI, respectively. Of the 72 conditions, 77.8%, 62.5%, and 83.3% yielded a CP within (.925, .975).