必要なパッケージ

require(corrplot)
require(factoextra)
require(Rtsne)
require(ape)
require(plotly)

背景

　農学分野では、同じ試料で複数の形質を測定することは珍しいことではない。例えば、圃場試験では、収量を評価することが目的であっても、収量に関連する様々な形質が同時に調べられることが多い。これらの計測された形質を散布図にすることで、何らかの知見を得ることができる場合が多い。しかし、測定される形質の数が多くなると、散布図からデータのばらつきを把握することが困難になる。散布図を用いて直感的に把握できる次元数はせいぜい数次元であり、例えば測定した特徴量が数十個に及ぶと、データのばらつきを把握することは容易でなくなる。本講演では、多次元データのバラツキを要約し、可視化するための手法について説明する。

主成分分析

主成分分析は、多次元データの変動を、できるだけ情報量を失わずに低次元データに要約する手法である。例えば、講義で例示するリモートセンシングデータに含まれる変動の要約では、12次元のデータの変動の90%以上を2次元で要約できる。主成分分析は、多数の変数に含まれる情報を効率的に抽出するのに非常に有効な手法である。

同一種類の複数の変数の解析

まず、Haoら（2015, PLoS ONE 10: e0137748）のリモートセンシングデータをもとに、主成分分析を行ってみよう。この論文では、Landsatなどの衛星から得られるリモートセンシング画像から、各畑で栽培されている作物を予測する機械学習モデルを構築することを目的としている。

リモートセンシングは年間11回実施し、そこから正規化差分植生指数（NDVI）が算出されている。NDVIは以下のように算出される。 \[ NDVI = \frac{\rho(NIR) - \rho(R)}{\rho(NIR) + \rho(R)} \] ここで\(\rho(NIR)\)と\(\rho(R)\)はそれぞれ近赤外線(NIR)強度と赤色(R)強度を表す。

植生は近赤外線を反射し、赤色を吸収するため、植生があるところではNDVIが高くなる。著者らは、栽培作物（綿、ブドウ、トウモロコシ、スイカ、小麦）とその作型（小麦は夏作物との輪作か否かを分ける）により6つに分類された多数の地点のNDVIを測定したデータを使用した。

本研究で使用したデータ（Bole地域の500サイト）のポートが「ndvi.csv」である。 このデータを読み込んでみる。

ndvi <- read.csv("ndvi.csv", row.names = 1, stringsAsFactors = T)
dim(ndvi) # size of data
## [1] 500  12

データの大きさは500×12（500サンプル、12変数）である。

データの最初の6行を確認する。

head(ndvi)  # the first six lines
##     Class  DOY96 DOY112 DOY130 DOY144 DOY159 DOY192 DOY208 DOY228 DOY256 DOY270
## 7  Cotton 0.1394 0.0731 0.0528 0.1085 0.2155 0.7805 0.8382 0.8104 0.6906 0.5407
## 12 Cotton 0.1501 0.0678 0.0561 0.1027 0.1970 0.7080 0.7807 0.7756 0.6272 0.4612
## 23 Cotton 0.1216 0.0680 0.0600 0.1139 0.1794 0.7612 0.8260 0.7988 0.6538 0.5415
## 30 Cotton 0.1304 0.0654 0.0716 0.1138 0.1592 0.7405 0.7987 0.7084 0.5952 0.4576
## 40 Cotton 0.1113 0.0589 0.0497 0.0868 0.1539 0.7246 0.8225 0.8243 0.6544 0.5754
## 44 Cotton 0.1389 0.0727 0.0626 0.0868 0.1526 0.7743 0.8165 0.8128 0.6632 0.5353
##    DOY289
## 7  0.3406
## 12 0.3439
## 23 0.3392
## 30 0.2713
## 40 0.3370
## 44 0.3729

折れ線グラフで変動を可視化する。

matplot(t(ndvi[, -1]), type = "l", col = ndvi$Class, lty = 1,
        xlab = "Month", ylab = "NDVI")
legend("topleft", lty = 1, levels(ndvi$Class), col = 1:nlevels(ndvi$Class))

変数間の相関を可視化する。

r <- cor(ndvi[, -1]) # calculate correlation among variables (without Class)
corrplot::corrplot(r) # use corrplot function in the corrplot package

ヒートマップを見ると、正の関係も負の関係もあるが、変数間の関係は強く、冗長性が高いことがわかる。

主成分分析を行ってみる。主成分分析には、Rに標準で実装されているprcomp関数を使用する。

