---
title: "CorSym: Correlation Estimation For Exchangeable/Symmetrical Variables"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{CorSym: Correlation Estimation For Exchangeable/Symmetrical Variables}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

# 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.
```{r}
# 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.
```{r}
library(corsym)
Xe <- partial_order( X )
```
We can also take biased data and "randomize" it, for further comparison:
```{r}
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).
```{r}
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:
```{r}
order_bias_test( data )
```
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:
```{r}
# Each function returns 3 values:
# estimate, and lower and upper confidence intervals (CIs)
pearson( X )
pearson( Xr )
pearson( Xe )
corsym( X )
corsym( Xr )
corsym( Xe )
```
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.
```{r}
# 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
```

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:
```{r}
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.
```{r}
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.
```{r}
# Bonferroni threshold
pcut <- 0.05 / nrow( data )
pcut

ggplot( data, aes( x = q, y = -log10( pval ) ) ) +
	geom_point() +
	theme_classic() +
	geom_hline( yintercept = -log10( pcut ), linetype = 'dashed', color = 'gray' )
```