Простая структура. Сколько факторов её составляет?

индуцировано обсуждением вопроса в  r-statistics

В СССР в 80-х сборники "братские могилы" по анализу данных буквально пронизывали две идеи статистиков:

1. Как сделать критерии для нелинейных задач.
2. Как сделать критерии для многомерного анализа данных.

Между тем господин Эфрон уже вовсю вел пропаганду применения бутстреп-оценки. Что признают даже его критики (наверное оставшись без двух "вечнозеленых" идей :) ), особенно полезен бутстреп именно в областях статистики, где нет разработанных методов получения классической оценки. На самом деле бутстреп полезен и для проверки наличия ошибок применения сложной (или новой для исследователя) процедуры получения классической оценки.

Попробуем применить бутстреп в многомерном анализе данных. Для анализа наличия простой структуры и её объема воспользуемся набором данных "массив Ревелле" из исходной дискуссии.

data(bfi)

# чистим данные от неполных случаев
data <- na.omit(bfi) 

# получаем коррелированные перевыборки из исходной 

data.boot<- replicate(1000, 
                      prcomp(data[sample(1:nrow(data), 
                                  size=nrow(data), 
                                  replace=T),],
                             center=T,
                             scale.=T)$sdev)


# получаем некоррелированные перевыборки из исходной

data.boot.random <- replicate(1000, 
                              prcomp(apply(data,2, 
                                           function(x) x[sample(1:length(x),
                                                                size=length(x),
                                                                replace=T)]),
                                     center=T,
                                     scale.=T)$sdev)

# объединяем оба набора собственных значений принципиальных компонент

data.boot<- cbind(data.boot, data.boot.random)


# преобразуем в "длинный" формат, для построения графика

library(reshape)
dimnames(data.boot)<-list(1:length(data.boot[,1]),
                          paste("v", 1:length(data.boot[1,]), sep=""))
data.melt<- melt.matrix(data.boot, 
                        varnames=names(dimnames(data.boot)))


# строим график всех перевыборок собственных значений

library(ggplot2)
pcp<- ggplot(data.melt, 
             aes(X1,
                 value, 
                 group=X2))
pcp+geom_line(colour=alpha("black",1/2))

Как видно 6-ю компоненту выделять уже нельзя.


TODO: провести рефакторинг кода; заменить prcomp() на gmodels::fast.prcomp если это актуально; replicate() на snow версию