必要なパッケージ

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

背景

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

主成分分析

主成分分析は、多次元データの変動を、できるだけ情報量を失わずに低次元データに要約する手法です。例えば、ブルーベリーデータに含まれる官能評価結果の要約では、5次元のデータの変動の90%以上を2次元で要約することができます。主成分分析は、多数の変数に含まれる情報を効率的に抽出するのに非常に有効な手法です。

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

まず、Gilbertら(2015）のブルーベリーデータの官能評価の結果について主成分分析を行ってみましょう。

前回と同じで、すべての環境で栽培されていた5品種について抜き出したデータを用います。

bb <- read.csv(here("data", "blueberry_5var.csv"), stringsAsFactors = T) 
dim(bb) # size of data
## [1] 133  19

変数名を確認します。

colnames(bb)  # variable names in bb
##  [1] "gID"            "Year"           "Location"       "Harvest"       
##  [5] "Genotype"       "Panelists"      "Overall.Liking" "Texture"       
##  [9] "Sweetness"      "Sourness"       "Flavor"         "SSC"           
## [13] "TA"             "SSC.TA"         "pH"             "Fructose"      
## [17] "Glucose"        "Sucrose"        "Total.Sugar"

官能評価（sensory evaluation）の結果を抜き出します。7列目から11列目に含まれているので、列番号を指定して抜き出します。

bb.sens <- bb[, 7:11]
colnames(bb.sens)
## [1] "Overall.Liking" "Texture"        "Sweetness"      "Sourness"      
## [5] "Flavor"

3次元の散布図を描くことで、3変数間の関係を把握することができます。また、5変数のpairwise plotを描けば、変数間の関係も把握できます。

pairs(bb.sens, col = bb$Genotype)

折れ線グラフで多次元の変動を可視化することもできます。

matplot(t(bb.sens), type = "l", 
        col = bb$Genotype, lty = 1,
        xlab = "Evaluation element", ylab = "Score")
legend("bottomleft", lty = 1, levels(bb$Genotype), col = 1:nlevels(bb$Genotype))

しかし、こうした図から、5つの変数に含まれる変動を把握することは簡単ではありません。

ここで、変数間の相関を可視化してみます。

r <- cor(bb.sens) # calculate correlation among variables
corrplot::corrplot(r) # use corrplot function in the corrplot package

相関係数のヒートマップを見ると、正の関係も負の関係もありますが、互いに関連をもっており冗長性があることがわかります。

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

res <- prcomp(bb.sens) # PCA

結果を要約します。

summary(res)  # summary of the analysis result
## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5
## Standard deviation     6.4709 5.1059 2.15648 1.28099 0.63042
## Proportion of Variance 0.5611 0.3493 0.06231 0.02199 0.00533
## Cumulative Proportion  0.5611 0.9104 0.97269 0.99467 1.00000

summary関数を使って表示した結果、第1主成分はデータの全変動のうち56％を説明し、第1第2主成分と第3主成分はそれぞれ91％以上、97％以上を説明することがわかります。このような説明量を寄与率（contribution）、累積寄与率（cummulative contribution）とよびます。

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

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 = bb$Genotype) # PC1 vs PC2 
legend("bottomright", legend = levels(bb$Genotype), 
    pch = 1, col = 1:nlevels(bb$Genotype))  # put a legend at the top right
plot(res$x[,2:3], col = bb$Genotype) # PC2 vs PC3

par(op)  #　reset the arrange

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

pairs(res$x[,1:3], col = bb$Genotype)

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

plotly::plot_ly(x = res$x[,1], y = res$x[,2], z = res$x[,3], color = bb$Genotype,
                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の得点が大きいサンプルは、酸っぱさSournessが高く、甘さSweetnessと総合評価Overall.Likingが低いことがわかります。一方、PC1の得点が小さいサンプルでは、その逆です。また、食感Textureと香りFlavorはPC1にはほとんど関与していません。また、PC2の得点が大きいサンプルでは、すべての要素が大きくなることがわかります。

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

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

様々な測定値が混在するデータの主成分分析には気をつけよう

次に、糖度や酸度などの計測値について、主成分分析を行ってみましょう。

colnames(bb)
##  [1] "gID"            "Year"           "Location"       "Harvest"       
##  [5] "Genotype"       "Panelists"      "Overall.Liking" "Texture"       
##  [9] "Sweetness"      "Sourness"       "Flavor"         "SSC"           
## [13] "TA"             "SSC.TA"         "pH"             "Fructose"      
## [17] "Glucose"        "Sucrose"        "Total.Sugar"

このうち、SSC.TAはSSCとTAの比、Total.Sugarは、3つの糖の和なのでこれらは除くことにします。

bb.chem <- bb[, c(12, 13, 15:18)]
dim(bb.chem)
## [1] 133   6

主成分分析を行う （注：この分析は誤りです）。

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

バイプロットを描くと、PC1もPC2も、FructoseとGlucoseが大きな影響をもっていることがわかります。なぜこのような結果が得られたのでしょうか。

まず、分散共分散行列をcorrplot関数で可視化してみます。

cv <- cov(bb.chem)
corrplot::corrplot(cv, is.corr = F)

すると、FructoseとGlucoseの分散が非常に大きいことがわかります。

主成分分析では、大きな分散を持つ変数は結果に大きな影響を与えます。ここでの大きな問題は、対象としている変数が異なる尺度で計測されている点です。例えば、Fructose、Glucose、Sucroseはmg/gFMを単位として測定されています。この単位を変えると、分散が大きくも、小さくもできます。それにより、主成分分析の結果に影響してしまいます。また、pHや糖度、酸度は全く異なる尺度で計測されていて、それをFructoseなどの含量と尺度を合わせることはできません。そのような場合、すべての変数の分散が1になるように標準化し、標準化された値をもとに主成分分析を行います。

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

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

因子負荷量も視覚化します。

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

スケーリングした変数で解析をすると、主成分は様々な変数の変動をうまく捉えられていることがわかります。PC1は、主に糖分に関係し、糖度が高いサンプルは大きな得点をもつ。酸度に関わるpHやTAも関連しているが、SSC、Glucose, Fructoseのように大きな寄与はない。PC2は、pHとTAに関連しており、TAが大きく、pHが低いサンプルが大きな得点をもつ。

クラスタ解析

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

階層的なクラスタリング

ここでは、Gilbertら(2015）のデータのうち、揮発性成分のデータを分析します。

まず、揮発性物質のデータを含めた全データを読みこみます（なお、このデータからは、X564.94.3（(1R)-(−)-Myrtenal）は多型が極めて少なかったため除いてある）。

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

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

vc <- bb[, 19: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: 44

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(1, 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

以下は、上記で用いたサンプル間の距離の定義です。なお、各サンプルは\(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| \]

これまで、サンプル間の距離の定義について述べてきました。階層型クラスタ解析では、互いに近いサンプルは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") # Median 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法を紹介する。

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

kms <- kmeans(data, centers = 5, 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
##   1       0      0        0          9          0
##   2       0      4        0          0          0
##   3       0      0        9          0          0
##   4       0      5        0          0          5
##   5       9      0        0          0          3

EnduraとPrimadonnaの品種は、明確にクラスタリングされていないためか、2つのクラスタに分かれており、その他の品種はそれぞれ1つのクラスタに分類されています。なお、ここで、nstart = 20というオプションは非常に重要です。なぜなら、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))