Канонический корреляционный анализ в R

Question

Я использую R (и пакет CCA) и пытаюсь выполнить регуляризованный канонический корреляционный анализ с двумя переменными наборами (численность видов и численность пищи, сохраняемые как две матрицы Y и X соответственно), в которых число единиц (N= 15) меньше, чем число переменных в матрицах, которое составляет> 400 (большинство из них являются потенциальными «объяснительными» переменными, причем только 12–13 «ответных» переменных).Гонсалес и соавт.(2008, http://www.jstatsoft.org/v23/i12/paper) обратите внимание, что пакет «включает в себя упорядоченную версию CCA для работы с наборами данных с большим количеством переменных, чем единиц», что, безусловно, то, что у меня есть только с 15 «единицами». Таким образом, япытаясь выполнить регуляризованный канонический корреляционный анализ с использованием пакета CCA, чтобы посмотреть на отношения в моих наборах переменных. Я следил за процессом, который Гонсалес и др. (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() занимаетдолго (~ 30-40 минут), чтобы завершить, когда вы используете образец данных нутримы из CCA).

Любой совет? Кто-нибудь знает, что делать с этой ошибкой? Это потому, что некоторые из моих столбцов имеютNA в них? Может ли это быть из-за столбцов со слишком многими нулями? Заранее благодарим за любую помощь, которую вы можете предложить этой объединенной статистике и новичку R.

Answer 1

Фон

Канонический корреляционный анализ (CCA) - это метод анализа поисковых данных (EDA), обеспечивающий оценки корреляционных отношений между двумя наборами переменных, собранных в одних и тех же экспериментальных единицах. Как правило, пользователи имеют две матрицы данных, X и Y, где строки представляют экспериментальные единицы, nrow (X) == nrow (Y).

В R базовый пакет предоставляет функцию cancor () для включения CCA. Это ограничено случаями, когда количество наблюдений больше, чем количество переменных (признаков), nrow (X)> ncol (X).

Пакет CCA R - один из нескольких, которые предоставляют расширенные функциональные возможности CCA. Пакет CCA предлагает набор функций-оболочек для cancor (), которые позволяют рассматривать случаи, когда число объектов превышает число экспериментальных единиц, ncol (X)> nrow (X). Gonzalez et al. (2008) CCA: пакет R для расширения канонического корреляционного анализа , описывает работу в некоторых деталях. Версия 1.2 пакета CCA (опубликована 2014-07-02) актуальна на момент написания.

Вероятно, также стоит упомянуть, что пакеты kinship и accuracy, упомянутые в предыдущем ответе, больше не размещаются в CRAN.

Диагностика

Прежде чем переходить к другим пакетам или применять неизвестные методы к вашим (предположительно с трудом завоеванным!) Данным, возможно, полезно попытаться диагностировать, в чем может быть проблема с данными.

Матрицы, передаваемые любой из подпрограмм CCA, упомянутых здесь, в идеале должны быть численно полными (без пропущенных значений). Матрицы, передаваемые любой из подпрограмм CCA, упомянутых здесь, в идеале должны быть численно полными (без пропущенных значений). Количество канонических коррелятов, оцененных по этой процедуре, будет равно минимальному рангу столбца X и Y, то есть <= min (ncol (X), ncol (Y)). В идеале столбцы каждой матрицы должны быть линейно независимыми (а не линейными комбинациями других). </p> * * 1 022 Пример: * 1 023 *

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) 

является отличным выбором для исследовательских целей. Если вы собираетесь делать ошибки, вы можете сделать их как можно быстрее.

Answer 2

Ваша X-дисперсионно-ковариационная матрица не является положительно определенной, поэтому ошибка при внутреннем вызове fda::geigen.

Существует аналогичная функция для регуляризованного CCA в пакете mixOmics , но я предполагаю, что это приведет к тому же сообщению об ошибке, потому что оно в основном использует тот же подход (за исключением того, что они подключили функцию geigen напрямую в функцию rcc). Я не могу вспомнить, как заставить его работать с моими данными для связанной проблемы (но я посмотрю свой старый код, как только найду его снова: -)

Одним из решений было бы использование обобщенного разложения Холецкого. Есть один в родстве (gchol; будьте осторожны, он возвращает нижнюю треугольную матрицу) или точность (sechol). Конечно, это означает изменение кода внутри функции, но это не проблема, IMO. Или вы можете попытаться сделать Var (X) PD с make.positive.definite из пакета corpcor .

В качестве альтернативы вы можете рассмотреть возможность использования пакета RGCCA , который предлагает унифицированный интерфейс для методов PLS (моделирование пути) и CCA с k блоками.