Параметры форматирования рисунка PCA в R - изменение отдельных точек на групповые средние/SE

Я попробовал этот вопрос в stats.stackexchange, и кто-то предложил мне попробовать его здесь, так что вот:

Я завершил PCA-анализ в R с пакетом VEGAN некоторых экологических данных о здоровье деревьев. Всего насчитывается 80 деревьев (то есть 80 «участков»), разделенных на четыре категории обработки. У меня есть данные, представленные точками с цветовой кодировкой — цвета в соответствии с группами лечения. Вместо того, чтобы строить отдельные сайты/деревья на двойном графике PCA, я хотел бы сделать что-то вроде графика в виде прямоугольника и усов, который имеет четыре «креста», которые показывают центр тяжести для каждой группы и SE в обоих измерениях PCA. Я видел подобные цифры в газетах, но не могу найти R-скрипт для такого построения. Какие-либо предложения? (Я хотел бы опубликовать здесь пример изображения того, что я ищу, но те, которые я могу найти, все платные, извините).

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

Код, который я запускал, очень прост:

p1<-princomp(scale(health, scale=T))
summary(p1)
scores(p1)
plot(p1)
loadings(p1)
biplot(p1, xlab = "PC 1 (38%)", ylab = "PC 2 (22%)",cex=0.6)
plot(p1$scores[,1],p1$scores[,2])
names(p1)

plot(p1$scores[,1],p1$scores[,2], type='n', xlab="PC I", ylab="PC II")
text(p1$scores[,1],p1$scores[,2] labels=Can$tree)

r visualization pca vegan
person Jas Max    schedule 21.03.2014    source источник
comment
Вы понимаете, что princomp() не в веганском, да? и даже не рекомендуется его использовать (предпочтительнее prcomp(), но даже это не в веганском).   -  person Gavin Simpson    schedule 21.03.2014
comment
Да, извините за путаницу — я просто показываю код, с которым я играл сегодня утром. Любые предложения по лучшему пакету/скрипту для того, что я пытаюсь сделать? Или, в том же духе, ресурс для различий между princomp() и prcomp()? Я разбираюсь в статистике вещей, но мои знания R несовершенны и в основном состоят из вещей, которые мне показывали другие люди. Я могу найти много документации, чтобы прочитать об отдельных пакетах, но мне бы хотелось получить хороший ресурс для «общей картины», например, когда и почему использовать один сценарий по сравнению с другим. Спасибо за ваше время - жаль тратить его.   -  person Jas Max    schedule 21.03.2014
comment
Спасибо за разъяснения; Я почти закончил ответ, используя веганские функции.   -  person Gavin Simpson    schedule 21.03.2014
comment
Что касается princomp и prcomp, то в первом используется разложение по собственным числам, а во втором - разложение по сингулярным числам (SVD). Последний, как правило, более надежен и также предлагает решение, когда переменных/столбцов больше, чем выборок/строк. Функция rda() в vegan также использует SVD для вычисления PCA.   -  person Gavin Simpson    schedule 21.03.2014

Ответы (1)


arrow_upward
0
arrow_downward

Я бы, наверное, начал с ordiellipse и посмотрел, подходит ли это вам.

### Reproducible example
require("vegan")
data(varespec)
data(varechem)

pca <- rda(varespec, scale = TRUE)
grp <- with(varechem, cut(Baresoil, 4, labels = 1:4))
cols <- c("red","orange","blue","forestgreen")
scl <- 1 ## scaling

plot(pca, display = "sites", scaling = scl, type = "n")
points(pca, display = "sites", scaling = scl, col = cols[grp], pch = 16)
lev <- levels(grp)
for (i in seq_along(lev)) { ## draw ellipse per group
  ordiellipse(pca, display = "sites", kind = "se", scaling = scl,
              groups = grp, col = cols[i], show.groups = lev[i])
}
## centroids
scrs <- as.data.frame(scores(pca, display = "sites", scaling = scl, 
                             choices = 1:2))
cent <- do.call(rbind, lapply(split(scrs, grp), colMeans))
points(cent, col = cols, pch = 3, cex = 1.1)

Это производит

введите здесь описание изображения

Вы можете удалить строку points() из приведенного выше кода, чтобы остановить отрисовку фактических образцов, но я подумал, что здесь это поучительно с точки зрения понимания того, что делает ordiellipse.

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

Однако при этом не будут учитываться ковариации оценок по двум осям. Обратите внимание на приведенный ниже пример, как в тех эллипсах, которые ориентированы под углами к осям, планки стандартных погрешностей не пересекаются с эллипсом в своих экстремальных точках. Если бы вы нарисовали рамку, содержащую область, определяемую планками погрешностей, она содержала бы эллипс, но это дает совсем другое представление о неопределенности в центроиде.

serrFun <- function(df) {
  apply(df, 2, function(x) sd(x) / sqrt(length(x)))
}
serr <- do.call(rbind, lapply(split(scrs, grp), serrFun))
for (i in seq_along(lev)) {
    arrows(cent[i, 1] - serr[i, 1], cent[i, 2],
           cent[i, 1] + serr[i, 1], cent[i, 2],
           col = cols[i], code = 3, angle = 90, length = 0.05)
    arrows(cent[i, 1], cent[i, 2] - serr[i, 2],
           cent[i, 1], cent[i, 2] + serr[i, 2], 
           col = cols[i], code = 3, angle = 90, length = 0.05)
}

введите здесь описание изображения

person Gavin Simpson    schedule 21.03.2014
comment
Это здорово, спасибо за подробный ответ. Я ценю то, что вы говорите об эллипсе, лучше представляющем ковариационную структуру. - person Jas Max; 21.03.2014