res <- prcomp(ndvi[, -1]) # PCA without Class variable

結果を要約する。

summary(res)  # summary of the analysis result
## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     0.6399 0.2476 0.15680 0.08867 0.06230 0.05626 0.04036
## Proportion of Variance 0.7943 0.1189 0.04769 0.01525 0.00753 0.00614 0.00316
## Cumulative Proportion  0.7943 0.9132 0.96094 0.97619 0.98372 0.98986 0.99302
##                            PC8     PC9    PC10    PC11
## Standard deviation     0.03644 0.03404 0.02632 0.02054
## Proportion of Variance 0.00258 0.00225 0.00134 0.00082
## Cumulative Proportion  0.99559 0.99784 0.99918 1.00000

summary関数を使って表示した結果、第1主成分はデータの全変動のうち79％を説明でき、第1第2主成分と第3主成分はそれぞれ91％と96％を説明できることがわかる。

主成分の寄与率を図示するには以下のようにする。

plot(res) # barplot of eigenvalues

図示されるのは固有値であり、この図では各主成分の寄与がわかりにくい。そこで、factoextraパッケージの関数を用いて図示する。

fviz_screeplot(res, addlabels = TRUE)

次に、主成分得点の散布図を描いてみる。まず、第1主成分と第2主成分で散布図を作成する。作物や作物種の区分による違いを見るために、それらで色分けしたグラフを描く。

op <- par(mfrow = c(1,2)) #　arrange graphs in 1 x 2
plot(res$x[,1:2], col = ndvi$Class) # PC1 vs PC2 
plot(res$x[,2:3], col = ndvi$Class) # PC2 vs PC3
legend("topright", legend = levels(ndvi$Class), 
    pch = 1, col = 1:nlevels(ndvi$Class))  # put a legend at the top right

par(op)  #　reset the arrange

plotlyで2Dプロットを描画する。

plotly::plot_ly(x = res$x[,1], y = res$x[,2], 
                type = "scatter", mode = "markers", color = ndvi$Class)

PC1〜3を同時に確認するには、pairs関数を使えば良い。

pairs(res$x[,1:3], col = ndvi$Class)

あるいは、plotlyで3Dプロットを描画することもできる。

plotly::plot_ly(x = res$x[,1], y = res$x[,2], z = res$x[,3], color = ndvi$Class,
                type = "scatter3d", mode = "markers")

主成分スコアの散布図から、冬作物のコムギは夏作物とは異なる位置にあることがわかる。主成分得点の散布図は、冬作コムギが夏作コムギと異なる位置にあり、同じカテゴリーに属する夏作のサンプルはいくつかのクラスターを形成するが、それらは互いに重なり合い、冬作-夏作間のように違いが明確でないことを示している。

これら最初の3つの主成分は、元の変数のどのような変動を捉えているのか。固有ベクトルを確認することで把握できるが、このデータのように変数が多い場合は把握することは難しい。バイプロットグラフで可視化して見る。

op <- par(mfrow = c(1,2))   
biplot(res, choices = 1:2)  #　PC1 vs. PC2
biplot(res, choices = 2:3)  #　PC2 vs. PC3

par(op)

変数の名前がついた矢印の大きさと向きは、各主成分とその変数との関係の強さを表す。例えば、左図では、PC1（夏作物）の値が大きいサンプルは、夏期（DOY192、DOY208、DOY228）にNDVIが高く、春期（DOY130、144、159）にNDVIが低くなることがわかる。一方、PC1の値が小さいサンプル（冬作物）では、前者が低く、後者が高い傾向がある。PC2の値が大きいサンプル（夏作物、輪作コムギ）では、秋にNDVI（DOY256、DOY270、DOY289）が高くなり、値が小さいサンプル（スイカ）ではその逆である。

なお、因子負荷量（変数と主成分間の相関）を図示して、主成分の意味づけをする場合もある。 ここでは、factoextraパッケージの関数を用いて図示する。各変数の主成分への寄与をもとに色付けしてある。

fviz_pca_var(res, col.var="contrib", repel = T)

なお、最近、多次元データの視覚化法として、t-SNEという方法がよく用いられるようになっている。ここでは、パッケージRtsneの関数を用いて、解析をしてみよう。

res <- Rtsne(ndvi, perplexity = 30)
plot(res$Y, col = ndvi$Class)

