Partial correlation

Computation

Using linear regression

A simple way to compute the sample partial correlation for some data is to solve the two associated linear regression problems and calculate the correlation between the residuals. Let X and Y be random variables taking real values, and let Z be the n-dimensional vector-valued random variable. Let x_i, y_i and z_i denote the ith of ${\displaystyle N}$ i.i.d. observations from some joint probability distribution over real random variables X, Y, and Z, with z_i having been augmented with a 1 to allow for a constant term in the regression. Solving the linear regression problem amounts to finding (n+1)-dimensional regression coefficient vectors ${\displaystyle {\mathbf{w}}_X^{\ast }}$ and ${\displaystyle {\mathbf{w}}_Y^{\ast }}$ such that

${\displaystyle {\mathbf{w}}_X^{\ast }=\arg \underset{\mathbf{w}}{min}\left\{\sum _{i=1}^N(x_i-⟨\mathbf{w},{\mathbf{z}}_i⟩{)}^2\right\}}$

${\displaystyle {\mathbf{w}}_Y^{\ast }=\arg \underset{\mathbf{w}}{min}\left\{\sum _{i=1}^N(y_i-⟨\mathbf{w},{\mathbf{z}}_i⟩{)}^2\right\}}$

where ${\displaystyle N}$ is the number of observations, and ${\displaystyle ⟨\mathbf{w},{\mathbf{z}}_i⟩}$ is the scalar product between the vectors ${\displaystyle \mathbf{w}}$ and ${\displaystyle {\mathbf{z}}_i}$.

The residuals are then

${\displaystyle e_{X,i}=x_i-⟨{\mathbf{w}}_X^{\ast },{\mathbf{z}}_i⟩}$

${\displaystyle e_{Y,i}=y_i-⟨{\mathbf{w}}_Y^{\ast },{\mathbf{z}}_i⟩}$

and the sample partial correlation is then given by the usual formula for sample correlation, but between these new derived values:

${\displaystyle \begin{array}{cc}{\hat{\rho }}_{XY\cdot \mathbf{Z}}& =\frac{N\sum _{i=1}^N{e}_{X,i}{e}_{Y,i}-\sum _{i=1}^N{e}_{X,i}\sum _{i=1}^N{e}_{Y,i}}{\sqrt{N\sum _{i=1}^N{e}_{X,i}^2-{\left(\sum _{i=1}^N{e}_{X,i}\right)}^2}\sqrt{N\sum _{i=1}^N{e}_{Y,i}^2-{\left(\sum _{i=1}^N{e}_{Y,i}\right)}^2}}\\ & =\frac{N\sum _{i=1}^N{e}_{X,i}{e}_{Y,i}}{\sqrt{N\sum _{i=1}^N{e}_{X,i}^2}\sqrt{N\sum _{i=1}^N{e}_{Y,i}^2}}.\end{array}}$

In the first expression the three terms after minus signs all equal 0 since each contains the sum of residuals from an ordinary least squares regression.

Example

Consider the following data on three variables, X, Y, and Z:

Computing the Pearson correlation coefficient between variables X and Y results in approximately 0.970, while computing the partial correlation between X and Y, using the formula given above, gives a partial correlation of 0.919. The computations were done using R with the following code.

