Как использовать функцию princomp() в R, когда ковариационная матрица имеет нули?

При использовании функции princomp() в R возникает следующая ошибка: "covariance matrix is not non-negative definite".

Я думаю, это связано с тем, что некоторые значения равны нулю (фактически близки к нулю, но становятся равными нулю при округлении) в ковариационной матрице.

Есть ли способ продолжить работу с PCA, когда ковариационная матрица содержит нули?

[К вашему сведению: получение ковариационной матрицы является промежуточным этапом внутри вызова princomp(). Файл данных для воспроизведения этой ошибки можно скачать отсюда - http://tinyurl.com/6rtxrc3]


r statistics pca eigenvector princomp
person 384X21    schedule 19.12.2011    source источник
comment
Добавление образца ввода для воспроизводимости проблемы полезно для отвечающих.   -  person Richie Cotton    schedule 19.12.2011
comment
Если вы посмотрите на stats:::princomp.default, вы увидите, что ошибка возникает, когда у вас есть отрицательные собственные значения в ковариационной матрице.   -  person Richie Cotton    schedule 19.12.2011
comment
@ Ричи Коттон: Я хотел бы предоставить. Мои данные огромны (10K x 10K), и я не понял, что вызывает ошибку. Я буду рад узнать, есть ли способ извлечь тревожную часть данных и опубликовать ее здесь!   -  person 384X21    schedule 19.12.2011
comment
cv <- matrix(c(1, 2, 2, 1), nrow = 2); princomp(covmat = cv) воспроизводит ошибку. Не знаю, насколько это актуально для вашего набора данных.   -  person Richie Cotton    schedule 19.12.2011
comment
Спасибо ! Я смог воспроизвести его с небольшого фрагмента исходной матрицы, 1K x 1K (размер файла 5,5 МБ). Мне было интересно, как опубликовать это. Я уверен, что некоторые элементы ковариационной матрицы будут равны нулю (или близки к этому), так как мои входные данные содержат большие куски одинаковых значений.   -  person 384X21    schedule 19.12.2011
comment
Вы можете закинуть его на github или другой бесплатный хостинг данных. Google выдает многих, кто утверждает, что предлагает эту услугу, но у меня нет опыта работы с ними, поэтому я не могу так или иначе поручиться.   -  person Chase    schedule 19.12.2011
comment
@Chase: Спасибо, я выложил файл онлайн.   -  person 384X21    schedule 19.12.2011
comment
@user1029725 user1029725, прочитать файл с помощью read.table непросто. Сделайте что-то вроде «m = matrix(runif(4), nrow = 2, ncol = 2)», а затем «write.table(m, file = bla)».   -  person Paul Hiemstra    schedule 19.12.2011
comment
Обновлен входной файл. Это симметричная матрица, поэтому загружается только нижний треугольник!   -  person 384X21    schedule 19.12.2011
comment
Чтобы преобразовать треугольную матрицу в полную матрицу (как требуется для функции princomp()), я обычно сначала заменяю NA's на 0's, а затем добавляю транспонирование матрицы в исходную матрицу (что-то вроде этого - matrix <- matrix + t(matrix))   -  person 384X21    schedule 19.12.2011
comment
Похоже, ваша матрица может быть единственной, т. е. какой-то столбец представляет собой линейную комбинацию других. Нули в ковариационной матрице — это хорошо, но эти нули возникают из-за того, что столбцы данных независимы. Если, как вы говорите, ваши данные имеют большие участки одинаковых значений, то ваша ковариационная матрица будет иметь большие записи, а не маленькие. Возможно, вам придется удалить избыточные столбцы, прежде чем вы сможете использовать princomp.   -  person Wesley    schedule 19.12.2011
comment
@Wesley: Вполне возможно, что некоторые столбцы представляют собой линейную комбинацию других. Но удаление избыточных столбцов равносильно удалению переменной из анализа данных, что, я думаю, нежелательно. В любом случае, я попробую. Кстати, как удалить лишние столбцы?   -  person 384X21    schedule 19.12.2011
comment
У @DWin есть несколько хороших предложений в его ответе ниже.   -  person Wesley    schedule 19.12.2011

Ответы (1)


arrow_upward
9
arrow_downward

Первая стратегия может заключаться в уменьшении аргумента допуска. Мне кажется, что princomp не будет передавать аргумент допуска, но prcomp принимает аргумент 'tol'. Если это неэффективно, это должно идентифицировать векторы с почти нулевой ковариацией:

nr0=0.001
which(abs(cov(M)) < nr0, arr.ind=TRUE)

И это будет идентифицировать векторы с отрицательными собственными значениями:

which(eigen(M)$values < 0)

Используя пример h9 на странице справки (qr):

> which(abs(cov(h9)) < .001, arr.ind=TRUE)
      row col
 [1,]   9   4
 [2,]   8   5
 [3,]   9   5
 [4,]   7   6
 [5,]   8   6
 [6,]   9   6
 [7,]   6   7
 [8,]   7   7
 [9,]   8   7
[10,]   9   7
[11,]   5   8
[12,]   6   8
[13,]   7   8
[14,]   8   8
[15,]   9   8
[16,]   4   9
[17,]   5   9
[18,]   6   9
[19,]   7   9
[20,]   8   9
[21,]   9   9
> qr(h9[-9,-9])$rank  
[1] 7                  # rank deficient, at least at the default tolerance
> qr(h9[-(8:9),-(8:9)])$ take out only the vector  with the most dependencies
[1] 6                   #Still rank deficient
> qr(h9[-(7:9),-(7:9)])$rank
[1] 6

Другой подход может заключаться в использовании функции alias:

alias( lm( rnorm(NROW(dfrm)) ~ dfrm) )
person IRTFM    schedule 19.12.2011
comment
Хороший. Я не встречал alias раньше. - person Richie Cotton; 22.12.2011