2016-05-28 6 views
2

Я пытаюсь визуализировать результаты PCoA{ape} путем создания биполя в R. Теперь оси получают метки меток по умолчанию 1 и ось 2, но я хочу отредактировать это.Изменение меток осей для biplot() в R

Это код, который я пробовал:

biplot(pcoa.ntK, Y=NULL, plot.axes=c(1,2), rn=ntnames, 
          xlabs="PC1 (%)", ylabs="PC2 (%)") 

Но метки не изменяются. Может кто-нибудь сказать мне, что я делаю неправильно здесь? И я также хотел бы отредактировать название, все советы для этого?

Мои данные:

ntK <- matrix( 
    c(0.00000, 0.01500, 0.01832, 0.02061, 0.01902, 0.01270, 0.02111, 0.01655, 0.01520, 0.01691, 
    0.01667, 0.00000, 0.01175, 0.01911, 0.01759, 0.01127, 0.01854, 0.01041, 0.00741, 0.02007, 
    0.02432, 0.01404, 0.00000, 0.02551, 0.01972, 0.01838, 0.02505, 0.01484, 0.01391, 0.02687, 
    0.01501, 0.01252, 0.01399, 0.00000, 0.01442, 0.01294, 0.01402, 0.01132, 0., 0.01455, 
    0.02343, 0.01951, 0.01830, 0.02440, 0.00000, 0.01727, 0.02470, 0.02021, 0.01699, 0.02482, 
    0.01320, 0.01054, 0.01439, 0.01847, 0.01457, 0.00000, 0.01818, 0.01366, 0.00977, 0.01394, 
    0.02468, 0.01950, 0.02206, 0.02251, 0.02343, 0.02040, 0.00000, 0.02028, 0.01875, 0.02558, 
    0.02254, 0.01276, 0.01522, 0.02117, 0.02234, 0.01790, 0.02363, 0.00000, 0.01152, 0.02557, 
    0.01804, 0.00792, 0.01244, 0.02019, 0.01637, 0.01116, 0.01904, 0.01004, 0.00000, 0.02099, 
    0.01862, 0.01988, 0.02227, 0.02200, 0.02218, 0.01476, 0.02408, 0.02066, 0.01947, 0.00000), 
    nrow=10, 
    ncol=10) 

library(ape) 
ntnames <- c("A","B","C","D","E","F","G","H","I","J") 
pcoa.ntK <- pcoa(ntK) 
+0

как насчет данных 'rednt'? – 989

+1

Кажется, что 'rednt' должен быть' pcoa.ntK', а 'namesaa' должен быть' ntnames' на вызове 'biplot'. – Molx

+0

Как было сказано ниже, заголовок оси жестко закодирован, где они берут свои имена из названий столбцов 'pcoa.ntK $ vectors'. Поэтому быстрое решение - переименовать это. Так что если вы планируете первые два компьютера. используйте 'colnames (векторы pcoa.ntK $) [1: 2] <- c (" Mytitle1 "," Mytitle2 ")', затем зарисуйте, как раньше. – user20650

ответ

0

biplot является родовой функцией. Метод по умолчанию и метод для использования с объектами, которые исходят от использования функции prcomp в пакете stats, позволяют вам указывать метки осей и заголовок, но по какой-то причине человек, который написал метод, который вызывается с объектами класса pcoa не дал вам указать их. Я думаю, что ваш единственный вариант - написать свою собственную версию biplot.pcoa (или попросить сопровождающего пакета добавить эту опцию).

Это очень быстрый и грязный взлом функции в пакете ape, который может делать то, что вы хотите, но никаких обещаний, что он не сломал что-то еще!