<span id="mwpg"></span><span class="gp" id="mwpw">> </span><span class="n" id="mwqA">x</span><span class="w" id="mwqQ"> </span><span class="o" id="mwqg"><-</span><span class="w" id="mwqw"> </span><span class="nf" id="mwrA">c</span><span class="p" id="mwrQ">(</span><span class="m" id="mwrg">2</span><span class="p" id="mwrw">,</span><span class="w" id="mwsA"> </span><span class="m" id="mwsQ">4</span><span class="p" id="mwsg">,</span><span class="w" id="mwsw"> </span><span class="m" id="mwtA">15</span><span class="p" id="mwtQ">,</span><span class="w" id="mwtg"> </span><span class="m" id="mwtw">20</span><span class="p" id="mwuA">)</span> <span class="gp" id="mwuQ">> </span><span class="n" id="mwug">y</span><span class="w" id="mwuw"> </span><span class="o" id="mwvA"><-</span><span class="w" id="mwvQ"> </span><span class="nf" id="mwvg">c</span><span class="p" id="mwvw">(</span><span class="m" id="mwwA">1</span><span class="p" id="mwwQ">,</span><span class="w" id="mwwg"> </span><span class="m" id="mwww">2</span><span class="p" id="mwxA">,</span><span class="w" id="mwxQ"> </span><span class="m" id="mwxg">3</span><span class="p" id="mwxw">,</span><span class="w" id="mwyA"> </span><span class="m" id="mwyQ">4</span><span class="p" id="mwyg">)</span> <span class="gp" id="mwyw">> </span><span class="n" id="mwzA">z</span><span class="w" id="mwzQ"> </span><span class="o" id="mwzg"><-</span><span class="w" id="mwzw"> </span><span class="nf" id="mw0A">c</span><span class="p" id="mw0Q">(</span><span class="m" id="mw0g">0</span><span class="p" id="mw0w">,</span><span class="w" id="mw1A"> </span><span class="m" id="mw1Q">0</span><span class="p" id="mw1g">,</span><span class="w" id="mw1w"> </span><span class="m" id="mw2A">1</span><span class="p" id="mw2Q">,</span><span class="w" id="mw2g"> </span><span class="m" id="mw2w">1</span><span class="p" id="mw3A">)</span> <span class="go" id="mw3Q"># regress x onto z and compute residuals</span> <span class="gp" id="mw3g">> </span><span class="n" id="mw3w">res_x</span><span class="w" id="mw4A"> </span><span class="o" id="mw4Q"><-</span><span class="w" id="mw4g"> </span><span class="nf" id="mw4w">lm</span><span class="p" id="mw5A">(</span><span class="n" id="mw5Q">x</span><span class="w" id="mw5g"> </span><span class="o" id="mw5w">~</span><span class="w" id="mw6A"> </span><span class="n" id="mw6Q">z</span><span class="p" id="mw6g">)</span><span class="o" id="mw6w">$</span><span class="n" id="mw7A">residuals</span> <span class="go" id="mw7Q"># regress y onto z and compute residuals</span> <span class="gp" id="mw7g">> </span><span class="n" id="mw7w">res_y</span><span class="w" id="mw8A"> </span><span class="o" id="mw8Q"><-</span><span class="w" id="mw8g"> </span><span class="nf" id="mw8w">lm</span><span class="p" id="mw9A">(</span><span class="n" id="mw9Q">y</span><span class="w" id="mw9g"> </span><span class="o" id="mw9w">~</span><span class="w" id="mw-A"> </span><span class="n" id="mw-Q">z</span><span class="p" id="mw-g">)</span><span class="o" id="mw-w">$</span><span class="n" id="mw_A">residuals</span> <span class="go" id="mw_Q"># compute correlation of residuals</span> <span class="gp" id="mw_g">> </span><span class="nf" id="mw_w">cor</span><span class="p" id="mwAQA">(</span><span class="n" id="mwAQE">res_x</span><span class="p" id="mwAQI">,</span><span class="w" id="mwAQM"> </span><span class="n" id="mwAQQ">res_y</span><span class="p" id="mwAQU">)</span> <span class="go" id="mwAQY"># [1] 0.919145</span> <span class="go" id="mwAQc"># show this is distinct from the correlation between x and y</span> <span class="gp" id="mwAQg">> </span><span class="nf" id="mwAQk">cor</span><span class="p" id="mwAQo">(</span><span class="n" id="mwAQs">x</span><span class="p" id="mwAQw">,</span><span class="w" id="mwAQ0"> </span><span class="n" id="mwAQ4">y</span><span class="p" id="mwAQ8">)</span> <span class="go" id="mwARA"># [1] 0.9695016</span> <span class="go" id="mwARE"># compute generalized partial correlations</span> <span class="gp" id="mwARI">> </span><span class="n" id="mwARM">generalCorr</span><span class="o" id="mwARQ">::</span><span class="nf" id="mwARU">parcorMany</span><span class="p" id="mwARY">(</span><span class="w" id="mwARc"> </span><span class="nf" id="mwARg">cbind</span><span class="p" id="mwARk">(</span><span class="n" id="mwARo">x</span><span class="p" id="mwARs">,</span><span class="w" id="mwARw"> </span><span class="n" id="mwAR0">y</span><span class="p" id="mwAR4">,</span><span class="w" id="mwAR8"> </span><span class="n" id="mwASA">z</span><span class="p" id="mwASE">)</span><span class="w" id="mwASI"> </span><span class="p" id="mwASM">)</span> <span class="go" id="mwASQ"># nami namj partij partji rijMrji </span> <span class="go" id="mwASU"># [1,] "x" "y" "0.8844" "1" "-0.1156"</span> <span class="go" id="mwASY"># [2,] "x" "z" "0.1581" "1" "-0.8419"</span>

