协方差矩阵而非原始数据的典型相关分析

Canonical correlation analysis on covariance matrices instead of raw data

由于隐私问题,我没有原始原始数据矩阵,但我可以有 x 和 y (x'x, y'y, x'y) 数据集的协方差矩阵或它们之间的相关矩阵它们中的两个(或不是原始数据矩阵的任何其他类型的矩阵)。

我需要找到一种直接对这些矩阵应用典型相关分析的方法。浏览网络我没有找到任何解决我的问题的方法。我想问一下是否已经有一个实现的算法能够处理这些数据,R 最好,但其他语言也可以

cca 包 R 教程中的示例:
(https://stats.idre.ucla.edu/r/dae/canonical-correlation-analysis/)

mm <- read.csv("https://stats.idre.ucla.edu/stat/data/mmreg.csv")
colnames(mm) <- c("控制", "概念", "动机", "读", "写", "数学", “科学”、“性”)

您将数据集分为 x 和 y:
x <- 毫米 [ 1:3]
y <- mm[ 4:8]

然后该函数将这两个数据集作为输入:cc(x,y)(请注意,该函数会自行标准化数据)。

我想知道是否有一种方法可以通过以均值为中心的矩阵开始执行 cca: x = 比例尺(x,比例尺 = F) y = 比例尺(Y, 比例尺 = F)

然后计算协方差矩阵 x'x, y'y, xy'xy:
cvx = crossprod(x); cvy = crossprod(y); cvxy = crossprod(x,y)

算法应该输入这些矩阵来计算典型变量和相关系数 比如:f(cvx, cvy, cvxy)

这篇文章写了一个从协方差矩阵开始的解决方案,但我不知道它是否只是理论或有人实际实现了它 http://graphics.stanford.edu/courses/cs233-20-spring/ReferencedPapers/CCA_Weenik.pdf

希望够详尽!

简而言之:相关性在大多数(可能是所有)CCA 分析中都在内部使用。

总而言之:您需要根据具体情况了解如何做到这一点。下面我给大家举个例子。

什么是 Canonical-correlation 分析 (CCA)?

Canonical-correlation 分析 (CCA):帮助您确定可以在两个数据集之间创建的最佳线性关系。有关数据和使用库,请参阅 wikipedia. See references for examples. I will follow this post。 设置库、上传数据、select 一些变量、删除 nans、确定数据。

import pandas as pd
import numpy as np
df = pd.read_csv('2016 School Explorer.csv')
# choose relevant features
df = df[['Rigorous Instruction %',
      'Collaborative Teachers %',
     'Supportive Environment %',
       'Effective School Leadership %',
   'Strong Family-Community Ties %',
    'Trust %','Average ELA Proficiency',
       'Average Math Proficiency']]
df.corr()


# drop missing values
df = df.dropna()
# separate X and Y groups
X = df[['Rigorous Instruction %',
      'Collaborative Teachers %',
     'Supportive Environment %',
       'Effective School Leadership %',
   'Strong Family-Community Ties %',
    'Trust %'
      ]]
Y = df[['Average ELA Proficiency',
       'Average Math Proficiency']]


for col in X.columns:
    X[col] = X[col].str.strip('%')
    X[col] = X[col].astype('int')
# Standardise the data
from sklearn.preprocessing import StandardScaler
sc = StandardScaler(with_mean=True, with_std=True)
X_sc = sc.fit_transform(X)
Y_sc = sc.fit_transform(Y)

什么是相关性?

我在这里停下来谈谈想法和实施。 首先,CCA 分析自然是基于这个想法,但是对于数值分辨率,有不同的方法可以做到这一点。 来自 wikipedia 的定义。看图片: 我说这个是因为我要修改那个库的一个函数,我希望你真正注意它。 参见 Bilenko et al 2016 中的 Eq 4。但是你需要非常小心如何放置它。 请注意,严格来说您不需要相关性。

让我在 pyrrccahere

中展示计算该表达式的函数
def kcca(data, reg=0., numCC=None, kernelcca=True,
ktype='linear',
         gausigma=1.0, degree=2):
    """Set up and solve the kernel CCA eigenproblem
    """
    if kernelcca:
        kernel = [_make_kernel(d, ktype=ktype, gausigma=gausigma,
                               degree=degree) for d in data]
    else:
        kernel = [d.T for d in data]

    nDs = len(kernel)
    nFs = [k.shape[0] for k in kernel]
    numCC = min([k.shape[1] for k in kernel]) if numCC is None else numCC

    # Get the auto- and cross-covariance matrices
    crosscovs = [np.dot(ki, kj.T) for ki in kernel for kj in kernel]

    # Allocate left-hand side (LH) and right-hand side (RH):
    LH = np.zeros((sum(nFs), sum(nFs)))
    RH = np.zeros((sum(nFs), sum(nFs)))

    # Fill the left and right sides of the eigenvalue problem
    for i in range(nDs):
        RH[sum(nFs[:i]) : sum(nFs[:i+1]),
           sum(nFs[:i]) : sum(nFs[:i+1])] = (crosscovs[i * (nDs + 1)]
                                             + reg * np.eye(nFs[i]))

        for j in range(nDs):
            if i != j:
                LH[sum(nFs[:j]) : sum(nFs[:j+1]),
                   sum(nFs[:i]) : sum(nFs[:i+1])] = crosscovs[nDs * j + i]

    LH = (LH + LH.T) / 2.
    RH = (RH + RH.T) / 2.

    maxCC = LH.shape[0]
    r, Vs = eigh(LH, RH, eigvals=(maxCC - numCC, maxCC - 1))
    r[np.isnan(r)] = 0
    rindex = np.argsort(r)[::-1]
    comp = []
    Vs = Vs[:, rindex]
    for i in range(nDs):
        comp.append(Vs[sum(nFs[:i]):sum(nFs[:i + 1]), :numCC])
    return comp

此处的输出是典型协变量 (comp),它们是 Bilenko 等人 2016 年方程式 4 中的 a 和 b。

我只是想让你注意这个:

# Get the auto- and cross-covariance matrices
crosscovs = [np.dot(ki, kj.T) for ki in kernel for kj in kernel]

那正是那个操作发生的地方。请注意,这不完全是维基百科的定义,但在数学上是等价的。

相关性计算

我将按照维基百科中的方法计算相关性,但稍后我会修改该函数,因此它会涉及一些细节,以确保它清楚地回答了最初的问题。

# Get the auto- and cross-covariance matrices
crosscovs = [np.dot(ki, kj.T) for ki in kernel for kj in kernel]
print(crosscovs)
[array([[1217.        ,  746.04496925,  736.14178336,  575.21073838,
         517.52474332,  641.25363806],
       [ 746.04496925, 1217.        ,  732.6297358 , 1094.38480773,
         572.95747557, 1073.96490387],
       [ 736.14178336,  732.6297358 , 1217.        ,  559.5753228 ,
         682.15312862,  774.36607617],
       [ 575.21073838, 1094.38480773,  559.5753228 , 1217.        ,
         495.79248754, 1047.31981248],
       [ 517.52474332,  572.95747557,  682.15312862,  495.79248754,
        1217.        ,  632.75610906],
       [ 641.25363806, 1073.96490387,  774.36607617, 1047.31981248,
         632.75610906, 1217.        ]]), array([[367.74099904, 391.82683717],
       [348.78464015, 355.81358426],
       [440.88117453, 514.22183796],
       [326.32173163, 311.97282341],
       [216.32441793, 269.72859023],
       [288.27601974, 304.20209135]]), array([[367.74099904, 348.78464015, 440.88117453, 326.32173163,
        216.32441793, 288.27601974],
       [391.82683717, 355.81358426, 514.22183796, 311.97282341,
        269.72859023, 304.20209135]]), array([[1217.        , 1139.05867099],
       [1139.05867099, 1217.        ]])]

看看输出,我将在 -1 和 1 之间稍微更改一下。同样,这个修改很小。按照维基百科的定义,作者只关心分子,我现在只包括分母。

max_unit = 0
for crosscov in crosscovs:
    max_unit = np.max([max_unit,np.max(crosscov)])
"""I normalice"""
crosscovs_new = []
for crosscov in crosscovs:
    crosscovs_new.append(crosscov/max_unit)
print(crosscovs_new)
[array([[1.        , 0.6130197 , 0.60488232, 0.47264646, 0.4252463 ,
        0.52691342],
       [0.6130197 , 1.        , 0.6019965 , 0.89924799, 0.47079497,
        0.88246911],
       [0.60488232, 0.6019965 , 1.        , 0.45979895, 0.56052024,
        0.63629094],
       [0.47264646, 0.89924799, 0.45979895, 1.        , 0.40738906,
        0.86057503],
       [0.4252463 , 0.47079497, 0.56052024, 0.40738906, 1.        ,
        0.51993107],
       [0.52691342, 0.88246911, 0.63629094, 0.86057503, 0.51993107,
        1.        ]]), array([[0.30217009, 0.32196125],
       [0.28659379, 0.29236942],
       [0.36226884, 0.42253232],
       [0.26813618, 0.25634579],
       [0.17775219, 0.22163401],
       [0.2368743 , 0.24996063]]), array([[0.30217009, 0.28659379, 0.36226884, 0.26813618, 0.17775219,
        0.2368743 ],
       [0.32196125, 0.29236942, 0.42253232, 0.25634579, 0.22163401,
        0.24996063]]), array([[1.        , 0.93595618],
       [0.93595618, 1.        ]])]

为清楚起见,我将以稍微不同的方式向您展示原始数据的数字和相关性。

df.corr()
                          Average ELA Proficiency  Average Math Proficiency
Average ELA Proficiency                  1.000000                  0.935956
Average Math Proficiency                 0.935956                  1.000000

这也是一种查看变量名称的方法。我只是想告诉你上面的数字是有道理的,就是你所说的相关性。

CCA 的计算

所以现在我将稍微修改 pyrrcca 中的函数 kcca。这个想法是让该函数接受先前计算的相关矩阵。

from rcca import _make_kernel
from scipy.linalg import eigh

def kcca_working(data, reg=0.,
    numCC=None,
    kernelcca=False,
    ktype='linear',
    gausigma=1.0,
    degree=2,
    crosscovs=None):
    """Set up and solve the kernel CCA eigenproblem
    """
    if kernelcca:
        kernel = [_make_kernel(d, ktype=ktype, gausigma=gausigma,
                               degree=degree) for d in data]
    else:
        kernel = [d.T for d in data]

    nDs = len(kernel)
    nFs = [k.shape[0] for k in kernel]
    numCC = min([k.shape[1] for k in kernel]) if numCC is None else numCC

    if crosscovs is None:
        # Get the auto- and cross-covariance matrices
        crosscovs = [np.dot(ki, kj.T) for ki in kernel for kj in kernel]

    # Allocate left-hand side (LH) and right-hand side (RH):
    LH = np.zeros((sum(nFs), sum(nFs)))
    RH = np.zeros((sum(nFs), sum(nFs)))

    # Fill the left and right sides of the eigenvalue problem
    for i in range(nDs):
        RH[sum(nFs[:i]) : sum(nFs[:i+1]),
           sum(nFs[:i]) : sum(nFs[:i+1])] = (crosscovs[i * (nDs + 1)]
                                             + reg * np.eye(nFs[i]))

        for j in range(nDs):
            if i != j:
                LH[sum(nFs[:j]) : sum(nFs[:j+1]),
                   sum(nFs[:i]) : sum(nFs[:i+1])] = crosscovs[nDs * j + i]

    LH = (LH + LH.T) / 2.
    RH = (RH + RH.T) / 2.

    maxCC = LH.shape[0]
    r, Vs = eigh(LH, RH, eigvals=(maxCC - numCC, maxCC - 1))
    r[np.isnan(r)] = 0
    rindex = np.argsort(r)[::-1]
    comp = []
    Vs = Vs[:, rindex]
    for i in range(nDs):
        comp.append(Vs[sum(nFs[:i]):sum(nFs[:i + 1]), :numCC])
    return comp, crosscovs

让运行函数:

comp, crosscovs = kcca_working([X_sc, Y_sc], reg=0.,
         numCC=2, kernelcca=False, ktype='linear',
         gausigma=1.0, degree=2, crosscovs = crosscovs_new)
print(comp)
[array([[-0.00375779,  0.0078263 ],
       [ 0.00061439, -0.00357358],
       [-0.02054012, -0.0083491 ],
       [-0.01252477,  0.02976148],
       [ 0.00046503, -0.00905069],
       [ 0.01415084, -0.01264106]]), array([[ 0.00632283,  0.05721601],
       [-0.02606459, -0.05132531]])]

所以我采用了原函数,并且可以引入相关性,我也输出它只是为了检查。

我打印了典型协变量 (comp),它们是 Bilenko 等人 2016 年方程式 4 中的 a 和 b。

比较结果

现在我要比较原始函数和修改函数的结果。我会告诉你结果是等价的。

我可以通过这种方式获得原始结果。加上crosscovs = None,所以按原样计算,不是我们介绍的:

comp, crosscovs = kcca_working([X_sc, Y_sc], reg=0.,
         numCC=2, kernelcca=False, ktype='linear',
         gausigma=1.0, degree=2, crosscovs = None)
print(comp)
[array([[-0.13109264,  0.27302457],
       [ 0.02143325, -0.12466608],
       [-0.71655285, -0.2912628 ],
       [-0.43693303,  1.03824477],
       [ 0.01622265, -0.31573818],
       [ 0.49365965, -0.44098996]]), array([[ 0.2205752 ,  1.99601077],
       [-0.90927705, -1.79051045]])]

我打印了典型协变量 (comp),它们是 Bilenko 等人 2016 年方程式 4 中的 a' 和 b'。

a, b 和 a', b' 不同,只是比例不同,所以就所有目的而言,它们是等价的。这是因为相关定义。

为了证明这一点,让我从每个案例中提取数字并计算比率:

print(0.00061439/-0.00375779)
-0.16349769412340764
print(0.02143325/-0.13109264)
-0.16349697435340382

它们是相同的结果。

修改后,您可以在其顶部构建。

参考文献:

  1. 很酷 post,示例和解释在 Python,使用库 pyrccahttps://towardsdatascience.com/understanding-how-schools-work-with-canonical-correlation-analysis-4c9a88c6b913

  2. Bilenko、Natalia Y. 和 Jack L. Gallant。 “Pyrcca:python 中的正则化内核典型相关分析及其在神经影像学中的应用。”神经信息学前沿 10 (2016):49。解释 pyrcca 的论文:https://www.frontiersin.org/articles/10.3389/fninf.2016.00049/full