主成分分析では、互いに重なってしまっていた作物が、離れて散布されていることがわかる。このよううに、t-SNE では近い関係がより近く、遠い関係がより遠く表されるため、よく似たデータのまとまりが明瞭になる場合がある。なお、t-SNE はハイパーパラメータperplexity に結果が強く依存する。

res <- Rtsne(ndvi, perplexity = 5)
plot(res$Y, col = ndvi$Class, asp = 1)

3 次元の布置を計算して、3 次元で確認してみる。

res <- Rtsne(ndvi, dims = 3, perplexity = 10)
plot_ly(x = res$Y[,1], y = res$Y[,2], z = res$Y[,3], color = ndvi$Class,
                type = "scatter3d", mode = "markers")

様々な測定値が混在するデータの解析

次に分析するデータは、先に分析したデータとは異なり、異なる種類のデータが混在している。ここで例として使用するのは、Tadesse et al. (2015, PLoS ONE 10: e0141339) で使用されたデータである。このデータは、2011年から2012年にかけて、乾燥地農業研究国際センター（ICARDA）のテルハディア（シリア）パイロットサイトで天水栽培と灌漑栽培の条件で収集されたコムギのデータで、表現型データは収量、品質、その他の農学形質について測定されたものである。この論文では、この表現型データを用いて、ゲノムデータとの関連性解析が行われた。ここでは、このデータをもとに主成分分析を行う。

データをロードしてみる。

wp <- read.csv("wheat.csv", row.names = 1) # read data
dim(wp)
## [1] 118  19

データを確認し、データのサブセット（最初の6変数）を抽出する。

head(wp)  # display first 6 rows
##    DH_rainfed DH_irrigated PlantHeight_rainfed PlanHeight_irrigated
## G1      149.0        154.5                67.5                 77.5
## G2      148.0        155.5                67.5                 77.5
## G3      148.0        157.0                70.0                 85.0
## G4      147.0        155.5                70.0                 87.5
## G5      147.5        156.0                75.0                 92.5
## G6      144.0        152.0                72.5                 92.5
##    Yield_rainfed Yield_irrigated Yield_average  TKW   TW  PSI Protein  FAB FDT
## G1      3691.923        5679.167      4685.545 29.0 74.2 45.1    15.7 63.0 4.0
## G2      3080.629        5081.250      4080.939 29.1 73.6 38.4    15.5 63.0 5.7
## G3      3150.219        6241.667      4695.943 29.1 73.2 42.9    15.7 63.0 4.0
## G4      2762.610        5132.292      3947.451 38.0 77.2 42.7    14.6 63.0 4.5
## G5      3368.640        5612.500      4490.570 33.3 75.2 47.1    14.9 60.0 5.2
## G6      3068.092        6566.667      4817.379 33.4 76.2 42.9    15.3 65.5 3.2
##    FST MTI  P   L P.L   W
## G1 9.0  40 57 147 0.4 253
## G2 7.3  37 56 139 0.4 247
## G3 7.5  46 54 170 0.3 254
## G4 7.6  43 47 158 0.3 192
## G5 8.0  44 48 161 0.3 233
## G6 4.2  50 64 129 0.5 207
data <- wp[, 1:6]   # analyze the first 6 variables

主成分分析を行う （注：この分析は不正確です）。

res <- prcomp(data)   # Principal component analysis. However, this is wrong!
biplot(res) # Display biplot

バイプロットを描くと、PC1が灌漑栽培の収量、PC2が天水栽培の収量に関係し、その他の変数はほとんど寄与していないことがわかる。

なぜこのような結果が得られたのだろうか。

まず、分散共分散行列を見てみる。

cv <- cov(data)
cv
##                        DH_rainfed DH_irrigated PlantHeight_rainfed
## DH_rainfed             15.5755650    14.251159           -5.479502
## DH_irrigated           14.2511589    15.038099           -2.676010
## PlantHeight_rainfed    -5.4795017    -2.676010           42.749529
## PlanHeight_irrigated   -0.6065298     1.272273           31.422932
## Yield_rainfed        -625.3846615  -591.914271            6.185322
## Yield_irrigated       304.9335936   399.925229         -575.042227
##                      PlanHeight_irrigated Yield_rainfed Yield_irrigated
## DH_rainfed                     -0.6065298   -625.384661        304.9336
## DH_irrigated                    1.2722729   -591.914271        399.9252
## PlantHeight_rainfed            31.4229321      6.185322       -575.0422
## PlanHeight_irrigated           53.0426083   -597.627073        468.3545
## Yield_rainfed                -597.6270725 127338.066374      24464.3820
## Yield_irrigated               468.3545146  24464.381974     431670.1230

