# Testing whether the conditional correlations/covariances differ between two groups

I have two samples of variables $$\{y_{1i},y_{2i},x_i,s_i\}$$. Where $$y_1$$ and $$y_2$$ are binary variables, $$x$$ is a continuous variable and $$s$$ is a sample indicator, taking the value 0 in one sample and the value 1 in the other.

I can test whether means are different for $$y_{1i}$$ and $$y_{2i}$$ separately using a simple t-test or equivalently by estimating this linear equation using OLS$$y_i=\alpha+\beta x_i+\tau s_i + \varepsilon_i$$The regression has the advantage of providing a test of differences in the conditional means $$E[y_i1|x_i,s_i=0]-E[y_i1|x_i,s_i=1]$$

Now I also want to test whether the relationship between $$y_1$$ and $$y_2$$ is the same in both samples. For example I want to test whether:$$Corr[y_{1i},y_{2i}|x_i,s_i=0]=Corr[y_{1i},y_{2i}|x_i,s_i=1]$$ To me this seems like it should be a textbook statistical test, but I can't seem to find what the default procedure/test for this would be.

• Similar to stats.stackexchange.com/questions/372385/…, but simpler. Oct 18 '18 at 4:33
• I see the similarity, though my intuition was that there should somehow be a very simple test for this, after all, I am just trying to compare one moment of a bi-variate distribution across two samples. Another difference is that the linked question is looking for a test in a paired sample, which is different in my case.
– sheß
Oct 18 '18 at 17:40
• Well, since the y's are binary, you can just present the data as two contingency tables for each sample. Then you could estimate/compare the log odds ratio for the tables. That can be calculated via a logistic regression Oct 18 '18 at 22:05
• Thanks for the suggestion. I am just not sure I understand which log odds ratio and how that helps me test the difference in the correlation. Keep in mind that I already know that there will be an effect on the means of y1 and y2. What I am interested in is figuring out if and how their relationship changed
– sheß
Oct 18 '18 at 22:22

Since you have two binary response variables in your analysis, there are four possible outcomes taken over those two variables. Hence, it makes sense to use some kind of regression model with a discrete four-category response ---e.g., the multinomial logistic regression model. The model form for this model can be written as:

\begin{aligned} \mathbb{P}(Y_{1i}=0,Y_{2i}=0 | x_i, s_i) &= \frac{1}{1 + e^{\beta_{01} x_i + \gamma_{01} s_i}+e^{\beta_{10} x_i + \gamma_{10} s_i}+e^{\beta_{11} x_i + \gamma_{11} s_i}}, \\[10pt] \mathbb{P}(Y_{1i}=0,Y_{2i}=1 | x_i, s_i) &= \frac{e^{\beta_{01} x_i + \gamma_{01} s_i}}{1 + e^{\beta_{01} x_i + \gamma_{01} s_i}+e^{\beta_{10} x_i + \gamma_{10} s_i}+e^{\beta_{11} x_i + \gamma_{11} s_i}}, \\[10pt] \mathbb{P}(Y_{1i}=1,Y_{2i}=0 | x_i, s_i) &= \frac{e^{\beta_{10} x_i + \gamma_{10} s_i}}{1 + e^{\beta_{01} x_i + \gamma_{01} s_i}+e^{\beta_{10} x_i + \gamma_{10} s_i}+e^{\beta_{11} x_i + \gamma_{11} s_i}}, \\[10pt] \mathbb{P}(Y_{1i}=1,Y_{2i}=1 | x_i, s_i) &= \frac{e^{\beta_{11} x_i + \gamma_{11} s_i}}{1 + e^{\beta_{01} x_i + \gamma_{01} s_i}+e^{\beta_{10} x_i + \gamma_{10} s_i}+e^{\beta_{11} x_i + \gamma_{11} s_i}}. \\[10pt] \end{aligned}

In this model form you have unknown coefficients $$\beta_{01}, \beta_{10}, \beta_{11}$$ that describe the effect of $$x_i$$ and you have coefficients $$\gamma_{01}, \gamma_{10}, \gamma_{11}$$ that describe the effect of $$s_i$$. This gives you a basic model you can fit to your data. If you are interested in testing a hypothesis relating to the correlation between the response variables then you will need to frame this hypothesis in terms of the underlying parameters of the model. To do this you will need to use some algebra to derive the relevant correlation expressions in terms of these parameters.

Finding the correlation formulae: With a bit of algebra, we can show that the true covariance between the two response variables is:

