主成分凸多角形

Question

とRにおける解析バイプロットIは、次のコードを使用して取り付けられたバイプロットを生成しました：

dd = data.frame(x = runif(10), y=runif(10)) 
pcr = prcomp(~x + y, data=dd) 
biplot(pcr) 

これはX及びYと10個のデータ点の各々に対して軸を示すバイプロットを生成します。

10のデータポイントは、2つの異なるグループ（1つのグループに5つ、他のグループに5つ）で構成されています。 2つのグループの分割を表示するために、各グループの周りに最小凸多角形を持つバイプロットを作成するにはどうすればよいですか？

Answer 1

私は統計::: biplot.defaultと統計::: biplot.prcompに見て、私はあなたが望むものに近いです。このコードは必要に応じて変更できます。 （私は虹彩のデータセットを使用しました）

require(plyr) 

data(iris) 

pcr <- prcomp(~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width, data = iris) 

indiv <- data.frame(pcr$x[,1:2]) 

indiv$species <- iris$Species 

column <- data.frame(pcr$rotation[ ,1:2]) 

n <- nrow(indiv) 

eigenval <- pcr$sdev[1:2] 

eigenval <- eigenval * sqrt(n) 

indiv <- transform(indiv, pc1 = PC1/eigenval[1], pc2 = PC2/eigenval[2]) 

column <- transform(column, pc1 = PC1 * eigenval[1], pc2 = PC2 * eigenval[2]) 

### based on stats:::biplot.default 

unsigned.range <- function(x) c(-abs(min(x, na.rm = TRUE)), abs(max(x, na.rm = TRUE))) 

rangx1 <- unsigned.range(indiv[, 4]) 
rangx2 <- unsigned.range(indiv[, 5]) 
rangy1 <- unsigned.range(column[, 3]) 
rangy2 <- unsigned.range(column[, 4]) 

mylim <- range(rangx1, rangx2) 
ratio <- max(rangy1/rangx1, rangy2/rangx2) 

nspecies <- table(iris$Species) 

# compute the convex hull for each species 
hull <- dlply(indiv[,1:3], .(species), chull) 

# get points connected 
hull <- llply(hull, function(x) c(x, x[1])) 


plot(pc2 ~ pc1, data = indiv, cex = 0.5, col = c("blue", "yellow", "green")[iris$Species], xlim = mylim, ylim = mylim) 

lines(indiv$pc1[hull$setosa], indiv$pc2[hull$setosa] , col = "blue") 

lines(indiv$pc1[cumsum(nspecies)[1] + hull$versicolor], indiv$pc2[cumsum(nspecies)[1] + hull$versicolor], col = "yellow") 

lines(indiv$pc1[cumsum(nspecies)[2] + hull$virginica], indiv$pc2[cumsum(nspecies)[2] + hull$virginica], col = "green") 

par(new = TRUE) 

plot(pc1 ~ pc2, data = column, axes = FALSE, type = "n", xlim = mylim * ratio, ylim = mylim * ratio, xlab = "", ylab = "") 

text(column$pc1, column$pc2, labels = rownames(column), cex = 0.5, col = "red") 

arrows(0, 0, column$pc1 * 0.8, column$pc2 * 0.8, col = "red", length = 0.2) 

axis(3, col = "red") 

axis(4, col = "red")