факторный анализ - 2

Для лучшего понимания факторного анализа очень рекомендуется владеть матричной алгеброй. Вручную выполнять матричные операции вряд ли под силу гуманитарию без ущерба для душевного здоровья, поэтому на помощь могут придти такие программы, как R. Поясню, что означают выражения, упомянутые в конце предыдущего поста. Функция c() ("concatenate"), объединяет несколько элементов в единую структуру, например, вектор. Функция cbind ("column bind") объединяет переданные ей векторы в матрицу. Так что выражение A <- cbind(c(0.6, 0.7, 0.8, 0.1, 0.2, 0.2), c(0.1, 0.1, 0.2, 0.8, 0.7, 0.6)) сначала создает два вектора действительных чисел, записывает их как столбцы в матрицу и присваивает результат матрице с именем A.
Второе выражение R <- A %*% t(A) использует функцию t() ("transpose"), которая переворачивает матрицу нагрузок А. Оператор %*% выполняет матричное умножение А на перевернутую А. В результате получается воспроизведенная матрица корреляций.

Теперь наступает ответственный момент: в реальном исследовании нам известна корреляционная матрица, а матрицу нагрузок (т.е. матрицу коэффициентов влияния факторов на переменные) необходимо найти на ее основе. Однако тут возникает проблема: обычная корреляционная матрица немного не похожа на R в нашем уравнении: на ее диагонали находятся единицы (корреляция переменной с самой собой), а в матрице R -- так называемые "общности". Общностью называется доля дисперсии переменной, которая объясняется общими с другими переменными факторами. Вторая часть дисперсии переменной, дополняющая общность до 1, называется в факторном анализе "уникальностью" -- это часть, объясняемая специфическим для данной переменной фактором и ошибкой ее измерения. Таким образом, чтобы приступить к факторному анализу, придется сначала решить проблему оценки общностей. Обычно стат. пакеты скрывают от пользователя этот этап, предлагая решить проблему методом по умолчанию (чаще всего в качестве оценок вычисляется квадрат множественной корреляции каждой переменной со всеми остальными в анализе). Есть и другой подход: факторизовать обычную, "нередуцированную", корреляционную матрицу, оставить некоторое число факторов и оценить общности по матрице нагрузок. Для вычисления оценки общности следует взять сумму квадратов нагрузок переменной на выделенные факторы. Затем оценки общностей подставляют на диагональ корреляционной матрицы и вновь выделяют факторы и т.д. Процесс завершается тогда (если завершается, конечно), когда оценки общностей перестанут изменяться с заданной точностью.
Каким же образом происходит собственно выделение факторов? Для этого чаще всего используется процедура вычисления так называемых "собственных значений" и "собственных векторов" корреляционной матрицы. Вектор и собственное число идут в паре и имеют интересное свойство: если матрицу умножить на вектор, то результат равне произведению собственного числа на вектор. Получается, что собственное число матрицы в каком-то смысле эквивалентно самой матрице, что удивительно. Одна из теорем матричной алгебры утверждает, что для симметричной матрицы выполняется такое условие: ее можно полностью воспроизвести на основе собственных векторов и собственных чисел. Если каждый элемент собственного вектора умножить на квадратный корень из собственного числа (обозначим такую матрицу буквой A), то исходную матрицу R можно будет получить djn nfr^ R <- A %*% t(A)
В R за извлечение собственных чисел и векторов отвечает функция eigen: ей передают матрицу, а она возвращает список из двух объектов -- матрицы vectors и вектора values. Итак, поместим на диагональ матрицы R единицы, чтобы она не отличалась от обычной корреляционной матрицы, и вычислим ее собственные числа и собственные вектора:

diag(R) <- 1
E <- eigen(R)

Теперь в Е хранится вот что:

$values
[1] 2.6561454 1.3963903 0.5810826 0.5552224 0.4113824 0.3997770

$vectors
           [,1]       [,2]        [,3]        [,4]        [,5]
[1,] -0.3653408 -0.4289664  0.79780245  0.03791813 -0.06249185
[2,] -0.3925694 -0.4407696 -0.56652385 -0.04075957 -0.19122656
[3,] -0.4518788 -0.3642986 -0.19311020 -0.03404948  0.15524607
[4,] -0.4034098  0.4665709  0.05410010 -0.24441286 -0.73414856
[5,] -0.4270119  0.3805622  0.02850586 -0.50727816  0.60964745
[6,] -0.4039088  0.3556410 -0.03912838  0.82381595  0.15742206
            [,6]
[1,]  0.20169830
[2,]  0.54077092
[3,] -0.77494505
[4,] -0.13391607
[5,]  0.20741148
[6,]  0.07342852

Умножим каждый собственный вектор на корень из соответствующего ему собственного числа, а результат умножения запишем в матрицу X:

X <- E$vectors %*% diag(sqrt(E$values))
V <- X[,c(1,2)]

Матрица Х и есть матрицей нагрузок, только факторов в ней столько же, сколько исходных переменных. Какой смысл оставлять такое количество факторов, если нашей целью было сокращение информации? Оставим только самые важные факторы. На важность фактора указывает величина собственного числа: чем оно больше, тем важнее фактор для объяснения корреляций. Во второй строке приведенного выше набора команд как раз и выполняется "обрезание" неважных факторов: в матрицу V копируется только два первых столбца. Обратите внимание, как хитро записывается это действие -- в квадратных скобках после имени матрицы указываются индексы элементов, которые нам нужны. Индекс для строк пропущен, поэтому берутся все строки матрицы, а индекс для столбца -- это объединение чисел 1 и 2. Теперь можно вычислить оценки общностей и результат занести на диагональ корреляционной матрицы:

  com = rowSums(V^2)
  diag(R) = com

Чтобы эту последовательность действий не пришлось выполнять много раз подряд вручную, поместим ее внутри цикла, выполняющегося, скажем, сто раз (для контроля в цикле выводится вектор общностей):

for (i in 1:100) {
  E = eigen(R)
  X = E$vectors %*% diag(sqrt(E$values))
  V = X[,c(1,2)]
  com = rowSums(V^2)
  print(com)
  diag(R) = com
}

При выполнении этого блока в R на каком-то шаге замечаем, что оценки общностей перестали меняться и равны 0.37 0.50 0.68 0.65 0.53 0.40. Как видим, нам удалось в точности воспроизвести диагональ нашей первоначальной матрицы R. Матрица V после выполнения цикла будет содержать нагрузки переменных на два выделенных фактора. Только здесь возникает новая проблема: матрица V совсем не похожа на матрицу A, с которой мы начинали!

           [,1]       [,2]
[1,] -0.4991478 -0.3476370
[2,] -0.5706950 -0.4175011
[3,] -0.7121064 -0.4158179
[4,] -0.6304597  0.5025142
[5,] -0.6321430  0.3611029
[6,] -0.5622789  0.2895556