2013-08-04 20 views
7

私はprcompggbiplotを使って主成分分析をプロットしようとしています。私は単位円の外にデータ値を取得しており、prcompを呼び出す前にデータを再スケールすることができず、データを単位円に拘束することができます。いずれかのコードチャンクのPCAスケーリングggbiplot

require(devtools) 
install_github('ggbiplot','vqv') 

出力を:次のように

wine2=wine[,1:3] 
mean=apply(wine2,2,mean) 
sd=apply(wine2,2,mean) 
for(i in 1:ncol(wine2)){ 
    wine2[,i]=(wine2[,i]-mean[i])/sd[i] 
} 
wine2.pca=prcomp(wine2,scale.=TRUE) 
ggbiplot(wine2.pca,obs.scale=1, 
     var.scale=1,groups=wine.class,ellipse=TRUE,circle=TRUE) 

ggbiplotパッケージがインストール:

data(wine) 
require(ggbiplot) 
wine.pca=prcomp(wine[,1:3],scale.=TRUE) 
ggbiplot(wine.pca,obs.scale = 1, 
     var.scale=1,groups=wine.class,ellipse=TRUE,circle=TRUE) 

私はprcompを呼び出す前に平均値を減算し、標準偏差で割ることによってスケーリングする試み

enter image description here

下記の@Brian Hansonのコメントごとに、私が得ようとしている出力を反映した追加画像を追加しています。

enter image description here

+0

さて、二つのアプローチが同じ答えを与える理由は、デフォルトでprcomp'センターは(prcomp'? '参照)'ということです。あなたの2回目の試みはあなたがする必要のない余分な仕事です。さて、もう一つの問題について:なぜ値は単位円の中にあるべきだと思いますか?そして、あなたは何の価値観を持っていますか?彼らがなぜユニットサークルになるべきなのかを説明することができますか? –

+0

私はその投稿を編集し、予想される結果に匹敵する公開されたプロットを追加します。あなたは何もしていない2番目のチャンクで正しいです。私は変数(矢印)が正しくスケールされていることに満足しています、それは私が単位円の中で拘束しようとしている観測です。 データを中央に置くことに加えてスケーリングした場合、観測値が単位円の外にある可能性があるのはなぜか分かりません。 – scs217

+0

変数をスケーリングすると、それらは分散1と平均0になるように個別にスケールされます。結果の得点が単位円内にある理論的な理由はないと思います。私はちょうどスケーリングとヘルプページの多くの例を見て回って、それらのどれも単位円の中に入るスコアを生成しません。 pcaのことをより深く理解するためには、私はバイプロットから離れていくことをお勧めします。良いプロットの重要な原則の1つを破ります。同じプロットに2つのスケールがありません。 –

答えて

7

私は、プロット関数のコードを編集し、私が望んでいた機能を得ることができました。

