Licenses

CC BY 4.0

多変量解析とは？

多くの変数の関連性を分析し、その特徴を要約・予測するための解析手法のことである．ある特定の解析手法を示すものではなく、下記に示す階層的クラスター解析や主成分分析など多くの解析方法を含んでいる．

パッケージのロード

library関数を用いて、本稿で利用するパッケージをロードしておく．

library(tidyverse)

階層的クラスター解析 (HCA)

クラスター解析は、データ群の中から似たもの同士を集めて、群 (クラスター) に分類する多変量解析の一つの手法である．クラスタリングは、階層的手法と非階層的手法に大別される．階層的手法では、サンプル間の類似度あるいは非類似度 (距離) に基づいて、近接するデータ点同士を逐次的につなぎ、デンドログラム (樹形図) を作成する．非階層的手法では階層構造を形成せずデータを分割する．

ここでは、前者の階層的クラスター解析 (HCA; Hierarchical Cluster Analysis) について、Rで実施する方法について説明する．

Rに標準で格納されているirisデータをデモとして扱う． irisデータとは、3種のアヤメ (setosa、versicolor、virginica) のがくの長さと幅 (Sepal.length、Sepal.Width)、花びらの長さと幅 (Petal.Length、Petal.Width) を50サンプルずつ調査したデータセットである．

irisデータは1列目から4列目が数字データである． 以降の多変量解析のため、数字データだけを抽出して、iris.value変数に入力しておく．

#数値データのみの取得
iris.value <- iris[,1:4]

多変量解析を行う際は、各変数を標準化 (平均0、分散1) しておく必要がある．scale関数で標準化を行っている． また、後のクラスター解析の結果アウトプットである系統樹のラベルのために、列名を変更しておく． ここでは、irisのSpecies列を指定する．

#標準化
iris.value.s<-scale(iris.value)

#列名の変更
rownames(iris.value.s)<-iris$Species

クラスター解析では、各サンプルの多変数の値をもとにサンプル間の距離を計算し、その距離から階層的にクラスターを作成していく． まず、dist関数でサンプル間の距離を計算する．距離の定義にもさまざまあるが、ここでは最も一般的なユークリッド (Euclidean) 距離を用いている．

#dist関数でサンプル間の距離 (ユークリッド距離) を計算する
d <- dist(iris.value.s, method="euclidean")

hclust関数でHCAを実施する．上記で求めた距離行列 (ここでは変数d) を用いる．また、HCAでは、クラスター間の距離をどのように求めるかで様々な手法があるが、ここでは、一般的なウォード (Ward) 法を用いる (method="ward.D2"と設定)．

#hca
hca <- hclust(d, method="ward.D2")

hca変数にHCAの結果が格納された． デンドログラム (樹形図) を描いて、HCAの結果を確認してみる．

# デンドログラムの作成
plot(hca)

大別すると、3つのクラスターに分かれていることがわかる．では、この3つのクラスターにどのように各サンプルが分類されているのかを確認する． cutree関数を用いて、各サンプルが3つのクラスターのどこに分類されているのかを調べる．その結果をclに代入する． クロス集計を行うtable関数を用いて、アヤメの3種類とクラスター分類がどのような対応関係になっているのかを調べる．

#各サンプルがどのクラスターに分類されているのかを調べる
cl <- cutree(hca, k=3)
#table関数でクロス集計
table(cl,iris$Species)
##    
## cl  setosa versicolor virginica
##   1     49          0         0
##   2      1         27         2
##   3      0         23        48

setosa (No.1-50) はクラスター1、versicolor (No.51-100) はクラスター2と3、virginica (No.101-150) はクラスター3に多くが分類されていることがわかる．

したがって、アヤメの種類によって花の形が異なることがわかる．

実習．データフレームdemoを対象に、HCAを行い、系統樹を描出せよ．

ここでは、read.csv関数を用いて、csvファイルを読みこむ． デモデータとして、コマツナを様々な栄養欠乏条件で栽培した際の多元素濃度 (イオノーム) データがまとめられたcsvファイル demo.csv を読み込む．

demo <- read.csv("demo.csv", header=TRUE)

数値データだけを抽出する

#数値データのみの取得

#
#
#

標準化および列名の変更

#標準化

#
#
#

#列名の変更

#
#
#

dist関数でサンプル間の距離を計算する．

#dist関数でサンプル間の距離 (ユークリッド距離) を計算する

#
#
#

hclust関数でHCAを実施する．

#hca
#
#
#

デンドログラム (樹形図) を描出する．

# デンドログラムの作成

#

主成分分析 (PCA)