\begin{aligned} \mathbb{Cov}(Y_{1i}, Y_{2i} | x_i,s_i) &= \mathbb{E}(Y_{1i} \cdot Y_{2i} | x_i,s_i) - \mathbb{E}(Y_{1i} | x_i,s_i) \cdot \mathbb{E}(Y_{1i} | x_i,s_i) \\[10pt] &= \mathbb{P}(Y_{1i} =1, Y_{2i}=1 | x_i,s_i) - \mathbb{P}(Y_{1i} =1 | x_i,s_i) \cdot \mathbb{P}(Y_{1i}=1 | x_i,s_i) \\[10pt] &= \frac{e^{\beta_{11} x_i + \gamma_{11} s_i} (1- e^{\beta_{11} x_i + \gamma_{11} s_i})}{(1 + e^{\beta_{01} x_i + \gamma_{01} s_i}+e^{\beta_{10} x_i + \gamma_{10} s_i}+e^{\beta_{11} x_i + \gamma_{11} s_i})^2}. \\[10pt] \end{aligned}

For the variances, with a bit more algebra we get:

\begin{aligned} \mathbb{V}(Y_{1i} | x_i,s_i) &= \frac{(e^{\beta_{10} x_i + \gamma_{10} s_i} + e^{\beta_{11} x_i + \gamma_{11} s_i} )(1- e^{\beta_{10} x_i + \gamma_{10} s_i}- e^{\beta_{11} x_i + \gamma_{11} s_i})}{(1 + e^{\beta_{01} x_i + \gamma_{01} s_i}+e^{\beta_{10} x_i + \gamma_{10} s_i}+e^{\beta_{11} x_i + \gamma_{11} s_i})^2}, \\[10pt] \mathbb{V}(Y_{2i} | x_i,s_i) &= \frac{(e^{\beta_{01} x_i + \gamma_{01} s_i} + e^{\beta_{11} x_i + \gamma_{11} s_i} )(1- e^{\beta_{01} x_i + \gamma_{01} s_i}- e^{\beta_{11} x_i + \gamma_{11} s_i})}{(1 + e^{\beta_{01} x_i + \gamma_{01} s_i}+e^{\beta_{10} x_i + \gamma_{10} s_i}+e^{\beta_{11} x_i + \gamma_{11} s_i})^2}. \\[10pt] \end{aligned}

Hence, we obtain the correlation formula:

\begin{aligned} \mathbb{Corr}(Y_{1i},Y_{2i} | x_i,s_i) &= \frac{(e^{\beta_{10} x_i + \gamma_{10} s_i} + e^{\beta_{11} x_i + \gamma_{11} s_i} )(1- e^{\beta_{10} x_i + \gamma_{10} s_i}- e^{\beta_{11} x_i + \gamma_{11} s_i})}{\sqrt{(e^{\beta_{10} x_i + \gamma_{10} s_i} + e^{\beta_{11} x_i + \gamma_{11} s_i} )(1- e^{\beta_{10} x_i + \gamma_{10} s_i}- e^{\beta_{11} x_i + \gamma_{11} s_i}) \\ \cdot (e^{\beta_{01} x_i + \gamma_{01} s_i} + e^{\beta_{11} x_i + \gamma_{11} s_i} )(1- e^{\beta_{01} x_i + \gamma_{01} s_i}- e^{\beta_{11} x_i + \gamma_{11} s_i})}}. \\[10pt] \end{aligned}

So the particular correlation formulae for $$s_i=0$$ and $$s_i=1$$ are:

\begin{aligned} \mathbb{Corr}(Y_{1i},Y_{2i} | x_i,s_i=0) &= \frac{(e^{\beta_{10} x_i} + e^{\beta_{11} x_i})(1- e^{\beta_{10} x_i}- e^{\beta_{11} x_i})}{\sqrt{(e^{\beta_{10} x_i} + e^{\beta_{11} x_i} )(1- e^{\beta_{10} x_i}- e^{\beta_{11} x_i}) (e^{\beta_{01} x_i} + e^{\beta_{11} x_i})(1- e^{\beta_{01} x_i}- e^{\beta_{11} x_i})}}, \\[10pt] \mathbb{Corr}(Y_{1i},Y_{2i} | x_i,s_i=1) &= \frac{(e^{\beta_{10} x_i + \gamma_{10}} + e^{\beta_{11} x_i + \gamma_{11}} )(1- e^{\beta_{10} x_i + \gamma_{10}}- e^{\beta_{11} x_i + \gamma_{11}})}{\sqrt{(e^{\beta_{10} x_i + \gamma_{10}} + e^{\beta_{11} x_i + \gamma_{11}} )(1- e^{\beta_{10} x_i + \gamma_{10}}- e^{\beta_{11} x_i + \gamma_{11}}) \\ \cdot (e^{\beta_{01} x_i + \gamma_{01}} + e^{\beta_{11} x_i + \gamma_{11}} )(1- e^{\beta_{01} x_i + \gamma_{01}}- e^{\beta_{11} x_i + \gamma_{11}})}}. \\[10pt] \end{aligned}

