从（缩小的）协方差矩阵计算（部分）相关性（帮助将 R 代码移植到 Python）

Question

有一篇论文我觉得很有趣，想使用Python中的一些方法。 Erb 等人。 2020 对成分数据实现部分相关，Jin 等人。 2022 在名为 Propr 的 R 包中实现它。

我找到了我在下面简化的函数

bShrink

：

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

现在我正在尝试将其移植到 Python 中。

Cb

之间存在一些细微差异，这是预期的，但

PC

中存在主要差异，即偏相关（cor2pcor 函数）。

我尝试使用这里的答案，但无法让它工作：Python 中的部分相关性

这是我的Python代码：

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

我试图避免使用 Pinouin 或除 numpy、pandas 和 sklearn 之外的任何其他包。

如何从缩小的协方差矩阵创建偏相关矩阵？

Answer 1

您似乎正在尝试将

bShrink

函数从 R 移植到 Python，并且在部分相关计算中遇到了问题。 R 代码使用

cor2pcor

函数来计算偏相关矩阵，该矩阵在 Python 标准库中可能没有直接等效项。但是，您可以使用 numpy 和 scipy 实现偏相关计算。

以下是如何修改 Python 代码来计算偏相关矩阵：

import numpy as np
import pandas as pd
from sklearn.covariance import LedoitWolf
from scipy.linalg import sqrtm

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 / np.sqrt(Z)

    return pd.DataFrame(partial, index=components, columns=components)

X = pd.read_csv("https://pastebin.com/raw/e3BSEZiK", sep="\t", index_col=0)

partial_corr_matrix = bShrink(X)
print(partial_corr_matrix)

在此修改后的代码中，在计算部分相关矩阵时，我将除法替换为

np.sqrt(Z)

而不是仅使用

Z

。这应该可以帮助您获得更接近在 R 中获得的结果。

请注意，虽然这种方法可以提供与 R 代码类似的结果，但由于底层数值计算，直接移植仍可能导致轻微的差异。如有必要，您可能需要微调代码或咨询数值线性代数专家进行进一步调整。