When to use CorSym
CorSym is a correlation estimator designed for exchangeable pair
data, where the two values you want the correlation of come in an
arbitrary/random order. This occurs frequently in the study of
assortative pairing, where we wish to desire if individuals that are
paired (meaning dating, co-habitation, marriage, or sharing children)
have more similar characteristics (for example, socioeconomics or
genetics) than random pairs in the population. More broadly, the
exchangeable variables that CorSym requires are essentially separate
measures of the same variable (with the same distribution) in two cases
(like two people, in assortative pairing). Two variables \(X\) and \(Y\) are exchangeable if the joint
probability distribution of \((X,Y)\)
and \((Y,X)\) are the same.
We developed this estimator because we found that the obvious way to
estimate correlations, using the Pearson estimator, is not only
sensitive to the order (the estimates change if order is changed) but
also statistically biased if the order is biased, meaning that the first
value tends to be smaller than the second, or viceversa. CorSym
estimates correlations for data like this without bias, and the answer
does not depend on the order of the data.
However, most applications requiring correlation estimates do not
satisfy these assumptions! For example, if you want to estimate the
correlation between age and smoking, you cannot randomly switch the age
and smoking status of one observation without doing it to all other
observations. In such a case, CorSym is not the right tool and its
estimates will be nonsense.
Simulated example: extreme ordering
Realistic data ideally is truly correlated, but here for simplicity
we only simulate uncorrelated data. You can easily generate correlated
data with the mvtnorm and other packages, as we do in our
paper, but here we don’t use those to reduce package dependencies.
Nevertheless, you will see how data with zero correlation will have
large Pearson estimates when the order of the data is biased. In
general, order biases result in upward Pearson correlation biases.
Let’s simulate simple independent (uncorrelated) data for these
examples. To be exchangeable, it is vital that the marginal
distributions of both columns be the same! (Here both are Standard
Normal.) Let’s simulate a relatively large number of pairs so estimates
don’t have much noise, and the biases are evident.
# number of pairs
n <- 1000
# matrix with different pairs in rows
X <- matrix( rnorm( n * 2 ), nrow = n, ncol = 2 )
The above data is in “random” order: although the proportion of pairs
where the value in the first column is less than that of the second
column is rarely exactly 0.5, it will differ from it by small random
amounts, without a preference (bias) for lower or higher values than
0.5.
To contrast, we create a version of this data with “extreme” order,
where every pair lists the minimum value first, and maximum value
last.
library(corsym)
Xe <- partial_order( X )
We can also take biased data and “randomize” it, for further
comparison:
Xr <- randomize_order( Xe )
Let’s quantify the order biases by counting the number of times a
given order was observed (say x<y where x
is the first column and y is the second).
data <- data.frame(
type = c('Random 1', 'Random 2', 'Extreme'),
x = c(
sum( X[,1] < X[,2] ),
sum( Xr[,1] < Xr[,2] ),
sum( Xe[,1] < Xe[,2] )
),
n = n
)
With real data it may not be clear if order is truly biased or if the
order proportions are very different than 0.5 by chance, which occurs
more often at small sample sizes. This function takes a table such as
the one above, and adds columns with the empirical order proportions and
the p-value from a two-sided Binomial test under the null hypothesis
that the proportion is 0.5:
order_bias_test( data )
#> type x n f pval
#> 1 Random 1 479 1000 0.479 1.947663e-01
#> 2 Random 2 500 1000 0.500 1.000000e+00
#> 3 Extreme 1000 1000 1.000 1.866527e-301
Since sample sizes are huge, you’ll see insignificant p-values for
the first two cases, and an extremely significant p-value in the last
one.
Now that we confirmed that the initial and randomized data did not
have an order bias, while the “extreme” data has the strongest bias
possible, let’s see the effect on the correlation estimates. Here is the
combinatorial experiment, two estimators times three orders:
# Each function returns 3 values:
# estimate, and lower and upper confidence intervals (CIs)
pearson( X )
#> [1] 0.03029735 -0.03175538 0.09211741
pearson( Xr )
#> [1] 0.02878629 -0.03326616 0.09061766
pearson( Xe )
#> [1] 0.4881381 0.4394431 0.5339725
corsym( X )
#> [1] 0.02917747 -0.03282849 0.09095968
corsym( Xr )
#> [1] 0.02917747 -0.03282849 0.09095968
corsym( Xe )
#> [1] 0.02917747 -0.03282849 0.09095968
When the data is in random order, Pearson and CorSym are both
unbiased and agree approximately. They do not give exactly the same
estimate, but they will be close to each other and to zero (because the
data was uncorrelated), and their CIs should overlap zero and each other
most of the time. Also note that Pearson gives slightly different
estimates for each random order. However, when the data has the extreme
order, Pearson gives a much bigger answer, significantly different than
zero! In contrast, CorSym gives exactly the same answer for all
orders.
Simulated example: partial orderings
Since CorSym always gives the same answer for all orders, we will
ignore it for now. Let’s focus on establishing the effect of partial
data ordering on Pearson estimates.
# sequence of forced order proportions between 0 and 1
# (q=0 is original data, q=1 is extreme order)
qs <- ( 0:10 ) / 10
# data frame to grow
data <- NULL
for ( q in qs ) {
# force a proportion `q` of rows to have extreme order
Xp <- partial_order( X, q )
# pearson estimates
out <- pearson( Xp )
# gather into data frame for plotting
data_q <- data.frame(
n = n,
q = q,
# counts for biased order test
x = sum( Xp[,1] < Xp[,2] ),
# Pearson estimates
r = out[1],
CIL = out[2],
CIU = out[3]
)
data <- rbind( data, data_q )
}
# perform order bias tests on these data
data <- order_bias_test( data )
# inspect table
data
#> n q x r CIL CIU f pval
#> 1 1000 0.0 479 0.03029735 -0.031755379 0.09211741 0.479 1.947663e-01
#> 2 1000 0.1 527 0.02885815 -0.033194316 0.09068899 0.527 9.368729e-02
#> 3 1000 0.2 592 0.03815082 -0.023898787 0.09990761 0.592 6.520547e-09
#> 4 1000 0.3 637 0.05323657 -0.008785503 0.11485061 0.637 3.708563e-18
#> 5 1000 0.4 691 0.07578979 0.013861832 0.13713853 0.691 3.063375e-34
#> 6 1000 0.5 749 0.11842491 0.056849188 0.17910310 0.749 4.037455e-58
#> 7 1000 0.6 781 0.13934254 0.078023445 0.19961132 0.781 1.512936e-74
#> 8 1000 0.7 847 0.19808971 0.137788700 0.25692767 0.847 5.328759e-117
#> 9 1000 0.8 892 0.25097069 0.191964270 0.30816914 0.892 3.984665e-154
#> 10 1000 0.9 950 0.37634334 0.321859457 0.42834289 0.950 1.863693e-216
#> 11 1000 1.0 1000 0.48813810 0.439443084 0.53397254 1.000 1.866527e-301
We can establish a lot of facts about this data now. A simple one
first, the “forced order proportion” \(q\) is not equal to the order proportion
\(f\), because random data (\(q=0\)) has an expected order proportion of
\(\text{E}[f] = 0.5\). Thus, when \(q\) fraction of the data is forced to have
extreme order, the expected order proportion becomes \(\text{E}[f] = ( 1 - q ) ( 0.5 ) + q ( 1 ) = 0.5(1
+ q)\). The data shows agreement with this formula:
library(ggplot2)
ggplot( data, aes( x = q, y = f ) ) +
geom_point() +
theme_classic() +
geom_abline( intercept = 0.5, slope = 0.5, linetype = 'dashed', color = 'gray' )
Now we can see the relationship between the forced order proportion
\(q\) and the bias. Since the true
correlation is zero, the estimate \(r\)
is the bias in this case. The data shows a non-linear relationship,
where \(q \le 0.5\) CIs are likely to
contain the true value of zero, while \(q >
0.5\) cases are more likely to be significantly non-zero and
quickly increase.
ggplot( data, aes( x = q, y = r ) ) +
geom_point() +
theme_classic() +
geom_errorbar( aes( ymin = CIL, ymax = CIU ) ) +
geom_hline( yintercept = 0, linetype = 'dashed', color = 'gray' )
Lastly, we can look directly at significance of the order bias. We
have a lot of power to detect differences, so even at \(q=0.1\) we tend to get significant results.
Comparing to the previous figure, we see that a significant order bias
does not directly result in significantly biased Pearson estimates.
# Bonferroni threshold
pcut <- 0.05 / nrow( data )
pcut
#> [1] 0.004545455
ggplot( data, aes( x = q, y = -log10( pval ) ) ) +
geom_point() +
theme_classic() +
geom_hline( yintercept = -log10( pcut ), linetype = 'dashed', color = 'gray' )