The lower section of the code reports the generalized nonlinear partial correlation coefficient between X and Y after removing the nonlinear effect of Z (0.8844). It also reports the generalized nonlinear partial correlation coefficient between X and Z after removing the nonlinear effect of Y (0.1581). Additional documentation is available in the R package generalCorr, including package vignettes.

Using recursive formula

It can be computationally expensive to solve the linear regression problems. Actually, the nth-order partial correlation (i.e., with |Z| = n) can be easily computed from three (n - 1)th-order partial correlations. The zeroth-order partial correlation ρ_XY·Ø is defined to be the regular correlation coefficient ρ_XY.

It holds, for any ${\displaystyle Z_0\in \mathbf{Z},}$ that

${\displaystyle {\rho }_{XY\cdot \mathbf{Z}}=\frac{{\rho }_{XY\cdot \mathbf{Z}\setminus \{Z_0\}}-{\rho }_{X{Z}_0\cdot \mathbf{Z}\setminus \{Z_0\}}{\rho }_{Z_{0}Y\cdot \mathbf{Z}\setminus \{Z_0\}}}{\sqrt{1-{\rho }_{X{Z}_0\cdot \mathbf{Z}\setminus \{Z_0\}}^2}\sqrt{1-{\rho }_{Z_{0}Y\cdot \mathbf{Z}\setminus \{Z_0\}}^2}}}$

Naïvely implementing this computation as a recursive algorithm yields an exponential time complexity. However, this computation has the overlapping subproblems property, such that using dynamic programming or simply caching the results of the recursive calls yields a complexity of ${\displaystyle \mathcal{O}(n^3)}$.

Note in the case where Z is a single variable, this reduces to:

${\displaystyle {\rho }_{XY\cdot Z}=\frac{{\rho }_{XY}-{\rho }_{XZ}{\rho }_{ZY}}{\sqrt{1-{\rho }_{XZ}^2}\sqrt{1-{\rho }_{ZY}^2}}}$

Using matrix inversion

The partial correlation can also be written in terms of the joint precision matrix. Consider a set of random variables, ${\displaystyle \mathbf{V}=X_1,\dots ,X_n}$ of cardinality n. We want the partial correlation between two variables ${\displaystyle X_i}$ and ${\displaystyle X_j}$ given all others, i.e., ${\displaystyle \mathbf{V}\setminus \{X_i,X_j\}}$. Suppose the (joint/full) covariance matrix ${\displaystyle \mathrm{\Sigma }=({\sigma }_{ij})}$ is positive definite and therefore invertible. If the precision matrix is defined as ${\displaystyle \mathrm{\Omega }=(p_{ij})={\mathrm{\Sigma }}^{-1}}$, then

Computing this requires ${\displaystyle {\mathrm{\Sigma }}^{-1}}$, the inverse of the covariance matrix ${\displaystyle \mathrm{\Sigma }}$ which runs in ${\displaystyle \mathcal{O}(n^3)}$ time (using the sample covariance matrix to obtain a sample partial correlation). Note that only a single matrix inversion is required to give all the partial correlations between pairs of variables in ${\displaystyle \mathbf{V}}$.