主成分分析 (PCA; Principal Component Analysis) とは、多次元データのもつ情報をできるだけ損わずに低次元空間に情報を縮約する多変量解析手法の一つである．ビッグデータは、多変量・多次元であるため、その特徴を理解することは難しいですが、PCAによりデータの情報をできる限り損なわずに低次元化することで、そのプロファイル理解を助けることができる．

以下では、RによるPCAの方法について説明する．

上述のクラスター解析と同様に、数値データだけを抽出して、iris.value変数に入力しておく．

iris.value <- iris[,1:4]

prcomp()関数でPCAを行う．多変量解析を行う際は、各変数を標準化(平均0、分散1)しておく必要がある．scale=TRUEで標準化を行っている．

結果がpcaに格納される．

#PCA
pca <- prcomp(iris.value, scale=TRUE)

summary()関数でPCAの結果を確認する．

#summary of PCA results
summary(pca)
## Importance of components:
##                           PC1    PC2     PC3     PC4
## Standard deviation     1.7084 0.9560 0.38309 0.14393
## Proportion of Variance 0.7296 0.2285 0.03669 0.00518
## Cumulative Proportion  0.7296 0.9581 0.99482 1.00000

この結果では、各主成分 (PC) が全変数のもつ変動のうち何%説明できているかを示している．Proportion of Varianceは寄与率を意味し、第一主成分 (PC1)は67.4%、第二主成分 (PC2)は30.3%を説明している．Cumulative Proportionは累積寄与率を意味している．つまり、このPCA結果では、第二主成分までで、全変数の変動の97.6%を説明できていることがわかる．

各サンプルが持つ変動は、主成分スコアとして表現される．主成分スコアは、pca$xに格納されている． head関数で主成分スコアを確認する．

#各サンプルの主成分スコア
head(pca$x)
##            PC1        PC2         PC3          PC4
## [1,] -2.257141 -0.4784238  0.12727962  0.024087508
## [2,] -2.074013  0.6718827  0.23382552  0.102662845
## [3,] -2.356335  0.3407664 -0.04405390  0.028282305
## [4,] -2.291707  0.5953999 -0.09098530 -0.065735340
## [5,] -2.381863 -0.6446757 -0.01568565 -0.035802870
## [6,] -2.068701 -1.4842053 -0.02687825  0.006586116

主成分スコアとirisデータと統合する．

#主成分スコアとirisの情報を統合
pca.data <- cbind.data.frame(pca$x,iris)

#データの中身を確認
head(pca.data, n=5)
##         PC1        PC2         PC3         PC4 Sepal.Length Sepal.Width
## 1 -2.257141 -0.4784238  0.12727962  0.02408751          5.1         3.5
## 2 -2.074013  0.6718827  0.23382552  0.10266284          4.9         3.0
## 3 -2.356335  0.3407664 -0.04405390  0.02828231          4.7         3.2
## 4 -2.291707  0.5953999 -0.09098530 -0.06573534          4.6         3.1
## 5 -2.381863 -0.6446757 -0.01568565 -0.03580287          5.0         3.6
##   Petal.Length Petal.Width Species
## 1          1.4         0.2  setosa
## 2          1.4         0.2  setosa
## 3          1.3         0.2  setosa
## 4          1.5         0.2  setosa
## 5          1.4         0.2  setosa

PC1とPC2をプロットする． ggplot関数で実施する．colorオプションにSpeciesを指定することで、アヤメの種類毎に色分けを行う．

# がくの幅
ggplot(pca.data, aes(x = PC1, y = PC2, color = Species)) + geom_point(size = 3) +
    theme(text = element_text(size = 20))

PC1とPC2の値が近いほど、近くにプロットされる．つまり、アヤメの花の形が類似していることを意味する．irisデータのPCA結果では、versicolorとvirginicaのデータは近くにプロットされており、setosaのデータは少し離れて分布している．setosaと比較して、versicolorとvirginicaは花の形が似ていることが読み取れる．

各主成分が、もとの変数とどのような関係にあるのかを確認してみる．この関係を図示したものを「biplot」と呼ぶ．biplot関数にPCAの結果 (ここではpca変数に格納) を指定することで実行できる．

#biplot
biplot(pca)

この結果では、各変数の矢印の向きと長さが各主成分に対する関連度合いを示している．PC1が大きいサンプルは、花の長さと幅 (petal.Langth, Petal.Width)、がくの長さ (Sepal.Length) の値が大きく、PC2が小さいサンプルは、がくの幅 (Sepal.Width) の値が大きいことを意味している．

この結果は、はじめに生データを可視化したときにも観察されており、このパターンをPCAがうまく捉えていることがわかる．今回デモで扱ったirisデータは次元数がそこまで多くないので、生データから結果は想定できた．が、より多次元・多変量なデータの場合、PCAなどの多変量解析を適用しないと全体の特徴を把握するのは難しい．