Как рассчитать процентные доли вклада видов для веганских rda / cca объектов? - PullRequest
0 голосов
/ 04 мая 2018

Я пытаюсь воспроизвести столбец («переменная» в FactoMineR::PCA, «вид» в vegan::rda) процентное отношение к осям из пакета FactoMineR в vegan. Вклады кодируются в FactoMiner::PCA объектах:

library(FactoMineR)
library(vegan)

data(dune)

fm <- FactoMineR::PCA(dune, scale.unit = FALSE, graph = FALSE)

head(round(sort(fm$var$contrib[,1], decreasing = TRUE), 3))
# Lolipere Agrostol Eleopalu Planlanc  Poaprat  Poatriv 
#  17.990   16.020   13.866    7.088    6.861    4.850 

Посмотрев на код для FactoMiner::PCA, я обнаружил, что вклады рассчитываются как квадратные координаты оси, разделенные на собственное значение оси и умноженные на 100%:

head(round(sort(100*fm$var$coord[,1]^2/fm$eig[1], decreasing = TRUE), 3))
# Lolipere Agrostol Eleopalu Planlanc  Poaprat  Poatriv 
#  17.990   16.020   13.866    7.088    6.861    4.850 

Я не могу воспроизвести вышеупомянутые вычисления, используя vegan::rda объект:

vg <- rda(dune)

head(round(sort(100*scores(vg, choices = 1, display = "sp", 
scaling = 0)[,1]^2/vg$CA$eig[1], decreasing = TRUE), 3))
# Lolipere Agrostol Eleopalu Planlanc  Poaprat  Poatriv 
#   0.726    0.646    0.559    0.286    0.277    0.196

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

# vegan eigenvalue for axis 1
vg$CA$eig[1]
#     PC1 
# 24.79532 

# FactoMineR eigenvalue for axis 1
fm$eig[1]
# [1] 23.55555

# vegan column coordinates for axis 1
head(round(scores(vg, choices = 1, display = "sp", scaling = 0)[,1], 3))
# Achimill Agrostol Airaprae Alopgeni Anthodor Bellpere 
#  -0.176    0.400    0.007    0.155   -0.163   -0.097 

#FactoMineR, column coordinates for axis 1
head(round(fm$var$coord[,1], 3))
# Achimill Agrostol Airaprae Alopgeni Anthodor Bellpere 
#   0.854   -1.943   -0.033   -0.751    0.791    0.472 

# Sum of column coordinates for vegan axis 1 to illustrate the difference
sum(scores(vg, choices = 1, display = "sp", scaling = 0)[,1])
# [1] -0.796912

# Sum of column coordinates for FactoMineR axis 1 to illustrate the difference
sum(fm$var$coord[,1])
# [1] 3.867738

Как рассчитать процентное отношение столбцов / видов к осям ординации с использованием объекта vegan rda ?

Ответы [ 2 ]

0 голосов
/ 05 мая 2018

Немасштабированные оценки в vegan не масштабируются в (нормальном) смысле, что их сумма квадратов равна 1 - независимо от собственных значений:

> colSums(scores(vg, choices=1:4,dis="sp", scaling=0)^2)
  PC1 PC2 PC3 PC4 
    1   1   1   1 

Я думаю, что это задокументировано. Если вы хотите назвать эти квадраты в качестве вклада, это нормально для меня. То же самое относится к cca, но там вам нужно изучить взвешенные суммы квадратов. Более того, немасштабированные баллы для сайтов (dis = "si") будут иметь одинаковую единицу суммы квадратов: это идея быть немасштабированной. Если вы масштабируете либо виды, либо участки, то те же отношения больше не действуют для другого набора баллов. Как правило, отсутствие немасштабирования означает, что оценки являются ортонормированными, так что их перекрестный продукт представляет собой единичную матрицу (диагональ или сумма квадратов 1 и недиагональных элементов 0). Для масштабированных баллов эти суммы квадратов пропорциональны собственным значениям (но может быть полезно прочитать веганский виньетка по проектным решениям для const муравейного масштабирования баллов).

веганский функции goodness и inertcomp могут (или не могут) дать вам информацию, которую вы ищете.

0 голосов
/ 04 мая 2018
Координаты столбца немасштабированного (scaling = 0)

vegan имеют равные суммы квадратов (то есть 1 для каждой оси). Вы можете получить «вклад столбцов», просто возведя в квадрат немасштабированные координаты:

head(sort(round(100*scores(vg, display = "sp", scaling = 0)[,1]^2, 3), decreasing = TRUE))
# Lolipere Agrostol Eleopalu Planlanc  Poaprat  Poatriv 
#  17.990   16.020   13.866    7.088    6.861    4.850

Сумма этих «взносов» равна 100%, как указано выше:

sum(100*scores(vg, display = "sp", scaling = 0)[,1]^2)
# [1] 100
...