【问题标题】:Canonical correlation analysis on covariance matrices instead of raw data协方差矩阵而非原始数据的典型相关分析
【发布时间】:2021-06-13 22:23:38
【问题描述】:

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

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

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

mm colnames(mm)

您将数据集分为 x 和 y :
x y

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

我想知道是否有一种方法可以从以均值为中心的矩阵开始执行 cca: x = 比例(x,比例 = F) y = scale(Y, scale = 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

我希望足够详尽!

【问题讨论】:

  • 你想做什么?特别喜欢任何示例或链接?
  • 我用一个例子更新了原来的问题!
  • 真实的。也许这有帮助:towardsdatascience.com/…
  • 谢谢,但这只是一个使用普通算法处理原始数据的教程。我正在寻找适用于协方差矩阵的版本

标签: statistics correlation analysis


【解决方案1】:

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

总之:您需要根据具体情况制定一些方法。让我在下面为您展示一个示例。

什么是典型相关分析 (CCA)?

典型相关分析 (CCA):帮助您确定可以在两个数据集之间创建的最佳线性关系。见wikipedia。有关示例,请参见参考资料。我将按照post 获取数据并使用库。 设置库,上传数据,选择一些变量,移除 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 年 Eq4 中的 a 和 b。

我只是想让你注意这一点:

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

这正是该操作发生的地方。请注意,这并不完全是 Wikipedia 中的定义,但在数学上是等价的。

相关性的计算

我将像在维基百科中一样计算相关性,但稍后我将修改该函数,因此将提供一些细节,以确保这清楚地回答了原始问题。

# 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 的计算

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

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 年 Eq4 中的 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 年 Eq4 中的 a' 和 b'。

a, b 和 a', b' 是不同的,但它们只是在尺度上不同,所以它们是等价的。这是因为相关性定义。

为了表明让我从每个案例中提取数字并计算比率:

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

它们是相同的结果。

当它被修改后,你可以在它的顶部构建。

参考资料:

  1. 使用库pyrccahttps://towardsdatascience.com/understanding-how-schools-work-with-canonical-correlation-analysis-4c9a88c6b913

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

【讨论】:

  • 非常感谢您的回答。同时我自己在R中找到了解决方案,但是您的回答将来会对我有所帮助!
  • @cpt69 也可以发布您的答案
  • 当然,我只是使用了 R 库 CCA 中已经实现的函数 geigen。它接受输入协方差矩阵并给出两个系数作为输出。 Fronìm,这只是计算规范变量和载荷的问题
  • 主要问题是我认为 RGCCA 会使用相同的算法,但显然事实并非如此。因此,我遇到了与以前相同的问题。你知道在 python 中是否有可能从协方差矩阵中应用 RGCCA?再次感谢您的支持和帮助! ——
  • 我不知道RGCCA,但思路应该差不多。
猜你喜欢
  • 2019-04-15
  • 1970-01-01
  • 2019-04-10
  • 2016-05-26
  • 2017-02-12
  • 2020-07-04
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多