共分散行列をcorrplot関数で可視化する。

corrplot::corrplot(cv, is.corr = F)

すると、収量の分散が非常に大きく、他の値を支配していることがわかる。

主成分分析では、大きな分散を持つ変数は結果に大きな影響を与える。ここでの大きな問題は、見出しまでの日数（days）、植物の高さ（cm）、収量（kg/ha）がすべて異なるスケールである点である。例えば、植物の高さをmm単位で測定した場合、高さの分散が大きくなり、主成分分析の結果に影響する。また、収量をt/haで測定すれば、収量の分散を小さくすることができる。しかし，長さと重さのスケールを合わせる標準はない．そのような場合，すべての変数の分散が1になるように標準化しなければならない（すなわち，「スケールされた」変数による分析）．

これを行うには、prcomp関数の引数にscale = Tを設定する。

res <- prcomp(data, scale = T) # Principal component analysis with standardized data
biplot(res)

スケーリングした変数で解析をすると、主成分は様々な変数の変動をうまく捉えられている。PC1は、主に出穂日数と天水栽培の収量（以下、天水収量）に関係し、出穂日が大きく（遅く）、天水収量が少ない系統は値が大きく、逆の傾向の系統は値が小さい。PC2は主に草丈に関係し、値が大きいほど草丈が低く、値が小さいほど草丈が高いことを示す。例えば、G55系統は、出穂までの日数が長く、草丈が高く、天水収量が低い。一方、灌漑収量は、第1成分、第2成分との関連性が低い。

クラスター解析

多くの対象を、その多次元的な特徴に基づいてグループ（クラスター）に分類することが有効な場合がある。例えば、遺伝資源に含まれる品種をDNA多型データに基づいてグループ化できれば、グループの情報をもとに遺伝資源の変動を体系化することができる。クラスター解析は、多数のデータを少数のグループにまとめることで、データの変動を要約しようとする方法である。本講義では、多数のデータを階層的なグループに分類する「階層的クラスター分析」と、決められた数のデータを決められた数のグループに分類する「非階層的クラスター分析」の2種類のクラスター分析について概説する。

階層的なクラスタリング

ここでは、Gilbertら（2015, PLoS ONE 10: e0138494）のブルーベリーデータを使用する。この研究は、複数のブルーベリー品種を用いた多環境試験（複数の試験場、複数の栽培年による栽培試験）において、嗜好性などの官能評価結果と、酸度、糖度、多数の揮発成分などの生化学的特性との関係を分析したものである。ここでは、揮発性成分に関するデータを分析する。

まず、ブルーベリーのデータを読みこむ。

bb <- read.csv("berry.csv", stringsAsFactor = T,
               row.names = 1) # load data from the csv file
#head(bb)  # first six rows of data

次に、読み込んだデータから揮発性成分のデータ（17列目から最後の列まで）を抽出する。なお、データ量が多いため、2013年のデータのみを分析する。また、品種名（Genotype）を抽出する。

vc <- bb[, 17:ncol(bb)]   # extract the volatile data
data <- vc[bb$Year == 2013, ]       # extract data in 2013
variety <- bb$Genotype[bb$Year == 2013]   # extract variety name

次に、データを平均0、分散1にスケーリングしてく。

data <- scale(data)   # standardizing
d <- dist(data)  # calculate Euclid distance

hclust関数でクラスター分析を行う。

tre <- hclust(d)   #  hierarchical clustering
tre      # show the result
## 
## Call:
## hclust(d = d)
## 
## Cluster method   : complete 
## Distance         : euclidean 
## Number of objects: 50

Callには、実行されたコマンドが表示される。“Cluster method”はクラスター分析の方法（クラスター間の距離の定義）、“Distance”は距離の計算方法を表しています。“Number of Objects”は分類されたオブジェクト（この場合はサンプル）の数です。

クラスター分析の結果をデンドログラム（樹形図）で表示してみる。

plot(tre)     # draw a dendrogram

上の図は、関数 hclust で直接得られた樹形図である。パッケージ ape を用いると、樹形図を様々な表現方法で描くことができる。そのためには、関数 hclust で得られた結果を、パッケージ ape で定義されている “phylo” というクラスに変換する必要がある。

phy <- ape::as.phylo(tre)  # convert the result of hclust to phylo object