To prove Equation (1), return to the previous notation (i.e. ${\displaystyle X,Y,\mathbf{Z}\leftrightarrow X_i,X_j,\mathbf{V}\setminus \{X_i,X_j\}}$) and start with the definition of partial correlation: ρ_XY·Z is the correlation between the residuals e_X and e_Y resulting from the linear regression of X with Z and of Y with Z, respectively.

First, suppose ${\displaystyle \beta ,\gamma }$ are the coefficients for linear regression fit; that is,

${\displaystyle \beta ={\mathrm{argmin}}_{\beta }\mathbb{E}\Vert X-{\beta }^{T}Z{\Vert }^2}$

${\displaystyle \gamma ={\mathrm{argmin}}_{\gamma }\mathbb{E}\Vert Y-{\gamma }^{T}Z{\Vert }^2}$

Write the joint covariance matrix for the vector ${\displaystyle (X,Y,Z^T{)}^T}$ as

${\displaystyle \mathrm{\Sigma }=\left[\begin{array}{ccc}{\mathrm{\Sigma }}_{XX}& {\mathrm{\Sigma }}_{XY}& {\mathrm{\Sigma }}_{XZ}\\ {\mathrm{\Sigma }}_{YX}& {\mathrm{\Sigma }}_{YY}& {\mathrm{\Sigma }}_{YZ}\\ {\mathrm{\Sigma }}_{ZX}& {\mathrm{\Sigma }}_{ZY}& {\mathrm{\Sigma }}_{ZZ}\end{array}\right]=\left[\begin{array}{cc}{C}_{11}& C_{12}\\ C_{21}& C_{22}\end{array}\right]}$

where$${\displaystyle C_{11}=\left[\begin{array}{cc}{\mathrm{\Sigma }}_{XX}& {\mathrm{\Sigma }}_{XY}\\ {\mathrm{\Sigma }}_{YX}& {\mathrm{\Sigma }}_{YY}\end{array}\right],C_{12}=\left[\begin{array}{c}{\mathrm{\Sigma }}_{XZ}\\ {\mathrm{\Sigma }}_{YZ}\end{array}\right],C_{21}=\left[\begin{array}{cc}{\mathrm{\Sigma }}_{ZX}& {\mathrm{\Sigma }}_{ZY}\end{array}\right],C_{22}={\mathrm{\Sigma }}_{ZZ}}$$Then the standard formula for linear regression gives

${\displaystyle \beta ={\left({\mathrm{\Sigma }}_{ZZ}\right)}^{-1}{\mathrm{\Sigma }}_{ZX}}$

Hence, the residuals can be written as

${\displaystyle R_X=X-{\beta }^{T}Z=X-{\mathrm{\Sigma }}_{XZ}{\left({\mathrm{\Sigma }}_{ZZ}\right)}^{-1}Z}$

Note that ${\displaystyle R_X}$ has expectation zero because of the inclusion of an intercept term in ${\displaystyle Z}$. Computing the covariance now gives

Next, write the precision matrix ${\displaystyle \mathrm{\Omega }={\mathrm{\Sigma }}^{-1}}$ in a similar block form:

${\displaystyle \mathrm{\Omega }=\left[\begin{array}{ccc}{\mathrm{\Omega }}_{XX}& {\mathrm{\Omega }}_{XY}& {\mathrm{\Omega }}_{XZ}\\ {\mathrm{\Omega }}_{YX}& {\mathrm{\Omega }}_{YY}& {\mathrm{\Omega }}_{YZ}\\ {\mathrm{\Omega }}_{ZX}& {\mathrm{\Omega }}_{ZY}& {\mathrm{\Omega }}_{ZZ}\end{array}\right]=\left[\begin{array}{cc}{P}_{11}& P_{12}\\ P_{21}& P_{22}\end{array}\right]}$

Then, by Schur's formula for block-matrix inversion,

