Question: There's a paper that I found interesting and would like to use some of the methods in Python. Erb et al. 2020 implements partial correlation on compositional data and Jin et al. 2022 implements it in an R package called Propr. I found the function bShrink that I'm simplifying below: library(corpcor) # Load iris dataset (not compositional but will work for this case) X = read.table("https://pastebin.com/raw/e3BSEZiK", sep = "\t", row.names = 1, header = TRUE, check.names = FALSE) bShrink <- function(M){ # transform counts to log proportions P <- M / rowSums(M) B <- log(P) # covariance shrinkage D <- ncol(M) Cb <- cov.shrink(B,verbose=FALSE) G <- diag(rep(1,D))-matrix(1/D,D,D) Cov <- G%*%Cb%*%G # partial correlation PC <- cor2pcor(Cov) return(PC) } > bShrink(X) [,1] [,2] [,3] [,4] [1,] 1.0000000 0.96409509 0.6647093 -0.23827651 [2,] 0.9640951 1.00000000 -0.4585507 0.02735205 [3,] 0.6647093 -0.45855072 1.0000000 0.85903005 [4,] -0.2382765 0.02735205 0.8590301 1.00000000  Now I'm trying to port this in Python. Getting some small differences between Cb which is expected but major differences in PC which is the partial correlation (cor2pcor function). I tried using the answers from here but couldn't get it to work: Partial Correlation in Python Here's my Python code: import numpy as np import pandas as pd from sklearn.covariance import LedoitWolf def bShrink(M:pd.DataFrame): components = M.columns M = M.values P = M/M.sum(axis=1).reshape(-1,1) B = np.log(P) D = M.shape[1] lw_model = LedoitWolf() lw_model.fit(B) Cb = lw_model.covariance_ G = np.eye(D) - np.ones((D, D)) / D Cov = G @ Cb @ G precision = np.linalg.inv(Cov) diag = np.diag(precision) Z = np.outer(diag, diag) partial = -precision / Z return pd.DataFrame(partial, index=components, columns=components) X = pd.read_csv("https://pastebin.com/raw/e3BSEZiK", sep="\t", index_col=0) bShrink(X) # sepal_length sepal_width petal_length petal_width # sepal_length -5.551115e-17 -5.551115e-17 -5.551115e-17 -5.551115e-17 # sepal_width -5.551115e-17 -5.551115e-17 -5.551115e-17 -5.551115e-17 # petal_length -5.551115e-17 -5.551115e-17 -5.551115e-17 -5.551115e-17 # petal_width -5.551115e-17 -5.551115e-17 -5.551115e-17 -5.551115e-17  I'm trying to avoid using Pingouin or any other packages than numpy, pandas, and sklearn. How can I create a partial correlation matrix from a shrunken covariance matrix?

Answer: The problems There are two issues here, the estimation of the covariate matrices are different and the translation of the cor2pcor step is wrong. Let's start with cor2pcor. I honestly did not understand very well how you are computing it so I decided to check the implementation in R. If we go on the corpcor implementation we get the function: cor2pcor = function(m, tol) { # invert, then negate off-diagonal entries m = -pseudoinverse(m, tol=tol) diag(m) = -diag(m) # standardize and return return(cov2cor(m)) }  We can also check the implementation of the R stats package to see the implementation of cov2cor: cov2cor <- function(V) { ## Purpose: Covariance matrix |--> Correlation matrix -- efficiently ## ---------------------------------------------------------------------- ## Arguments: V: a covariance matrix (i.e. symmetric and positive definite) ## ---------------------------------------------------------------------- ## Author: Martin Maechler, Date: 12 Jun 2003, 11L:50 p <- (d <- dim(V))[1L] if(!is.numeric(V) || length(d) != 2L || p != d[2L]) stop("'V' is not a square numeric matrix") Is <- sqrt(1/diag(V)) # diag( 1/sigma_i ) if(any(!is.finite(Is))) warning("diag(.) had 0 or NA entries; non-finite result is doubtful") r <- V # keep dimnames r[] <- Is * V * rep(Is, each = p) ## == D %*% V %*% D where D = diag(Is) r[cbind(1L:p,1L:p)] <- 1 # exact in diagonal r }  With those two functions in mind, we can begin to translate them into Python, make sure to use pinv instead of inv because your matrix can be singular: def cov2cor(V): """ Convert Covariance matrix to Correlation matrix efficiently. Arguments: V: a covariance matrix (i.e. symmetric and positive definite) """ p, d = V.shape Is = np.sqrt(1 / np.diag(V)) # diag( 1/sigma_i ) r = V.copy() # keep dimnames r *= Is.reshape(-1, 1) * Is.reshape(1, -1) np.fill_diagonal(r, 1) # exact in diagonal return r def cor2pcor(m, tol=1e-15): """ Convert a correlation matrix to a partial correlation matrix efficiently. Arguments: m: a correlation matrix tol: tolerance for calculating the pseudo-inverse """ # Invert, then negate off-diagonal entries m = -np.linalg.pinv(m, rcond=tol) np.fill_diagonal(m, -np.diag(m)) # Standardize and return return cov2cor(m)  Finally, we can use those to implement b_shrink: def b_shrink(M: pd.DataFrame): components = M.columns M = M.values # transform counts to log proportions P = M / M.sum(axis=1).reshape(-1, 1) B = np.log(P) # Covariance shrinkage D = M.shape[1] Cb, shrinkage = ledoit_wolf(B) G = np.eye(D) - np.ones((D, D)) / D cov = G @ Cb @ G # Partial autocorrelation partial = cor2pcor(cov) return pd.DataFrame(partial, index=components, columns=components)  Doing this I got the following matrix, I replaced LedoitWolf with ledoit_wolf which is pretty much the same thing on this function:      sepal_length sepal_width petal_length petal_width     sepal_length 1 0.951383 0.606789 -0.188676   sepal_width 0.951383 1 -0.352785 -0.0561523   petal_length 0.606789 -0.352785 1 0.867358   petal_width -0.188676 -0.0561523 0.867358 1    The covariance matrix According to the documentation of cov.shrink the implementation in R is based on (Opgen-Rhein and Strimmer, 2007) and (Schäfer and Strimmer, 2005) which build on top of Ledoit-Wolf's work, while Sklearn's shrinkage implements only Ledoit-Wolf shrinkage. This gives slightly different matrices Cb in both languages: In R:  sepal_length sepal_width petal_length petal_width sepal_length 0.01384628 0.03049407 -0.03983544 -0.08249265 sepal_width 0.03049407 0.09092548 -0.10892279 -0.21334250 petal_length -0.03983544 -0.10892279 0.13876001 0.26559868 petal_width -0.08249265 -0.21334250 0.26559868 0.58653240  In Python: [[ 0.01425639 0.0286866 -0.03757796 -0.07812832] [ 0.0286866 0.09108648 -0.10768074 -0.21175137] [-0.03757796 -0.10768074 0.13876641 0.26434719] [-0.07812832 -0.21175137 0.26434719 0.58509169]]  We would get the same results if we replaced the Cb matrix from R in Python. As far as I know, there are no implementations of Rhein, Schäfer, and Strimmer's work in Python.