Can we isolate the effect of two different control variables in a semi-partial correlation?

Question

I know that we are able to use a partial correlation when we want to correlate X and Y but Z affects both of them and that we may use a semi-partial correlation when we want to correlate X and Y and we know that Z only affects X.

• In my case, I have two different "Zs"/third variables.The idea is the following: I want to correlate X1 and X2, but I know that Z1 is extremely correlated to X1 and that Z2 is extremely correlated to X2.

So, I'm not assuming that Z1 and Z2 are correlated to both X1 and X2 (as it would be via a partial correlation, right?) I'm not assuming that Z1 and Z2 only effects X1 or X2 (as it would be for a semi-partial)

• To put it simpler:

  Z1 is correlated to X1
  Z2 is correlated to X2
  I want to see if X1 and X2 are correlated controlling for the effects of Z1 on X1 and of Z2 on X2

Obs: all variables are continuos

• Edit1: reproducible data: HERE

Any ideas would be MUCH appreaciated! Thanks in advance!

• Edit2:

I was thinking about this solution after watching this video and from the kind comment below, I'm really hoping that this may be a good path :

If I apply the same logic as the video's for the partial corr, I guess that maybe it would be all right?

res1 <- rstandard(lm(X1 ~ Z1))
res 2 <- rstandard(lm(X2 ~ Z2))

so,

res1 = > the portion of X1 that is not explained by Z1
res2 = > the portion of X2 that is not explained by Z2 

so,

cor(res1,res2) = > should be the correlation between X1 (accouting for Z1's effect on it) and X2 (accounting for Z2's effect on it), right? 

• Question 1: I'm assuming this would be fine for a Pearson's correlation, but my data is non-parametric. Would it be okay to apply a Spearman's correlation cor(res1, res2) ?

• Question 2 : why rstandard ? Should we always apply this, even if data is non-parametric? (I guess so, since even Spearman corr is going to use Pearson's equation, after ranking the data, right?) but, still, I didn't quite understand why standartizing is a "must"
  
  

• Edit 3 : Using rstandard or not
  

So,

• When I standartize the residuals, I get a value of Pearson's rho cor(res1, res2) (0.95) and slope lm(res1 ~ res2) (0.95), which differs from Pearson's rho pcor.test(X1, X2, Z) (0.96) (I know the difference in this case seems minimal, but I'd want to understand why is there a difference in the first place...)

• when I don't standartize, the Pearson's rho(res1, res2) (0.96) is equal to Pearson's rho pcor.test(X1, X2, Z) (0.96) , but not to the slope lm(res1 ~ res2) (1.477e-04) anymore...why is that?

Answer 1

The problem

You say

I'm looking, ultimately [for] the relationship between $X_1$ and $X_2$

But since $Z_1$ and $Z_2$ are only related to one of both, you want to address this separately. Let's clarify first what concepts are we dealing with.

Definitions

Let us first recall the definition of correlation $\rho_{X_1,X_2}$:

$$ \rho_{X_1, X_2} = \frac{\mathbb{E}\bigg[(X_1 - \mathbb{E}[X_1])(X_2 - \mathbb{E}[X_2])\bigg]}{\sigma_{X_1}\sigma_{X_2}} $$

, where $\sigma_\cdot$ is the standard deviation of a specific variable.

Then, the semi-partial correlation you're interested in by removing the effects of $Z_1$ on $X_1$ and of $Z_2$ on $X_2$ is, that is $\rho_{X_1 | Z_1, X_2 | Z_2}$ is defined as:

$$ \rho_{X_1| Z_1, X_2 |Z_2} = \frac{\mathbb{E}\bigg[(\epsilon_{X_1 \sim Z_1})(\epsilon_{X_2 \sim Z_2})\bigg]}{\sigma_{\epsilon_{X_1 \sim Z_1}}\sigma_{\epsilon_{X_2 \sim Z_2}}} $$

, where $\epsilon_\cdot$ correspond to the errors from a model that regresses the variables in the subscript (e.g. the regression of $X_1$ by $Z_1$ if $\epsilon_{X_1 \sim Z_1}$). The differences from the mean disappear since errors are assumed to have mean zero: $\mathbb{E}[\epsilon_\cdot] = 0$

Lastly, let us recall that the slope of a simple linear regression of $Y$ on $X$ is:

$$ \beta_{Y \sim X} = \frac{\operatorname{Cov}[X, Y]}{\operatorname{Var}[X]} = \frac{\mathbb{E}\bigg[(X - \mathbb{E}[X])(Y - \mathbb{E}[Y]\bigg]}{\sigma_X \sigma_X} $$

We can see that the slope coefficient for the regression of the errors of $X_1$ on the errors of $X_2$, that is $\beta_{\epsilon_{X_1 \sim Z_1} \sim \epsilon_{X_2 \sim Z_2}}$, will be equal to our semi-partial correlation of interest only if $\sigma_X\sigma_X = \sigma_{\epsilon_{X_1 \sim Z_1}}\sigma_{\epsilon_{X_2 \sim Z_2}}$. With this, we can now answer your questions.

Answers

Answer 1

You say that your "data is non-parametric", which, although often described like this by applied researchers, does not make sense. A model can be parametric or non-parametric, data by itself has no parameters. Even more so, there is also the related common misconception that data must be normally distributed, rather than the errors from a model. What you mean by "okay" depends on what is your goal, and in this case whether the association as measured by ranks (which is what Spearman's correlation does) matters to you.

Answer 2

You don't necessarily need to standardize the residuals, as we see from the correlation definition that a correlation is a standardized covariance. Thus, computing a correlation will divide by the standard deviations anyways.

Follow-up answers

You then ask for why you get different answers from pcor.test(X1, X2, Z) than from the previous methods$^1$. Since you do not share your code, I will assume that $Z$ is a vector of both $Z_1$ and $Z_2$. Then, unless $Z_1$ is uncorrelated with $X_2$ and $Z_2$ is uncorrelated with $X_1$ (in which case doing a partial correlation on $Z_1$ and $Z_2$ would have been simpler), it is natural to expect the result to be different: $\rho_{X_1 | Z_1, X_2 | Z_2} \neq \rho_{X_1 | \{Z_1, Z_2\}, X_2 | \{Z_1, Z_2 \}}$.

Then, you ask why the correlation result is different if you do not standardize. Note our last equation in the definition section, the slope coefficient is the covariance divided by the standard deviation of your predictor squared (i.e. the variance). If you want the slope coefficient to equal the correlation, you need it to be divided by the product of standard deviations of both variables. This is exactly what happens when you standardize, since then the standard deviation (and variances) of the variables become one: $\sigma_{\epsilon^{\text{std}}_{X_1 | Z_1}} = \sigma_{\epsilon^{\text{std}}_{X_2 | Z_2}} = \sigma_{\epsilon^{\text{std}}_{X_1 | Z_1}} \times \sigma_{\epsilon^{\text{std}}_{X_1 | Z_1}} = \sigma_{\epsilon^{\text{std}}_{X_2 | Z_2}} \times \sigma_{\epsilon^{\text{std}}_{X_2 | Z_2}} = \sigma_{\epsilon^{\text{std}}_{X_1 | Z_1}} \times \sigma_{\epsilon^{\text{std}}_{X_2 | Z_2}} = 1$.

---

$^1$ This minimal difference could also just be approximation error from how different functions compute the same quantity.