パッケージapeの “phylo”クラスは、様々な表現方法で樹形図を描画できる。いろいろな樹形図を試してみる。

par(mfrow = c(2, 2), cex = 0.6) # draw figures in a 2x2 format and with a small font size (0.6)
# draw dendrograms in four different representation styles
plot(phy, type = "phylogram", tip.color = as.numeric(variety)) 
plot(phy, type = "cladogram", tip.color = as.numeric(variety))
plot(phy, type = "fan", tip.color = as.numeric(variety))
plot(phy, type = "unrooted", tip.color = as.numeric(variety))

クラスター分析は、サンプルとクラスター間の距離を計算し、計算された距離に基づいてクラスタリングを行うものである。そのため、距離の定義が異なれば、得られる結果も異なる。ここでは、サンプル間とクラスタ間の距離の定義について説明する。

サンプル間の距離を計算するためには、様々な定義がある。異なる定義の距離に基づいて、樹形図を描いてみる。

# try different methods for calculating distance
par(mfrow = c(2, 2), cex = 0.6) # draw figures in a 2x2 format and with a small font size (0.6)
d <- dist(data, method = "euclidean") # Euclid distance
plot(ape::as.phylo(hclust(d)), tip.color = as.numeric(variety)) # draw the dendrogram
d <- dist(data, method = "manhattan") # Manhattan distance
plot(ape::as.phylo(hclust(d)), tip.color = as.numeric(variety)) # draw the dendrogram
d <- dist(data, method = "minkowski", p = 1.5) # Minkowski distance
plot(ape::as.phylo(hclust(d)), tip.color = as.numeric(variety)) # draw the dendrogram
d <- as.dist(1 - cor(t(data))) # Correlation based distance
plot(ape::as.phylo(hclust(d)), tip.color = as.numeric(variety)) # draw the dendrogram

以下は、上記で用いたサンプル間の距離の定義である。なお、各サンプルは\(q\)個の特徴量で記述され、\(i\)番目のサンプルのデータベクトルを\(\mathbf{x}_i=(x_{i1},...,x_{iq} )^T\)、\(j\)番目のサンプルのデータベクトルを\(\mathbf{x}_j=(x_{j1}/...、x_{jq} )^T\)で示すとすると このとき、サンプル\(i\)と\(j\)の距離\(d(\mathbf{x}_i,\mathbf{x}_j)\) は次のように定義される。

  • ユークリッド距離 \[ d(\mathbf{x}_i,\mathbf{x}_j) = \sqrt{\sum_{k=1}^q (x_{ik} - x_{jk})^2} \]
  • マハラノビス（Manhattan）距離 \[ d(\mathbf{x}_i,\mathbf{x}_j) = \sum_{k=1}^q \left|x_{ik} - x_{jk}\right| \]
  • ミンコフスキー（Minkowski）距離 \[ d(\mathbf{x}_i,\mathbf{x}_j) = \left(\sum_{k=1}^q \left|x_{ik} - x_{jk}\right|^p\right)^{1/p} \]
  • 相関に基づく距離 \[ d(\mathbf{x}_i,\mathbf{x}_j) = r_{ij} = 1 - \frac{\sum_{k=1}^q (x_{ik} - \bar x_i)(x_{jk} - \bar x_j)} {\sqrt{\sum_{k=1}^q (x_{ik} - \bar x_i)^2}\sqrt{\sum_{k=1}^q (x_{jk} - \bar x_j)^2}} \] ここで、\(\bar x_i = \frac{1}{n} \sum_{k=1}^q x_{ik}\)

これまで、サンプル間の距離の定義について述べてきた。階層型クラスター分析では、互いに近いサンプルは1つのクラスターにまとめられ、サンプルとクラスター、または、複数のクラスターが、さらに上位のクラスターにまとめられる。そのため、サンプル間の距離だけでなく、サンプルとクラスタ、あるいはクラスタ間の距離も定義する必要がある。 ま ず、クラスタ間距離の様々な定義に基づいて樹形図を描いてみよう。hclust 関数では、クラスタ間距離の計算方法（定義）をオプションの method で指定することができる。