${\displaystyle P_{11}^{-1}=C_{11}-C_{12}{C}_{22}^{-1}{C}_{21}}$

The entries of the right-hand-side matrix are precisely the covariances previously computed in (2), giving

${\displaystyle P_{11}^{-1}=\left[\begin{array}{cc}\mathrm{Cov}(R_X,R_X)& \mathrm{Cov}(R_X,R_Y)\\ \mathrm{Cov}(R_Y,R_X)& \mathrm{Cov}(R_Y,R_Y)\end{array}\right]}$

Using the formula for the inverse of a 2×2 matrix gives

${\displaystyle \begin{array}{cc}{P}_{11}^{-1}& =\frac{1}{\text{det}{P}_{11}}\left(\begin{array}{cc}[P_{11}{]}_{22}& -[P_{11}{]}_{12}\\ -[P_{11}{]}_{21}& [P_{11}{]}_{11}\end{array}\right)\\ & =\frac{1}{\text{det}{P}_{11}}\left(\begin{array}{cc}{p}_{YY}& -p_{XY}\\ -p_{YX}& p_{XX}\end{array}\right)\end{array}}$

So indeed, the partial correlation is

${\displaystyle {\rho }_{XY\cdot Z}=\frac{\mathrm{Cov}(R_X,R_Y)}{\sqrt{\mathrm{Cov}(R_X,R_X)\mathrm{Cov}(R_Y,R_Y)}}=\frac{-\frac{1}{\text{det}{P}_{11}}{p}_{XY}}{\sqrt{\frac{1}{\text{det}{P}_{11}}{p}_{XX}\frac{1}{\text{det}{P}_{11}}{p}_{YY}}}=-\frac{p_{XY}}{\sqrt{p_{XX}{p}_{YY}}}}$

as claimed in (1).

Interpretation

Geometrical

Let three variables X, Y, Z (where Z is the "control" or "extra variable") be chosen from a joint probability distribution over n variables V. Further, let v_i, 1 ≤ i ≤ N, be N n-dimensional i.i.d. observations taken from the joint probability distribution over V. The geometrical interpretation comes from considering the N-dimensional vectors x (formed by the successive values of X over the observations), y (formed by the values of Y), and z (formed by the values of Z).

It can be shown that the residuals e_X,i coming from the linear regression of X on Z, if also considered as an N-dimensional vector e_X (denoted r_X in the accompanying graph), have a zero scalar product with the vector z generated by Z. This means that the residuals vector lies on an (N–1)-dimensional hyperplane S_z that is perpendicular to z.

The same also applies to the residuals e_Y,i generating a vector e_Y. The desired partial correlation is then the cosine of the angle φ between the projections e_X and e_Y of x and y, respectively, onto the hyperplane perpendicular to z.

As conditional independence test

With the assumption that all involved variables are multivariate Gaussian, the partial correlation ρ_XY·Z is zero if and only if X is conditionally independent from Y given Z. This property does not hold in the general case.

To test if a sample partial correlation ${\displaystyle {\hat{\rho }}_{XY\cdot \mathbf{Z}}}$ implies that the true population partial correlation differs from 0, Fisher's z-transform of the partial correlation can be used:

${\displaystyle z({\hat{\rho }}_{XY\cdot \mathbf{Z}})=\frac{1}{2}\ln \left(\frac{1+{\hat{\rho }}_{XY\cdot \mathbf{Z}}}{1-{\hat{\rho }}_{XY\cdot \mathbf{Z}}}\right)}$

The null hypothesis is ${\displaystyle H_0:{\rho }_{XY\cdot \mathbf{Z}}=0}$, to be tested against the two-tail alternative ${\displaystyle H_A:{\rho }_{XY\cdot \mathbf{Z}}\ne 0}$. ${\displaystyle H_0}$ can be rejected if

${\displaystyle \sqrt{N-|\mathbf{Z}|-3}\cdot |z({\hat{\rho }}_{XY\cdot \mathbf{Z}})|>{\mathrm{\Phi }}^{-1}(1-\alpha /2)}$