ggbiplot2=function(pcobj, choices = 1:2, scale = 1, pc.biplot = TRUE, 
     obs.scale = 1 - scale, var.scale = scale, 
     groups = NULL, ellipse = FALSE, ellipse.prob = 0.68, 
     labels = NULL, labels.size = 3, alpha = 1, 
     var.axes = TRUE, 
     circle = FALSE, circle.prob = 0.69, 
     varname.size = 3, varname.adjust = 1.5, 
     varname.abbrev = FALSE, ...) 
{ 
    library(ggplot2) 
    library(plyr) 
    library(scales) 
    library(grid) 

    stopifnot(length(choices) == 2) 

    # Recover the SVD 
    if(inherits(pcobj, 'prcomp')){ 
    nobs.factor <- sqrt(nrow(pcobj$x) - 1) 
    d <- pcobj$sdev 
    u <- sweep(pcobj$x, 2, 1/(d * nobs.factor), FUN = '*') 
    v <- pcobj$rotation 
    } else if(inherits(pcobj, 'princomp')) { 
    nobs.factor <- sqrt(pcobj$n.obs) 
    d <- pcobj$sdev 
    u <- sweep(pcobj$scores, 2, 1/(d * nobs.factor), FUN = '*') 
    v <- pcobj$loadings 
    } else if(inherits(pcobj, 'PCA')) { 
    nobs.factor <- sqrt(nrow(pcobj$call$X)) 
    d <- unlist(sqrt(pcobj$eig)[1]) 
    u <- sweep(pcobj$ind$coord, 2, 1/(d * nobs.factor), FUN = '*') 
    v <- sweep(pcobj$var$coord,2,sqrt(pcobj$eig[1:ncol(pcobj$var$coord),1]),FUN="/") 
    } else { 
    stop('Expected a object of class prcomp, princomp or PCA') 
    } 

    # Scores 
    df.u <- as.data.frame(sweep(u[,choices], 2, d[choices]^obs.scale, FUN='*')) 

    # Directions 
    v <- sweep(v, 2, d^var.scale, FUN='*') 
    df.v <- as.data.frame(v[, choices]) 

    names(df.u) <- c('xvar', 'yvar') 
    names(df.v) <- names(df.u) 

    if(pc.biplot) { 
    df.u <- df.u * nobs.factor 
    } 

    # Scale the radius of the correlation circle so that it corresponds to 
    # a data ellipse for the standardized PC scores 
    r <- 1 

    # Scale directions 
    v.scale <- rowSums(v^2) 
    df.v <- df.v/sqrt(max(v.scale)) 

    ## Scale Scores 
    r.scale=sqrt(max(df.u[,1]^2+df.u[,2]^2)) 
    df.u=.99*df.u/r.scale 

    # Change the labels for the axes 
    if(obs.scale == 0) { 
    u.axis.labs <- paste('standardized PC', choices, sep='') 
    } else { 
    u.axis.labs <- paste('PC', choices, sep='') 
    } 

    # Append the proportion of explained variance to the axis labels 
    u.axis.labs <- paste(u.axis.labs, 
         sprintf('(%0.1f%% explained var.)', 
           100 * pcobj$sdev[choices]^2/sum(pcobj$sdev^2))) 

    # Score Labels 
    if(!is.null(labels)) { 
    df.u$labels <- labels 
    } 

    # Grouping variable 
    if(!is.null(groups)) { 
    df.u$groups <- groups 
    } 

    # Variable Names 
    if(varname.abbrev) { 
    df.v$varname <- abbreviate(rownames(v)) 
    } else { 
    df.v$varname <- rownames(v) 
    } 

    # Variables for text label placement 
    df.v$angle <- with(df.v, (180/pi) * atan(yvar/xvar)) 
    df.v$hjust = with(df.v, (1 - varname.adjust * sign(xvar))/2) 

    # Base plot 
    g <- ggplot(data = df.u, aes(x = xvar, y = yvar)) + 
    xlab(u.axis.labs[1]) + ylab(u.axis.labs[2]) + coord_equal() 

    if(var.axes) { 
    # Draw circle 
    if(circle) 
    { 
     theta <- c(seq(-pi, pi, length = 50), seq(pi, -pi, length = 50)) 
     circle <- data.frame(xvar = r * cos(theta), yvar = r * sin(theta)) 
     g <- g + geom_path(data = circle, color = muted('white'), 
         size = 1/2, alpha = 1/3) 
    } 

    # Draw directions 
    g <- g + 
     geom_segment(data = df.v, 
        aes(x = 0, y = 0, xend = xvar, yend = yvar), 
        arrow = arrow(length = unit(1/2, 'picas')), 
        color = muted('red')) 
    } 

    # Draw either labels or points 
    if(!is.null(df.u$labels)) { 
    if(!is.null(df.u$groups)) { 
     g <- g + geom_text(aes(label = labels, color = groups), 
         size = labels.size) 
    } else { 
     g <- g + geom_text(aes(label = labels), size = labels.size)  
    } 
    } else { 
    if(!is.null(df.u$groups)) { 
     g <- g + geom_point(aes(color = groups), alpha = alpha) 
    } else { 
     g <- g + geom_point(alpha = alpha)  
    } 
    } 

    # Overlay a concentration ellipse if there are groups 
    if(!is.null(df.u$groups) && ellipse) { 
    theta <- c(seq(-pi, pi, length = 50), seq(pi, -pi, length = 50)) 
    circle <- cbind(cos(theta), sin(theta)) 

    ell <- ddply(df.u, 'groups', function(x) { 
     if(nrow(x) < 2) { 
     return(NULL) 
     } else if(nrow(x) == 2) { 
     sigma <- var(cbind(x$xvar, x$yvar)) 
     } else { 
     sigma <- diag(c(var(x$xvar), var(x$yvar))) 
     } 
     mu <- c(mean(x$xvar), mean(x$yvar)) 
     ed <- sqrt(qchisq(ellipse.prob, df = 2)) 
     data.frame(sweep(circle %*% chol(sigma) * ed, 2, mu, FUN = '+'), 
       groups = x$groups[1]) 
    }) 
    names(ell)[1:2] <- c('xvar', 'yvar') 
    g <- g + geom_path(data = ell, aes(color = groups, group = groups)) 
    } 

    # Label the variable axes 
    if(var.axes) { 
    g <- g + 
     geom_text(data = df.v, 
       aes(label = varname, x = xvar, y = yvar, 
        angle = angle, hjust = hjust), 
       color = 'darkred', size = varname.size) 
    } 
    # Change the name of the legend for groups 
    # if(!is.null(groups)) { 
    # g <- g + scale_color_brewer(name = deparse(substitute(groups)), 
    #        palette = 'Dark2') 
    # } 

    # TODO: Add a second set of axes 

    return(g) 
} 

+1

+1、このfxは今では大幅に改善されていますが、私の質問は、この機能がpsychパッケージの主要なオブジェクトを受け入れることができますか?私はstats :: prcomp()を使って異なる数字になります。 – doctorate