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.
 A: 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:


*

*fitting a simple linear (or logit) model to estimate $E[y_{1i}|x_{i},s_i]$ and $E[y_{2i}|x_{i},s_i]$.

*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]})$$

*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.