where ${\displaystyle \mathrm{\Phi }}$ is the cumulative distribution function of a Gaussian distribution with zero mean and unit standard deviation, ${\displaystyle \alpha }$ is the significance level of ${\displaystyle H_0}$, and ${\displaystyle N}$ is the sample size. This z-transform is approximate, and the actual distribution of the sample (partial) correlation coefficient is not straightforward. However, an exact t-test based on a combination of the partial regression coefficient, the partial correlation coefficient, and the partial variances is available.

The distribution of the sample partial correlation was described by Fisher.

Semipartial correlation (part correlation)

The semipartial (or part) correlation statistic is similar to the partial correlation statistic; both compare variations of two variables after certain factors are controlled for. However, to calculate the semipartial correlation, one holds the third variable constant for either X or Y but not both; whereas for the partial correlation, one holds the third variable constant for both. The semipartial correlation compares the unique variation of one variable (having removed variation associated with the Z variable(s)) with the unfiltered variation of the other, while the partial correlation compares the unique variation of one variable to the unique variation of the other.

The semipartial correlation can be viewed as more practically relevant "because it is scaled to (i.e., relative to) the total variability in the dependent (response) variable." Conversely, it is less theoretically useful because it is less precise about the role of the unique contribution of the independent variable.

The absolute value of the semipartial correlation of X with Y is always less than or equal to that of the partial correlation of X with Y. The reason is this: Suppose the correlation of X with Z has been removed from X, giving the residual vector e_x . In computing the semipartial correlation, Y still contains both unique variance and variance due to its association with Z. But e_x , being uncorrelated with Z, can only explain some of the unique part of the variance of Y and not the part related to Z. In contrast, with the partial correlation, only e_y (the part of the variance of Y that is unrelated to Z) is to be explained, so there is less variance of the type that e_x cannot explain.

Use in time series analysis

In time series analysis, the partial autocorrelation function (sometimes "partial correlation function") of a time series is defined, for lag ${\displaystyle h}$, as

${\displaystyle \phi (h)={\rho }_{X_0{X}_h\cdot \{X_1,\dots ,X_{h-1}\}}}$

This function is used to determine the appropriate lag length for an autoregression.

Partial correlations with shrinkage

When the sample size is smaller than the number of variables, a.k.a. high-dimensional setting, estimating partial correlations can be challenging. In this scenario, the sample covariance ${\displaystyle \hat{\mathrm{\Sigma }}}$ is not well-conditioned, and finding its inverse ${\displaystyle \hat{\mathrm{\Omega }}}$ turns problematic.

Shrinkage estimation methods improve ${\displaystyle \hat{\mathrm{\Sigma }}}$ or ${\displaystyle \hat{\mathrm{\Omega }}}$ and produces more reliable partial correlation estimates. One example is the Ledoit-Wolf shrinkage estimator,

${\displaystyle {\hat{\mathrm{\Sigma }}}^{[\lambda ]}=\lambda T+(1-\lambda )\mathrm{\Sigma }}$

where ${\displaystyle \hat{\mathrm{\Sigma }}}$ is the sample covariance matrix, ${\displaystyle T}$ is a target matrix (e.g., a diagonal matrix), and the shrinkage intensity ${\displaystyle \lambda \in (0,1)}$.

The partial correlation under the Ledoit-Wolf shrinkage is then:

${\displaystyle {\hat{P}}_{ij}^{[\lambda ]}=\frac{{\hat{\mathrm{\Omega }}}_{ij}^{[\lambda ]}}{\sqrt{{\hat{\mathrm{\Omega }}}_{ii}^{[\lambda ]}{\hat{\mathrm{\Omega }}}_{jj}^{[\lambda ]}}}}$

where ${\displaystyle {\hat{\mathrm{\Omega }}}_{ij}^{[\lambda ]}}$ is the inverse of ${\displaystyle {\hat{\mathrm{\Sigma }}}_{ij}^{[\lambda ]}}$. This method is used in a variety of fields including finance and genomics.