Possible reasons for difference between pg.partial_corr and function from stack?

[leodtprojectsd]

Hi guys,

I am tyring to use pingouin.partial_corr to get the SPEARMAN partial correlation of all variables in a matrix.
I tried to compare the output of pingouin.partial_corr to one from a function from stack, however there are different resutls.

I was hoping you would know what the differences are attributed to (if its a bug on your side or not):

Here is the code:

#!pip install -q pingouin #uncomment if not installed 
import pingouin as pg

import pandas as pd
import numpy as np

from itertools import combinations
from sklearn import linear_model
from sklearn.datasets import load_iris

import warnings
warnings.simplefilter(action='ignore', category=FutureWarning) #too many warnings


#STACK PARTIAL CORRELATION https://stackoverflow.com/questions/73633787/get-partial-correlations-matrix-from-pandas-dataframe-using-spearman/73830455?noredirect=1#comment130427144_73830455
def part_corr_stack(xyz):
    var1, var2, rest = *xyz[0], xyz[1]
    var1_reg = linear_model.LinearRegression().fit(iris_data[rest], iris_data[var1])
    var2_reg = linear_model.LinearRegression().fit(iris_data[rest], iris_data[var2])
    var1_res = iris_data[var1] - var1_reg.predict(iris_data[rest])
    var2_res = iris_data[var2] - var2_reg.predict(iris_data[rest])
    part_corr_df = pd.concat([var1_res, var2_res], axis=1).corr(method='spearman')
    return part_corr_df.unstack()
#WITH PINGOUIN
def part_corr_ping(xyz,  df_='df'):
  var1, var2, rest = *xyz[0], xyz[1]
  part_corr_pingouin  = pg.partial_corr(data=df_, x=var1, y=var2, covar=rest, method='spearman').reset_index()
  xyrest_corr = pd.concat([pd.DataFrame({'x':var1, 'y':var2,'rest': [rest]}),part_corr_pingouin], axis=1 )
  return xyrest_corr

# MAIN 
 # import data
iris_data = load_iris()
iris_data = pd.DataFrame(iris_data['data'], columns=iris_data['feature_names'])

#get all the pairs of variables
xy_combinations = list(combinations(iris_data.columns, 2))
z = [[col for col in iris_data.columns if col not in xy] for xy in xy_combinations]
xyz_combinations = list(zip(xy_combinations, z))


corr_part_stack = pd.DataFrame(pd.concat(list(map(part_corr_stack, xyz_combinations))), columns=['part_corr']).reset_index().rename(columns={'level_0' :'x', 'level_1':'y', 'part_corr':'r_part_stack'})
 #modified version fo stack 
 # partial with pingou
with warnings.catch_warnings(): # ignore: RuntimeWarning:divide by zero encountered in true_divide in this function. 
    warnings.simplefilter("ignore")
    corr_part_pingouin_melt = pd.concat([part_corr_ping(p, df_=iris_data) for p in xyz_combinations]).rename(columns={'r' :'r_part_ping'})

corr_df = (
           corr_part_stack.set_index(['x','y'])
           .merge(corr_part_pingouin_melt.set_index(['x','y'])['r_part_ping'], left_index=True, right_index=True, how='right') #merge on index, keep only index of corr_part_pingouin_melt, that is the smallest. 
           .reset_index()
  )

OUT

x	y	r_part_stack	r_part_ping
petal length (cm)	petal width (cm)	0.862649	0.755538
sepal length (cm)	petal length (cm)	0.681566	0.545328
sepal length (cm)	petal width (cm)	-0.303597	0.045292
sepal length (cm)	sepal width (cm)	0.615629	0.237173
sepal width (cm)	petal length (cm)	-0.633985	-0.223860
sepal width (cm)	petal width (cm)	0.362407	-0.006916

[raphaelvallat]

Hi @leodtprojectsd,

Please see #147 (comment) and https://pingouin-stats.org/changelog.html#v0-4-0-august-2021

Briefly, since v0.4.0, Pingouin uses an inverse-covariance matrix approach to calculate the partial correlation, instead of a regression-residuals approach (which you implemented in part_corr_stack, and was the default prior to 0.4.0). Tho reasons for this:

1. This is the same behavior as R (ppcor package)
2. This method is more numerically stable than the residuals method for Spearman correlations, especially when there are ties in the data and/or one of the variable is binary.

In the iris_dataset, there are many duplicate values, which will lead to ties when calculating ranks. I would therefore recommend using a regular Pearson partial correlation:

In [8]: iris_data.shape
Out[8]: (150, 4)

In [6]: iris_data.nunique()
Out[6]: 
sepal length (cm)    35
sepal width (cm)     23
petal length (cm)    43
petal width (cm)     22
dtype: int64

[leodtprojectsd]

@raphaelvallat thanks!
just to be sure I understand, you do not reccomend using Pingouin.partial_corr in the case of high number of variables with few unique values, correct?

Do you have another option in Pingoung for high # of variables, low # of unique values?

[raphaelvallat]

I do not recommend using a Spearman partial correlation (pg.partial_corr(method="spearman")) when you have one or more variables with ties (= duplicates). The default Pearson partial correlation (pg.partial_corr()) should be fine though.