R中的典型相关分析

Question

我正在使用R(和CCA包)并试图用两个变量集(物种丰度和食物丰度分别存储为两个矩阵Y和X)进行正则化的典型相关分析,其中单位数(N = 15)小于矩阵中变量的数量,大于400(大多数是潜在的"解释"变量,只有12-13"响应"变量).冈萨雷斯等人.(2008年,http://www.jstatsoft.org/v23/i12/paper)请注意,该软件包"包含一个CCA的正则化版本,用于处理比单位更多变量的数据集",这当然只是我所拥有的15"单位." 因此,我正在尝试使用CCA包执行正则化的规范相关分析,以便查看我的变量集中的关系.我一直在关注Gonzalez等人(2008)在他们的论文中进行的过程.但是,我收到一条错误消息Error in chol.default(Bmat) : the leading minor of order 12 is not positive definite,我不知道这意味着什么或该做些什么.这是代码,任何有关该主题的想法或知识将不胜感激.

library(CCA)
correl <- matcor(X, Y)
img.matcor(correl, type = 2)
res.regul <- estim.regul(X, Y, plt = TRUE,
    grid1 = seq(0.0001, 0.2, l=51),
    grid2 = seq(0, 0.2, l=51))

Error in chol.default(Bmat) : the leading minor of order 12 is not positive definite

(注意:estim.regul()当您使用CCA的样本数据nutrimouse时,需要很长时间(~30-40分钟)才能完成).

有什么建议？有谁知道如何处理这个错误？是因为我的一些列中有一个NA吗？可能是因为有太多0的列？预先感谢您提供给这个综合统计数据和R新手的任何帮助.

Answer 1

背景

典型相关分析(CCA)是一种探索性数据分析(EDA)技术,提供在相同实验单元上收集的两组变量之间的相关关系的估计.通常,用户将有两个数据矩阵,X和Y,其中行代表实验单位,nrow(X)== nrow(Y).

在R中,基础包提供功能cancor()以启用CCA.这仅限于观察数量大于变量(特征)数量的情况,nrow(X)> ncol(X).

R包CCA是提供扩展CCA功能的几种之一.包CCA在cancor()周围提供了一组包装函数,可以考虑特征计数超过实验单位数ncol(X)> nrow(X)的情况.Gonzalez等人(2008)CCA:扩展典型相关分析的R包,详细描述了工作原理.CCA软件包 1.2版(2014-07-02发布)是撰写本文时的最新版本.

可能还值得一提的是kinship,accuracy早期答案中提到的软件包不再托管在CRAN上.

诊断

在跳转到其他软件包或对您的(可能是来之不易的!)数据应用未知方法之前,尝试和诊断数据问题可能是有益的.

传递给这里提到的任何CCA例程的矩阵理想情况下应该是数字完整的(没有缺失值).传递给这里提到的任何CCA例程的矩阵理想情况下应该是数字完整的(没有缺失值).通过该过程估计的规范相关的数量将等于X和Y的最小列等级,即<= min(ncol(X),ncol(Y)).理想情况下,每个矩阵的列将是线性独立的(而不是其他矩阵的线性组合).

例:

library(CCA)
data(nutrimouse)
X <- as.matrix(nutrimouse$gene[,1:10])
Y <- as.matrix(nutrimouse$lipid)

cc(X,Y) ## works

X[,1] <- 2 * X[,9] ## column 9 no longer provides unique information

cc(X,Y)

Error in chol.default(Bmat) :
  the leading minor of order 9 is not positive definite

这是原帖中的症状.一个简单的测试是在没有该列的情况下尝试拟合

cc(X[,-9],Y) ## works

因此,虽然从您从分析中删除数据的意义上这可能令人沮丧,但该数据无论如何都不提供信息.您的分析只能与您提供的数据一样好.

此外,有时可以通过?scale对输入矩阵中的一个(或两个)使用标准化(参见)变量来解决数值不稳定性:

X <- scale(X)

虽然我们在这里,但也许值得指出的是,正规化的CCA基本上是一个两步过程.进行交叉验证以估计正则化参数(使用estim.regul()),然后使用那些参数来执行正则化的CCA(with rcc()).

文中提供的示例(原始帖子中逐字逐句使用的参数)

res.regul <- estim.regul(X, Y, plt = TRUE,
                               grid1 = seq(0.0001, 0.2, l=51),
                               grid2 = seq(0, 0.2, l=51))

要求在51*51 = 2601单元网格上进行交叉验证.虽然这会为纸张生成漂亮的图形,但这些对于您自己的数据进行初始测试并不是明智的设置.正如作者所述,"计算要求不高.在51 x 51网格的'当前使用'计算机上,它持续不到一个半小时".自2008年以来,情况有所改善,但默认的5 x 5网格产生了

estim.regul(X,Y,plt=TRUE) 

是探索目的的绝佳选择.如果你要犯错误,你也可以尽快制作错误.

Answer 2

你的X方差 - 协方差矩阵不是正定的,因此内部调用时会出错fda::geigen.

有一个在一个类似的功能正则CCA mixOmics包,但我想它会导致同样的错误消息,因为它基本上使用相同的方法(除了它们插入的geigen直接功能,进入rcc功能).我实际上无法记住如何使用我的数据来解决相关问题(但是一旦我再次找到它,我会查看我的旧代码:-)

一种解决方案是使用广义Cholesky分解.在亲属关系中有一个(gchol;小心,它返回一个下三角矩阵)或precision(sechol)包.当然,这意味着修改函数内部的代码,但它确实不是问题,IMO.或者你可以尝试让VAR(X)PD与make.positive.definite来自corpcor包.

作为替代方案,您可以考虑使用RGCCA包,它提供了PLS(路径建模)和CCA方法与K个区块的统一接口.