Как добавить цвета и легенду к ordiellipse в базовом пакете в R?

Question

У меня есть датафрейм (сайт по видам матрицы), который выглядит так:

            SP1      SP2     SP3     SP4
    US       5       6       2       5
    US       5       6       2       5
    UK       5       6       2       5
   AUS       5       6       2       5

Я пытаюсь создать график PCoA (анализ основных координат) с полигонами / эллипсами с вероятностью 95%. Мне нужно однозначно раскрасить каждую страну (точки) цветом, а каждый эллипс должен иметь соответствующий цветовой код страны и легенды.

#My current code
#If you  need a dataframe
df <- t(data.frame(matrix(rexp(10000, rate=10),nrow=100,ncol=100)))
rownames(df) <- rep(c("UK", "US", "Aus","Spain"),length.out=100)#I dont know how to loop this over 100 times to make it the rownames


df <- as.matrix(df[,-1]) #Use this to convert dataframe to matrix
row.names(df) <- df[,1]#Use this to convert dataframe to matrix
dat <- df
dat.db <- vegdist(dat, method = "bray")
dat.pcoa <- cmdscale(dat.db, eig = TRUE, k = 3)
explainvar1 <- round(dat.pcoa$eig[1] / sum(dat.pcoa$eig), 3) * 100
explainvar2 <- round(dat.pcoa$eig[2] / sum(dat.pcoa$eig), 3) * 100
explainvar3 <- round(dat.pcoa$eig[3] / sum(dat.pcoa$eig), 3) * 100
sum.eig <- sum(explainvar1, explainvar2, explainvar3)

# Define Plot Parameters
par(mar = c(5, 5, 1, 2) + 0.1)
# Initiate Plot
plot(dat.pcoa$points[ ,1], dat.pcoa$points[ ,2],
     xlab = paste("PCoA 1 (", explainvar1, "%)", sep = ""),
     ylab = paste("PCoA 2 (", explainvar2, "%)", sep = ""),
     pch = 16, cex = 2.0, type = "n", cex.lab = 1.5, cex.axis = 1.2, axes = FALSE)
axis(side = 1, labels = T, lwd.ticks = 2, cex.axis = 1.2, las = 1)
axis(side = 2, labels = T, lwd.ticks = 2, cex.axis = 1.2, las = 1)
abline(h = 0, v = 0, lty = 3)
box(lwd = 2)
points(dat.pcoa$points[ ,1], dat.pcoa$points[ ,2],
       pch = 19, cex = 1, bg = "gray", col = "grey")
ordiellipse(ord = dat.pcoa, groups = rownames(dat), kind = "se",conf = .95,col = NULL)

Примечание: это не тот же вопрос, что и здесь. Этот вопрос только спрашивает, как выполнить ordiplot в базовом пакете (так как я ударил стену для ggplot2)

Answer 1

Для базового сюжета вы можете указать вектор цветов.

Следовательно, для каждой точки вы должны назначить цвет и использовать этот вектор цветов длиной n точек в качестве входных данных для параметра col в points, Самый простой способ сделать это - проанализировать имена строк стран в данных PCA и использовать целочисленные значения коэффициента для индексации определенной цветовой палитры.

Точно так же вы можете предоставить вектор цветов длины n групп для ordiellipse

Вот код, основанный на вашем примере df:

library(vegan)
df <- t(data.frame(matrix(rexp(10000, rate=10),nrow=100,ncol=100)))
rownames(df) <- rep(c("UK", "US", "Aus","Spain"), length.out=100)#I dont know how to loop this over 100 times to make it the rownames
colnames(df) <- paste0("SP", 1:ncol(df))

Теперь мы учитываем названия стран

# factor country and set levels
col_vector <- factor(rownames(dat.pcoa$points), levels=unique(rownames(dat.pcoa$points)))
col_vector

# see integer values of your factored countries and the given order
str(col_vector) 
# Factor w/ 4 levels "UK","US","Aus",..: 1 2 3 4 1 2 3 4 1 2 ...

Создайте цветовую палитру и индексируйте, используя факторные названия стран

# palette() is the default R color palette or you can specify a vector of 4 colors
col_palette <- palette()[col_vector]
col_palette

Затем используйте эти значения в своем сюжете

par(mar = c(5, 5, 1, 2) + 0.1)
# Initiate Plot
plot(dat.pcoa$points[ ,1], dat.pcoa$points[ ,2],
     xlab = paste("PCoA 1 (", explainvar1, "%)", sep = ""),
     ylab = paste("PCoA 2 (", explainvar2, "%)", sep = ""),
     pch = 16, cex = 2.0, type = "n", cex.lab = 1.5, cex.axis = 1.2, axes = FALSE)
axis(side = 1, labels = T, lwd.ticks = 2, cex.axis = 1.2, las = 1)
axis(side = 2, labels = T, lwd.ticks = 2, cex.axis = 1.2, las = 1)
abline(h = 0, v = 0, lty = 3)
box(lwd = 2)
points(dat.pcoa$points[ ,1], dat.pcoa$points[ ,2],
   pch = 19, cex = 1, bg = "gray", col = col_palette)
ordiellipse(ord = dat.pcoa, groups = col_vector, kind = "se", conf = .95, col = unique(col_palette)) 

Здесь мы проверяем, что метки соответствуют назначенному цвету

# sanity check that point colors match label
text(dat.pcoa$points[ ,1], dat.pcoa$points[ ,2], labels=rownames(dat.pcoa$points), col = col_palette) 
# sanity check that ellipse color match labels
ordiellipse(ord = dat.pcoa, groups = col_vector, kind = "se", conf = .95, col = unique(col_palette), label=TRUE) 

Наконец добавить легенду

legend('topright', legend=unique(col_vector), col=unique(col_palette), pch = 16)