Correlation between two vectors of binary values

Question

This may simply be a problem of me not knowing how to describe my issue properly, if so, apologies in advance for the duplicate.

Some context: I am studying the distribution of genetic mutations in a population. My data is organised like so:

mutationID | Patient 1 | Patient 2 | ... | Patient N
  mut01    |  0        |  0        | ... |  1
  mut02    |  1        |  0        | ... |  1
  mut03    |  0        |  0        | ... |  0
  mut04    |  1        |  1        | ... |  1

Where 0 indicates 'This patient does not have this mutation', and 1 indicates 'This patient has this mutation'.

I am looking to measures 'correlation' between mutations, not between patients. So these two mutations would show up as having a strong 'score':

mutA: 0,0,1,1,0,0,1,1,0,0
mutB: 0,0,1,1,0,0,1,1,0,1

Originally I had tried Pearson's R, but I soon realised that it's probably not the best option here. Currently I am using Jaccard similarity index, but it does not capture all the information that I would like it to. For example:

mutC: 0,0,0,1,1,1,0,0,0,1
mutD: 1,1,1,0,0,0,1,1,1,1

For these two mutations, the Jaccard similarity index would be next to 0, seeing as they share few occurrences. However, they are strongly anticorrelated, and that is some information that I would like to capture.

What would be the best measure of correlation to use here?

Answer 1

Short answer: You can simply substitute zeros by $-1$s, multiply the two vectors and divide by the vector length. Then a vector of all ones would be perfectly negatively "correlated" with a vector of minus ones.

Long answer: You are on the right path with the Jaccard index, because your data are not naturally numeric. You have coded them as numbers, but they are essentially nominal (Boolean, to be more precise). Correlation is probably not the right measure here.

Jaccard index, however, has the drawback that it only takes into account, as the normalising factor, the patients having at least one mutation. Patients with no mutation are ignored, which is also a loss of information.

What you might do instead is to calculate the Simple Matching Coefficient (SMC), which is one minus the Hamming distance between two vectors (the number of positions they differ) divided by the vector length. It is still not quite correlation-like, because two "anticorrelated" vectors would have a SMC of 0.

Luckily, it is trivially easy to transform SMC into something looking like "correlation" for binary vectors: multiply by two and subtract one. This is essentially what you get by the above short "algorithm". You can, of course, do the same with Jaccard index, but the results would differ due to different denominators.

For your first example

mutA: 0,0,1,1,0,0,1,1,0,0
mutB: 0,0,1,1,0,0,1,1,0,1

Jaccard and SMC produce 0.8 and 0.9, respectively, which translate into "correlations" of 0.6 and 0.8.

On the other example:

mutC: 0,0,0,1,1,1,0,0,0,1
mutD: 1,1,1,0,0,0,1,1,1,1

the Jaccard index and SMC would be identical, 0.1 and the "correlation" -0.8.

The difference between the two measures is more striking for the following example (mine):

mutX: 0,0,0,0,0,0,0,0,0,1
mutY: 0,0,0,0,0,0,0,0,1,1

Jaccard: 0.5 ("correlation" = 0) vs. SMC: 0.9 ("correlation" = 0.8). I believe the SMC-based "correlation" better captures the relationship you're after. It is identical to your first example and symmetrical to your second one:

• all bits except one are same $\Rightarrow$ correlation = +0.8
• all bits except one are different $\Rightarrow$ correlation = -0.8

Answer 2

In order to think about what is the "best" measure of association here, it is worth stepping back and looking at what you are trying to summarise. Since you have $N$ patients, each mutation vector is a binary vector of length $N$. If we assume that data for different patients is IID for each mutation, then the data can be summarised without loss of information by a $2 \times 2$ contingency table, which we can denote in matrix form as:

$$\mathbf{m} = \begin{bmatrix} m_{00} & m_{01} \\ m_{10} & m_{11} \end{bmatrix}.$$

The four values in the contingency tables give counts of the pairs of values between two mutations; the value $m_{ij}$ is the count of the number of pairs where the first mutation is equal to $i$ and the second is equal to $j$ (with these being binary indicators).

Now, a full summary of the underlying distribution would give the probabilities of each of the four outcomes in the table, which is three free parameters. By reducing to a measure of association you are reducing this to a single parameter. Any measure of association should have some basic monotonicity properties and give a fixed outcome value for independence. However, there are a range of possible measures that use different scales and have different rates of increase/decrease as you get more positive/negative association.

You can find a comprehensive examination of various alternative measures of association for binary data in Duan et al (2014). In that paper the authors examine several measures of association to see which desirable properties they have, and they conduct analysis using simulated data and real data. In particular, this gives you some information on the properties of different measures. As you can see from that paper, even the standard Pearson correlation coefficient has reasonble properties for binary data, including a basic symmetry property that is desirable when looking at association of values in a $2 \times 2$ contingency table.

Answer 3

I'll argue that the simple covariance can be enough, if we know how to interpret it.

---

Since mutations are binary, we can treat each within-subject mutation as Bernoulli-distributed.

Given $X \sim \mathcal B (p) $ and $Y \sim \mathcal B (q)$, we can calculate the covariance:

$$\operatorname{Cov}(X,Y) = E((X-E(X))(Y-E(Y))=E(XY)-E(X)E(Y)\\=\frac{1}{n}\sum(x_iy_i)-\frac{1}{n}\sum(x_i)\frac{1}{n}\sum(y_i)=\\k- \hat p\hat q$$

Where $k = \frac{1}{n}\sum(x_iy_i)$, i.e. the proportion of times where both x & y co-occur. This number is non-normalized, below I'll give hints towards interpreting it.

Notice that this has a very direct interpretation already: if the covariance is positive, then the co-occurrences are more frequent than expected by chance, and if its negative, then they are less frequent than expected by chance. In fact: $ \max(\hat p + \hat q - 1, 0)\leq k \leq \min(\hat p, \hat q)$ & $0 \leq\hat p\hat q\leq 1$. So this covariance necessarily has a tighter bound than $(-1,1)$, a nice property to have.

The only thing missing is the lack of treatment for the rates of occurrences, $\hat p$ and $\hat q$. They appear naturally in the covariance matrix though.

---

We could follow on from here to the sample correlation, noting that

$$\operatorname{Var}(X) = E(X^2)-E(X)^2=E(X)-E(X)^2=\hat p-\hat p^2=\hat p(1-\hat p)\\ \operatorname{Var}(Y) = \hat q(1-\hat q)$$

Then

$$\operatorname{Cor}(X,Y)=\frac{k - \hat p\hat q}{\hat p(1-\hat p)\hat q(1-\hat q)} = \left(\frac{k}{\hat p\hat q}-1\right)\frac{1}{(1-\hat p)(1-\hat q)}$$

The correlation is obviously bounded to $(-1,1)$, and it's a good idea to note when that occurs.

$$\begin{cases} \begin{aligned} &Y = 1 - X &\rightarrow \quad &\operatorname{Cor}(X,Y)= -1 \\ &k < \hat p \hat q &\rightarrow \quad &\operatorname{Cor}(X,Y)<0\\ &k = \hat p \hat q &\rightarrow \quad &\operatorname{Cor}(X,Y)=0 \quad \text{(chance)}\\ &k > \hat p \hat q &\rightarrow \quad &\operatorname{Cor}(X,Y)>0\\ &Y = X &\rightarrow \quad &\operatorname{Cor}(X,Y)= 1 \end{aligned} \end{cases}$$