biplot.pcoa <- function (x, Y = NULL, plot.axes = c(1, 2), dir.axis1 = 1, dir.axis2 = 1, 
      rn = NULL, xlabs = NULL, ylabs = NULL, main = NULL, ...) 
{ 
    k <- ncol(x$vectors) 
    if (k < 2) 
    stop("There is a single eigenvalue. No plot can be produced.") 
    if (k < plot.axes[1]) 
    stop("Axis", plot.axes[1], "does not exist.") 
    if (k < plot.axes[2]) 
    stop("Axis", plot.axes[2], "does not exist.") 
    if (!is.null(rn)) 
    rownames(x$vectors) <- rn 
    labels = colnames(x$vectors[, plot.axes]) 
    if (!is.null(xlabs)) labels[1] <- xlabs 
    if (!is.null(ylabs)) labels[2] <- ylabs 
    diag.dir <- diag(c(dir.axis1, dir.axis2)) 
    x$vectors[, plot.axes] <- x$vectors[, plot.axes] %*% diag.dir 
    if (is.null(Y)) { 
    limits <- apply(x$vectors[, plot.axes], 2, range) 
    ran.x <- limits[2, 1] - limits[1, 1] 
    ran.y <- limits[2, 2] - limits[1, 2] 
    xlim <- c((limits[1, 1] - ran.x/10), (limits[2, 1] + 
              ran.x/5)) 
    ylim <- c((limits[1, 2] - ran.y/10), (limits[2, 2] + 
              ran.y/10)) 
    par(mai = c(1, 1, 1, 0.5)) 
    plot(x$vectors[, plot.axes], xlab = labels[1], ylab = labels[2], 
     xlim = xlim, ylim = ylim, asp = 1) 
    text(x$vectors[, plot.axes], labels = rownames(x$vectors), 
     pos = 4, cex = 1, offset = 0.5) 
    if (is.null(main)){ 
     title(main = "PCoA ordination", line = 2.5) 
    } else title(main = main, line = 2.5) 
    } 
    else { 
    n <- nrow(Y) 
    points.stand <- scale(x$vectors[, plot.axes]) 
    S <- cov(Y, points.stand) 
    U <- S %*% diag((x$values$Eigenvalues[plot.axes]/(n - 
                 1))^(-0.5)) 
    colnames(U) <- colnames(x$vectors[, plot.axes]) 
    par(mai = c(1, 0.5, 1.4, 0)) 
    biplot(x$vectors[, plot.axes], U, xlab = labels[1], ylab = labels[2]) 
    if (is.null(main)) { 
    title(main = c("PCoA biplot", "Response variables projected", 
        "as in PCA with scaling 1"), line = 4) 
    } else title(main = main, line = 4) 
    } 
    invisible() 
} 

biplot(pcoa.ntK, xlabs = 'My x label', ylabs = 'My y label', main = 'My title') 
0

Вы можете проверить исходный код biplot.pcoa, и вы увидите, что это не так уж трудно изменить. Автор пакета решил жестко закодировать метки оси на основе ввода, а также основное название сюжета. Вот модифицированная версия, которая будет сначала проверить, если значения xlab, ylab и main использовались перед использованием предварительно определенных из них:

biplot.pcoa <- function (x, Y = NULL, plot.axes = c(1, 2), dir.axis1 = 1, dir.axis2 = 1, 
      rn = NULL, ...) 
{ 
    k <- ncol(x$vectors) 
    if (k < 2) 
    stop("There is a single eigenvalue. No plot can be produced.") 
    if (k < plot.axes[1]) 
    stop("Axis", plot.axes[1], "does not exist.") 
    if (k < plot.axes[2]) 
    stop("Axis", plot.axes[2], "does not exist.") 
    if (!is.null(rn)) 
    rownames(x$vectors) <- rn 
    args <- list(...) 
    labels = ifelse(c("xlab", "ylab") %in% names(args), c(args$xlab, args$ylab), colnames(x$vectors[, plot.axes])) 
    diag.dir <- diag(c(dir.axis1, dir.axis2)) 
    x$vectors[, plot.axes] <- x$vectors[, plot.axes] %*% diag.dir 
    if (is.null(Y)) { 
    limits <- apply(x$vectors[, plot.axes], 2, range) 
    ran.x <- limits[2, 1] - limits[1, 1] 
    ran.y <- limits[2, 2] - limits[1, 2] 
    xlim <- c((limits[1, 1] - ran.x/10), (limits[2, 1] + 
              ran.x/5)) 
    ylim <- c((limits[1, 2] - ran.y/10), (limits[2, 2] + 
              ran.y/10)) 
    par(mai = c(1, 1, 1, 0.5)) 
    title <- ifelse("main" %in% names(args), args$main, "PCoA ordination") 
    plot(x$vectors[, plot.axes], xlab = labels[1], ylab = labels[2], 
     xlim = xlim, ylim = ylim, asp = 1, 
     main = title) 
    text(x$vectors[, plot.axes], labels = rownames(x$vectors), 
     pos = 4, cex = 1, offset = 0.5) 
    #title(main = "PCoA ordination", line = 2.5) 
    } 
    else { 
    n <- nrow(Y) 
    points.stand <- scale(x$vectors[, plot.axes]) 
    S <- cov(Y, points.stand) 
    U <- S %*% diag((x$values$Eigenvalues[plot.axes]/(n - 
                 1))^(-0.5)) 
    colnames(U) <- colnames(x$vectors[, plot.axes]) 
    par(mai = c(1, 0.5, 1.4, 0)) 
    title <- ifelse("main" %in% names(args), args$main, c("PCoA biplot", "Response variables projected", 
                   "as in PCA with scaling 1")) 
    biplot(x$vectors[, plot.axes], U, xlab = labels[1], ylab = labels[2], main = title) 
    # title(main = c("PCoA biplot", "Response variables projected", 
    #    "as in PCA with scaling 1"), line = 4) 
    } 
    invisible() 
} 

Тогда:

biplot(pcoa.ntK, Y=NULL, plot.axes=c(1,2), rn=ntnames, 
    xlab="PC1 (%)", main = "Main Title") 

enter image description here

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

Смежные вопросы