我正在使用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
(注意:当您使用来自CCA的示例数据营养品时,
estim.regul()
需要很长时间(〜30-40分钟)才能完成)。有什么建议吗?有人知道该错误怎么办吗?是否因为我的某些专栏文章中包含NA?可能是因为列的0太多了吗?预先感谢您可以为该统计和R新手提供的任何帮助。
最佳答案
背景
规范相关分析(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: An R Package to Extend Canonical Correlation Analysis详细描述了工作原理。 package CCA的1.2版(于2014年7月2日发布)在撰写本文时为最新版本。
可能还值得一提的是,先前答案中提到的kinship
和accuracy
软件包不再托管在CRAN上。
诊断
在跳到其他软件包或对您的数据(可能来之不易!)应用未知方法之前,尝试诊断数据问题可能是有益的。
理想情况下,传递给此处提到的任何CCA例程的矩阵在数值上都应该完整(无缺失值)。理想情况下,传递给此处提到的任何CCA例程的矩阵在数值上都应该完整(无缺失值)。该过程估计的规范相关数将等于X和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(使用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)
是用于探索性目的的绝佳选择。如果您要犯错误,则最好尽快使它们出错。
关于r - R中的典型相关分析,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/5850763/