# Clustering analysis with different methods (definitions of between-cluster distances)
d <- dist(data) # calculate Euclid distance
par(mfrow = c(3, 2), cex = 0.6) # draw graph in a 3x2 format with a small font size (0.6)
tre <- hclust(d, method = "complete") # Complete linkage
plot(as.phylo(tre), tip.color = as.numeric(variety))
tre <- hclust(d, method = "single") # Single linkage
plot(as.phylo(tre), tip.color = as.numeric(variety))
tre <- hclust(d, method = "average") # average linkage
plot(as.phylo(tre), tip.color = as.numeric(variety))
tre <- hclust(d, method = "centroid") # Centroid method
plot(as.phylo(tre), tip.color = as.numeric(variety))
tre <- hclust(d, method = "median") # Cmedian method
plot(as.phylo(tre), tip.color = as.numeric(variety))
tre <- hclust(d, method = "ward.D2") # Ward's method
plot(as.phylo(tre), tip.color = as.numeric(variety))

関数hclustで指定できるクラスタ間距離の定義を示す。サンプル間の距離\(d( \mathbf{x}_i,\mathbf{x}_j)\) をもとに、クラスタAとBの間の距離\(d_{AB}\)を以下のように計算する。

  • 最大距離法（完全連結法）（method=“complete”とする） \[ d_{AB}=\underset{i\in A,j\in B}{\max} d(\mathbf{x}_i,\mathbf{x}_j) \]
  • 最短距離法（単結合法）（method=“single”とする） \[ d_{AB}=\underset{i\in A,j\in B}{\min} d(\mathbf{x}_i,\mathbf{x}_j) \]
  • 平均距離法（method=“average”とする) \[ d_{AB}=\frac{1}{n_A n_B} \sum_{i\in A}\sum_{j\in B}d(\mathbf{x}_i,\mathbf{x}_j) \]

以下の3つの定義では、クラスタAとBが合併して新しいクラスタCを形成するとき、新しいクラスタCとA、B以外のクラスタOとの距離\(d_{CO}\)は次のように定義される。クラスタAとBの距離を\(d_{AB}\)、クラスタAとOの距離を\(d_{AO}\)、クラスタBとOの距離を\(d_{BO}\)、クラスタA、B、Oに含まれるサンプル数を\(n_A\), \(n_B\), \(n_O\)と表記する。

  • 重心法 （method=“centroid”とする） \[ d_{CO}^2=\frac{n_A}{n_A+n_B} d_{AO}^2+\frac{n_B}{n_A+n_B} d_{BO}^2-\frac{n_A n_B}{(n_A+n_B)^2} d_{AB}^2 \]
  • メディアン法（method=“median”とする） \[ d_{CO}=\frac{1}{2} d_{AO}+\frac{1}{2} d_{BO}-\frac{1}{4} d_{AB} \]
  • ウォード法（method=“ward.D2”とする） \[ d_{CO}^2=\frac{n_A+n_O}{n_A+n_B+n_O} d_{AO}^2 +\frac{n_B+n_O}{n_A+n_B+n_O} d_{BO}^2-\frac{n_O}{n_A+n_B+n_O} d_{AB}^2 \]

非階層的クラスタリング

ある程度の数のグループに分類する場合、階層的に分類する必要はない。ここでは、非階層的なクラスタリング解析手法の一つであるk-means法を紹介する。

まず、サンプル間の距離を計算する。

d <- dist(data) # Euclid distance

k-means関数を使って5つのグループに分類してみる。サンプルを6（品種の数）にクラスタリングする。

kms <- kmeans(data, centers = 6, nstart = 20) # Repeat the analysis with 20 different starting points
table(kms$cluster, variety) # compare the clustering result with the varietal classes
##    variety
##     Emerald Endura Farthing Meadowlark Primadonna Scintilla
##   1       0      0        0          0          0         6
##   2       0      4        0          0          0         0
##   3       0      0        9          0          0         0
##   4       9      0        0          0          3         0
##   5       0      0        0          9          0         0
##   6       0      5        0          0          5         0

EnduraとPrimadonnaの品種は、明確にクラスタリングされていないためか、2つのクラスタに分かれており、その他の品種はそれぞれ1つのクラスタに分類されている。

最後に、クラスタリング結果をPCスコアのプロットで可視化する。

pca <- prcomp(data) # principal component analysis
par(mfrow = c(1,2)) # 1x2 format
plot(pca$x[,1:2], pch = kms$cluster, col = as.numeric(variety))
legend("topright", levels(variety), col = 1:nlevels(variety), pch = '-')
plot(pca$x[,3:4], pch = kms$cluster, col = as.numeric(variety))

3D散布図により変動を可視化する。

plotly::plot_ly(x = pca$x[,1], y = pca$x[,2], z = pca$x[,3], color = variety,
                type = "scatter3d", mode = "markers")