The requirement that $$\mathbb{Corr}(Y_{1i},Y_{2i} | x_i,s_i=0) = \mathbb{Corr}(Y_{1i},Y_{2i} | x_i,s_i=1)$$ means that you equate these two formulae and find the conditions on the parameters that lead to this equation. (Algebraically, it is easier to equate the squared-correlations, which simplifies the expression.) Since you want this equation to hold for all $$x_i$$ this will impose quite strict requirements on the parameters.

A sufficient condition for equivalence is to have $$\gamma_{01} = \gamma_{10} = \gamma_{11} = 0$$ so that the variable $$s_i$$ has no effect at all on the responses. (In other words, it is not really a distinct test in that case.) It is not obvious whether there is any weaker condition that would yield the required correlation equivalence.

• Seems to me like this reduces to $\gamma_{ij}=0$ Oct 23 '18 at 5:04
• That might well be the case - I'll leave open the possibility of a weaker condition for now, but it might be possible to show that that is a necessary condition.
– Ben
Oct 23 '18 at 5:13
• This is very useful thank you. Slightly off topic: Is there a direct way this extends to the more general case where $y_1$ and $y_2$ are not binary but continuous?
– sheß
Oct 23 '18 at 9:13
• If they are both continuous then some kind of multivariate regression would be a reasonable place to start.
– Ben
Oct 23 '18 at 21:55
• Thanks. @Ben, the answer I came up with (and posted below) seems a lot simpler and is not limited to the binary case.. Since, if something sounds to good to be true it probably isn't, can you tell me what's wrong with my answer?
– sheß
Oct 25 '18 at 0:30

EDIT: This is my attempt assuming the outcomes are joinly normal, instead of binary

First, I phrase my assumptions about the joint distribution as: $$\left(\begin{matrix}y_1\\y_2\end{matrix}\right) \sim N\left(\mu,\Sigma\right),$$ with $$\mu=\left(\begin{matrix}\mu_1+\gamma_1s+\delta_1x\\\mu_2+\gamma_2s+\delta_2x\end{matrix}\right)$$ and $$\Sigma=\left(\begin{matrix}\sigma^2_1&,& \alpha+\beta_1 x + \beta_2 s\\\cdot&,& \sigma^2_2\end{matrix}\right)$$

Testing whether the covariance differs by sample is equivalent to testing that $$\beta_2=0$$. Do do that, I estimate $$Cov[y_{1i},y_{2i}|x_{i},s_i] =\\ E[(y_{1i}-E[y_{1i}|x_{i},s_i])\cdot(y_{2i}-E[y_{2i}|x_{i},s_i])|x_{i},s_i]=\alpha+\beta_1x_i+\beta_2s_i \quad\quad\quad\quad\quad\quad$$ by:

1. fitting a simple linear (or logit) model to estimate $$E[y_{1i}|x_{i},s_i]$$ and $$E[y_{2i}|x_{i},s_i]$$.
2. obtain the fitted values from this and compute the expression: $$\widehat{cov_i}=(y_{1i}-\widehat{E[y_{1i}|x_{i},s_i]})\cdot(y_{2i}-\widehat{E[y_{2i}|x_{i},s_i]})$$
3. estimate the model: $$\widehat{cov}_i=\alpha+\beta_1x_i+\beta_2s_i + \varepsilon_i$$

This seems to work well. If correlations are of interest instead of covariances, step 2. can be augmented dividing the $$\widehat{cov_i}$$ by similarly obtained estimators for the two standardeviations.

• @Ben, given that this seems so much simpler than you answer, I assume it is somehow wrong, or makes outrageously stupid assumptions. Can you tell me what this is?
– sheß
Oct 25 '18 at 0:28
• Well, the first issue with it is that you are fitting linear regression to a binary response, which is generally going to be a bad model. Secondly, even doing this, the first equation you give does not appear to be correct. (Indeed, it cannot possibly be correct, because it purports to give the covariance, but the right-hand-side depends on the error term, which is random.) If you want to estimate the covariance using a linear regression then you will need to derive the result correctly. That the moment what you have appears to me to be incorrect.
– Ben
Oct 25 '18 at 1:06
• Thanks a lot for your answer. Changed by post, as I realized that the way I wrote it was confusing. I also realized that my solution is for Normally distributed variables. But I think the version as it is now, addresses your concerns at least for the Normal-case. Do you agree?
– sheß
Oct 25 '18 at 1:50
• ..and if this is correct for bivariate normals, then the same logic should also apply to binary outcomes,no?
– sheß
Oct 25 '18 at 1:59
• You are getting closer, but there are still a lot of problems. Firstly, the specified covariance matrix in your regression will not be a valid covariance if $\alpha+\beta_1x+\beta_2s > \sigma_1 \sigma_2$. So at the moment your model is invalidly specified. Secondly, your covariance result appears to actually be the covariance of the error terms, not the response terms, so that still looks wrong to me.
– Ben
Oct 25 '18 at 3:01