How to weight a Spearman rank correlation by statistical errors?

Question

I'm trying to evaluate whether two quantities, X and Y, are correlated or not. I have a sample of N items, for which I have measured X and Y, both with measurement errors X_err and Y_err. X and Y are not thought to follow a normal distribution, so Spearman's rank is preferential. However, as far as I understand, Spearman's rank is not designed to take error bars into account - or is there a way to weight the test? Or should I be using a completely different test? If it can be done with the scipy package, that would be great.

As an example to illustrate when the non-weighted Spearman's rank does not behave like I would prefer, here are two sets of made-up data:

This one has a correlation coefficient of 1, and a p-value of 0, i.e. a perfect correlation, because the values are monotonously increasing.

This one has a correlation coefficient of 0.9 and a p-value of 0.03, so definitely a significant correlation, but worse than the first.

So obviously I wouldn't want the test to tell me that the first one shows any significant correlation - since the errors are so big it was just pure luck that the values ended up being completely monotonously increasing. The second one is, on the other hand, a pretty good correlation. What test can account for this?

Answer 1

This paper might help you. Here its abstract:

This manuscript describes a number of easily implemented, Monte Carlo based methods to estimate the uncertainty on the Spearman’s rank correlation coefficient, or more precisely to estimate its probability distribution.

Basically, the idea is the following:

1. Simulate many samples from the original data, using the "error bars" in your data ($X_{err}$) to introduce noise to the samples.
2. Then, for each sample, take the Spearman correlation.
3. You get so many Spearman "values" as the number of samples you took on the first step.
4. Use these many points to calculate a "distribution" of Spearman, rather than a point estimate.

PS: @Cadnr has given the same advice in a previous answer.

Answer 2

I hate to be another brainless advocate of Monte Carlo methods, but one solution would be to build up a distribution of p values by taking a large number of samples of your data error distributions. For each data point, generate random errors in x and y (within the envelope defined by the measurement errors for that data point), and once that's been done for all the data points, generate the p value for the synthetic dataset. Hopefully, as you repeat this process many times, your distribution of synthetic p values will approach a well-defined functional form (such as a Gaussian), for which you can find the median and useful limits by doing a (Gaussian) fit or taking e.g. the median and +-67% percentiles. You'll then end up with a p value and +- errors, from which you'll be able to tell if the correlation is significant.

I'm not aware of any off-the-shelf software to help you accomplish this, but it shouldn't be hard to code.

Answer 3

You can construct a Spearman-like correlation that takes into account weights.

Let's say, we have two rankings, $Q$ and $R$ and two set of weights $W_q$ and $W_r$ (you can have one of these be all ones if you have only one set of weights). You would have to compute these from your errors. All of these have $n$ elements.

Now we want to compute the weighted rank correlation in a Spearman-like manner, i.e. the correlation coefficiont should be a function of $\sum_i^n w^Q_iw^R_iD_i^2=\sum_i^n w^Q_iw^R_i(R_i-Q_i)^2$.

We assume that this function will have the form $A+B\sum_i^nw^Q_iw^R_iD_i^2$, just like the original Spearman-function.

Since we want it to be equal to one if the rank orders are identical, i.e. $D_i = 0$ for all $i$, it follows that $A=1$.

If the rankings are reversed, we want the result to be $-1$. In this case the $D_i$ will be equal to $n-1, n-3,...,-(n-1)=n-2i+1$.

This allows us to compute $B$: $$ B\sum_i^nw^Q_iw^R_i(n-2i+1)^2 = -2 \\ B = \frac{-2}{\sum_i^nw^Q_iw^R_i(n-2i+1)^2} $$

Our function is therefore: $$ r(R,Q,W_r,W_q) = 1 - \frac{2\sum_i^nw^Q_iw^R_i(R_i-Q_i)^2}{\sum_i^nw^Q_iw^R_i(n-2i+1)^2} $$

I'm not 100% sure about this, but I believe your weights must be nonnegative, because otherwise your function might (as I said, I'm not completely sure) leave $[-1,1]$.

Also, of course, your weights must be on the same scale.

Also, for your the function to behave correctly, the weights need to be monotonous with increasing rank.

EDIT: I just realized that this answer probably won't help in your case, since the type of weights you are looking for are something different than what I wrote my function for. Sorry. I'll leave this answer here if it helps anybody else though.