#必要なパッケージ
require(tidyverse)
## 要求されたパッケージ tidyverse をロード中です
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6 ✔ purrr 0.3.4
## ✔ tibble 3.1.8 ✔ dplyr 1.0.9
## ✔ tidyr 1.2.0 ✔ stringr 1.4.0
## ✔ readr 2.1.2 ✔ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
require(dendextend)
## 要求されたパッケージ dendextend をロード中です
##
## ---------------------
## Welcome to dendextend version 1.16.0
## Type citation('dendextend') for how to cite the package.
##
## Type browseVignettes(package = 'dendextend') for the package vignette.
## The github page is: https://github.com/talgalili/dendextend/
##
## Suggestions and bug-reports can be submitted at: https://github.com/talgalili/dendextend/issues
## You may ask questions at stackoverflow, use the r and dendextend tags:
## https://stackoverflow.com/questions/tagged/dendextend
##
## To suppress this message use: suppressPackageStartupMessages(library(dendextend))
## ---------------------
##
##
## 次のパッケージを付け加えます: 'dendextend'
##
## 以下のオブジェクトは 'package:stats' からマスクされています:
##
## cutree
require(kohonen)
## 要求されたパッケージ kohonen をロード中です
##
## 次のパッケージを付け加えます: 'kohonen'
##
## 以下のオブジェクトは 'package:purrr' からマスクされています:
##
## map
require(Rtsne)
## 要求されたパッケージ Rtsne をロード中です
require(kernlab)
## 要求されたパッケージ kernlab をロード中です
##
## 次のパッケージを付け加えます: 'kernlab'
##
## 以下のオブジェクトは 'package:purrr' からマスクされています:
##
## cross
##
## 以下のオブジェクトは 'package:ggplot2' からマスクされています:
##
## alpha
ここでは、Gilbertら(2015, PLoS ONE 10: e0138494)のブルーベリーデータを解析対象データとして使用する。この研究は、複数のブルーベリー品種を用いた多環境(複数の試験場、複数の栽培年)試験(Multi-environmental trials: MET)において、嗜好性などの官能評価結果と、酸度、糖度、多数の揮発成分などの生化学的特性との関係を分析したものである。ここでは、これらの変数のうち、揮発性成分のデータについて解析を行う。
データを読み込む。
data <- read.csv("blueberry.csv", row.names = 1, stringsAsFactors = T)
最初の6行を表示する。
head(data)
## Year Location Harvest. Panelists Genotype SSC TA pH
## 2012 C1 Emerald 2012 C 1 95 Emerald 10.4 0.5110 3.43
## 2012 C1 Endura 2012 C 1 95 Endura 10.5 0.6199 3.36
## 2012 C1 Farthing 2012 C 1 95 Farthing 10.6 0.3064 3.55
## 2012 C1 Meadowlark 2012 C 1 95 Meadowlark 9.2 0.2631 3.78
## 2012 C1 Primadonna 2012 C 1 95 Primadonna 11.3 0.2346 3.95
## 2012 C1 Scintilla 2012 C 1 95 Scintilla 12.9 0.4719 3.40
## Overall.Liking Texture Sweetness Sourness Flavor Fructose
## 2012 C1 Emerald 23.5 24.9 22.4 10.7 25.0 48.74
## 2012 C1 Endura 17.5 25.1 19.2 22.6 27.8 47.79
## 2012 C1 Farthing 24.2 23.9 23.6 11.8 26.5 51.31
## 2012 C1 Meadowlark 22.3 28.7 21.5 11.2 24.4 42.57
## 2012 C1 Primadonna 24.0 23.6 25.5 10.4 25.8 56.83
## 2012 C1 Scintilla 24.3 23.6 24.6 14.4 24.9 62.40
## Glucose Sucrose X590.86.3 X616.25.1 X107.87.9 X110.62.3
## 2012 C1 Emerald 49.39 0.13 1.05 7.20 3.79 4.79
## 2012 C1 Endura 46.49 0.19 1.40 8.31 3.56 4.32
## 2012 C1 Farthing 51.66 0.18 0.82 5.27 21.69 4.53
## 2012 C1 Meadowlark 45.11 0.12 0.85 5.25 5.87 4.20
## 2012 C1 Primadonna 56.35 0.21 1.82 7.33 6.88 7.47
## 2012 C1 Scintilla 62.46 1.76 1.76 5.28 0.77 4.83
## X105.37.3 X623.42.7 X123.51.3 X1576.87.0 X71.41.0 X1576.95.0
## 2012 C1 Emerald 0 0.11 0.33 4.46 0.94 1.92
## 2012 C1 Endura 0 0.43 0.55 3.73 0.70 2.32
## 2012 C1 Farthing 0 0.13 0.59 3.05 0.87 1.21
## 2012 C1 Meadowlark 0 0.20 0.24 2.12 1.23 1.52
## 2012 C1 Primadonna 0 0.36 0.50 4.93 1.53 1.67
## 2012 C1 Scintilla 0 0.04 0.99 2.30 0.81 1.36
## X108.88.3 X556.24.1 X66.25.1 X6728.26.3 X928.95.0 X111.27.3
## 2012 C1 Emerald 0.03 0.09 1196.56 4345.51 355.57 190.77
## 2012 C1 Endura 0.08 8.70 907.13 5165.84 401.18 254.44
## 2012 C1 Farthing 0.00 7.65 1058.29 3146.05 358.01 166.66
## 2012 C1 Meadowlark 0.02 43.87 1380.10 4123.45 390.25 245.92
## 2012 C1 Primadonna 0.16 216.42 1766.28 4999.20 290.46 226.99
## 2012 C1 Scintilla 0.03 6.67 1379.83 3488.36 272.86 165.15
## X123.92.2 X110.43.0 X111.71.7 X1191.16.8 X106.70.7
## 2012 C1 Emerald 0 0.72 0.90 0.00 0.08
## 2012 C1 Endura 0 1.35 0.53 0.22 0.26
## 2012 C1 Farthing 0 1.61 0.35 0.28 0.05
## 2012 C1 Meadowlark 0 2.62 0.49 0.00 0.09
## 2012 C1 Primadonna 0 0.78 1.02 0.00 0.25
## 2012 C1 Scintilla 0 1.69 0.78 0.21 0.21
## X7785.70.8 X928.68.7 X142.62.1 X110.93.0 X111.13.7 X123.66.0
## 2012 C1 Emerald 0.00 0.32 0.42 2.42 0.01 1.54
## 2012 C1 Endura 0.05 0.43 0.33 2.54 0.16 1.42
## 2012 C1 Farthing 0.00 0.07 0.37 1.53 0.00 0.04
## 2012 C1 Meadowlark 0.00 0.46 0.33 12.18 0.15 3.00
## 2012 C1 Primadonna 0.00 0.31 0.65 5.21 0.12 1.96
## 2012 C1 Scintilla 0.00 0.10 0.50 1.31 0.00 0.00
## X3681.71.8 X142.92.7 X2497.18.9 X104.76.7 X5989.27.5
## 2012 C1 Emerald 0 4.27 0.00 0 0.79
## 2012 C1 Endura 0 10.75 20.05 0 4.94
## 2012 C1 Farthing 0 1.92 5.78 0 0.43
## 2012 C1 Meadowlark 0 6.63 0.00 0 4.22
## 2012 C1 Primadonna 0 7.46 0.00 0 2.51
## 2012 C1 Scintilla 0 2.10 5.75 0 0.00
## X470.82.6 X122.78.1 X821.55.6 X78.70.6 X124.19.6 X2639.63.6
## 2012 C1 Emerald 1.87 0.13 0.50 36.42 0.74 0.21
## 2012 C1 Endura 10.63 0.42 4.79 16.18 0.69 0.18
## 2012 C1 Farthing 3.29 0.07 0.32 10.86 0.26 0.29
## 2012 C1 Meadowlark 1.26 0.20 1.62 21.46 0.25 0.18
## 2012 C1 Primadonna 1.92 0.07 0.34 10.41 0.89 0.20
## 2012 C1 Scintilla 0.21 0.00 2.25 2.90 0.72 0.00
## X53398.83.7 X112.40.3 X119.36.8 X106.26.3 X141.27.5
## 2012 C1 Emerald 0.71 0.42 1.42 0.70 1.63
## 2012 C1 Endura 3.38 0.05 0.38 0.61 1.94
## 2012 C1 Farthing 1.41 0.01 0.11 0.39 0.52
## 2012 C1 Meadowlark 0.46 0.25 0.00 0.18 0.03
## 2012 C1 Primadonna 0.69 0.17 0.00 0.55 1.67
## 2012 C1 Scintilla 0.15 0.00 0.21 0.26 0.87
## X112.12.9 X629.50.5 X3879.26.3 X689.67.8 X1139.30.6
## 2012 C1 Emerald 0.38 0.11 0.00 2.32 1.78
## 2012 C1 Endura 1.39 0.11 0.14 3.30 5.46
## 2012 C1 Farthing 0.11 0.05 0.00 0.66 2.08
## 2012 C1 Meadowlark 0.85 0.09 0.42 14.76 4.74
## 2012 C1 Primadonna 0.05 0.06 0.00 1.52 1.65
## 2012 C1 Scintilla 1.67 0.06 0.07 2.14 4.64
## X5989.33.3 X43219.68.7 X582.16.1
## 2012 C1 Emerald 0.14 1.27 0.00
## 2012 C1 Endura 1.65 0.66 0.31
## 2012 C1 Farthing 0.15 0.20 0.22
## 2012 C1 Meadowlark 0.81 1.50 0.12
## 2012 C1 Primadonna 0.83 0.41 0.17
## 2012 C1 Scintilla 0.00 0.00 0.00
年(Year)のデータを確認する。全部で3年分のデータがある。
table(data$Year)
##
## 2012 2013 2014
## 51 50 50
2012年、13年のデータのみを解析する(2014年のデータは、後で分類精度を評価するのに用いる)。
met1213<- data %>% filter(Year == 2012 | Year == 2013) %>% select(starts_with("X"))
head(met1213)
## X590.86.3 X616.25.1 X107.87.9 X110.62.3 X105.37.3 X623.42.7
## 2012 C1 Emerald 1.05 7.20 3.79 4.79 0 0.11
## 2012 C1 Endura 1.40 8.31 3.56 4.32 0 0.43
## 2012 C1 Farthing 0.82 5.27 21.69 4.53 0 0.13
## 2012 C1 Meadowlark 0.85 5.25 5.87 4.20 0 0.20
## 2012 C1 Primadonna 1.82 7.33 6.88 7.47 0 0.36
## 2012 C1 Scintilla 1.76 5.28 0.77 4.83 0 0.04
## X123.51.3 X1576.87.0 X71.41.0 X1576.95.0 X108.88.3 X556.24.1
## 2012 C1 Emerald 0.33 4.46 0.94 1.92 0.03 0.09
## 2012 C1 Endura 0.55 3.73 0.70 2.32 0.08 8.70
## 2012 C1 Farthing 0.59 3.05 0.87 1.21 0.00 7.65
## 2012 C1 Meadowlark 0.24 2.12 1.23 1.52 0.02 43.87
## 2012 C1 Primadonna 0.50 4.93 1.53 1.67 0.16 216.42
## 2012 C1 Scintilla 0.99 2.30 0.81 1.36 0.03 6.67
## X66.25.1 X6728.26.3 X928.95.0 X111.27.3 X123.92.2 X110.43.0
## 2012 C1 Emerald 1196.56 4345.51 355.57 190.77 0 0.72
## 2012 C1 Endura 907.13 5165.84 401.18 254.44 0 1.35
## 2012 C1 Farthing 1058.29 3146.05 358.01 166.66 0 1.61
## 2012 C1 Meadowlark 1380.10 4123.45 390.25 245.92 0 2.62
## 2012 C1 Primadonna 1766.28 4999.20 290.46 226.99 0 0.78
## 2012 C1 Scintilla 1379.83 3488.36 272.86 165.15 0 1.69
## X111.71.7 X1191.16.8 X106.70.7 X7785.70.8 X928.68.7
## 2012 C1 Emerald 0.90 0.00 0.08 0.00 0.32
## 2012 C1 Endura 0.53 0.22 0.26 0.05 0.43
## 2012 C1 Farthing 0.35 0.28 0.05 0.00 0.07
## 2012 C1 Meadowlark 0.49 0.00 0.09 0.00 0.46
## 2012 C1 Primadonna 1.02 0.00 0.25 0.00 0.31
## 2012 C1 Scintilla 0.78 0.21 0.21 0.00 0.10
## X142.62.1 X110.93.0 X111.13.7 X123.66.0 X3681.71.8 X142.92.7
## 2012 C1 Emerald 0.42 2.42 0.01 1.54 0 4.27
## 2012 C1 Endura 0.33 2.54 0.16 1.42 0 10.75
## 2012 C1 Farthing 0.37 1.53 0.00 0.04 0 1.92
## 2012 C1 Meadowlark 0.33 12.18 0.15 3.00 0 6.63
## 2012 C1 Primadonna 0.65 5.21 0.12 1.96 0 7.46
## 2012 C1 Scintilla 0.50 1.31 0.00 0.00 0 2.10
## X2497.18.9 X104.76.7 X5989.27.5 X470.82.6 X122.78.1
## 2012 C1 Emerald 0.00 0 0.79 1.87 0.13
## 2012 C1 Endura 20.05 0 4.94 10.63 0.42
## 2012 C1 Farthing 5.78 0 0.43 3.29 0.07
## 2012 C1 Meadowlark 0.00 0 4.22 1.26 0.20
## 2012 C1 Primadonna 0.00 0 2.51 1.92 0.07
## 2012 C1 Scintilla 5.75 0 0.00 0.21 0.00
## X821.55.6 X78.70.6 X124.19.6 X2639.63.6 X53398.83.7
## 2012 C1 Emerald 0.50 36.42 0.74 0.21 0.71
## 2012 C1 Endura 4.79 16.18 0.69 0.18 3.38
## 2012 C1 Farthing 0.32 10.86 0.26 0.29 1.41
## 2012 C1 Meadowlark 1.62 21.46 0.25 0.18 0.46
## 2012 C1 Primadonna 0.34 10.41 0.89 0.20 0.69
## 2012 C1 Scintilla 2.25 2.90 0.72 0.00 0.15
## X112.40.3 X119.36.8 X106.26.3 X141.27.5 X112.12.9 X629.50.5
## 2012 C1 Emerald 0.42 1.42 0.70 1.63 0.38 0.11
## 2012 C1 Endura 0.05 0.38 0.61 1.94 1.39 0.11
## 2012 C1 Farthing 0.01 0.11 0.39 0.52 0.11 0.05
## 2012 C1 Meadowlark 0.25 0.00 0.18 0.03 0.85 0.09
## 2012 C1 Primadonna 0.17 0.00 0.55 1.67 0.05 0.06
## 2012 C1 Scintilla 0.00 0.21 0.26 0.87 1.67 0.06
## X3879.26.3 X689.67.8 X1139.30.6 X5989.33.3 X43219.68.7
## 2012 C1 Emerald 0.00 2.32 1.78 0.14 1.27
## 2012 C1 Endura 0.14 3.30 5.46 1.65 0.66
## 2012 C1 Farthing 0.00 0.66 2.08 0.15 0.20
## 2012 C1 Meadowlark 0.42 14.76 4.74 0.81 1.50
## 2012 C1 Primadonna 0.00 1.52 1.65 0.83 0.41
## 2012 C1 Scintilla 0.07 2.14 4.64 0.00 0.00
## X582.16.1
## 2012 C1 Emerald 0.00
## 2012 C1 Endura 0.31
## 2012 C1 Farthing 0.22
## 2012 C1 Meadowlark 0.12
## 2012 C1 Primadonna 0.17
## 2012 C1 Scintilla 0.00
最初の10成分の変数のサマリを確認する。
summary(met1213[,1:10])
## X590.86.3 X616.25.1 X107.87.9 X110.62.3
## Min. :0.360 Min. :1.530 Min. : 0.000 Min. :2.270
## 1st Qu.:0.840 1st Qu.:4.370 1st Qu.: 0.740 1st Qu.:3.870
## Median :1.080 Median :5.010 Median : 3.790 Median :4.530
## Mean :1.296 Mean :5.098 Mean : 5.181 Mean :4.661
## 3rd Qu.:1.560 3rd Qu.:5.720 3rd Qu.: 7.100 3rd Qu.:5.290
## Max. :3.800 Max. :8.550 Max. :28.680 Max. :7.550
## X105.37.3 X623.42.7 X123.51.3 X1576.87.0
## Min. :0.0000 Min. :0.0000 Min. :0.030 Min. :1.130
## 1st Qu.:0.0000 1st Qu.:0.0800 1st Qu.:0.150 1st Qu.:2.030
## Median :0.0600 Median :0.1500 Median :0.230 Median :3.050
## Mean :0.1227 Mean :0.2025 Mean :0.335 Mean :2.964
## 3rd Qu.:0.1500 3rd Qu.:0.2700 3rd Qu.:0.380 3rd Qu.:3.680
## Max. :1.6600 Max. :0.6600 Max. :2.310 Max. :5.400
## X71.41.0 X1576.95.0
## Min. :0.4400 Min. :0.370
## 1st Qu.:0.6700 1st Qu.:0.950
## Median :0.7500 Median :1.190
## Mean :0.9136 Mean :1.201
## 3rd Qu.:1.0700 3rd Qu.:1.360
## Max. :2.6900 Max. :2.320
成分によってとりうる値の範囲が様々である。そこで、平均0、分散1に基準化する。
met1213 <- scale(met1213)
場所と遺伝子型(品種)のデータについてもattribute(属性)を表す変数として抽出しておく。
att1213 <- data %>% filter(Year == 2012 | Year == 2013) %>% select(Location, Genotype)
head(att1213)
## Location Genotype
## 2012 C1 Emerald C Emerald
## 2012 C1 Endura C Endura
## 2012 C1 Farthing C Farthing
## 2012 C1 Meadowlark C Meadowlark
## 2012 C1 Primadonna C Primadonna
## 2012 C1 Scintilla C Scintilla
まずは、基準化されたデータから距離を計算する。最初に、マンハッタン距離を計算してみる。
d <- dist(met1213, method = "manhattan")
計算された距離をもとに、完全結法(最遠隣法)で階層的クラスタ解析を行う。
tre <- hclust(d, method = "complete")
plot(tre)
ラベルに色付けするために一手間かける。
dend <- as.dendrogram(tre)
col4lab <- as.numeric(att1213$Genotype)
col4lab <- col4lab[order.dendrogram(dend)]
labels_colors(dend) <- col4lab
par(mar = c(10,4,2,2)) # 下の余白を大きく取る
plot(dend)
ユークリッド距離を計算する。
d <- dist(met1213, method = "euclidean")
完全連結法で階層的クラスタリングを行う。
dend <- as.dendrogram(hclust(d, method = "complete"))
labels_colors(dend) <- as.numeric(att1213$Genotype)[order.dendrogram(dend)]
par(mar = c(10,4,2,2)) # 下の余白を大きく取る
plot(dend)
ウォード法を適用する。
dend <- as.dendrogram(hclust(d, method = "ward.D2"))
labels_colors(dend) <- as.numeric(att1213$Genotype)[order.dendrogram(dend)]
par(mar = c(10,4,2,2)) # 下の余白を大きく取る
plot(dend)
今度は、揮発性成分のほうをクラスタにまとめてみる。
d <- dist(t(met1213)) # 転置して、距離を計算する。
ウォード法を適用する。
tre <- hclust(d, method = "ward.D2")
plot(tre)
5つのグループに分けてみる。
plot(tre)
rect.hclust(tre, k = 5)
k-平均法を適用する。kは品種数の6とする。
res <- kmeans(met1213, centers = 6, nstart = 100)
最初の10サンプルのグループを示す。
head(res$cluster, 10)
## 2012 C1 Emerald 2012 C1 Endura 2012 C1 Farthing 2012 C1 Meadowlark
## 6 5 4 1
## 2012 C1 Primadonna 2012 C1 Scintilla 2012 C2 Emerald 2012 C2 Endura
## 6 3 6 5
## 2012 C2 Farthing 2012 C2 Meadowlark
## 6 1
table関数で数え上げる。
table(att1213$Genotype, res$cluster)
##
## 1 2 3 4 5 6
## Emerald 0 0 0 0 0 18
## Endura 0 9 0 3 6 0
## Farthing 0 1 0 15 0 2
## Meadowlark 17 0 1 0 0 0
## Primadonna 0 0 0 4 0 13
## Scintilla 0 0 12 0 0 0
同じ品種がまとまって割り付けられるクラスタもある。
table関数を用いて、場所との関係も確認してみる。
table(att1213$Location, res$cluster)
##
## 1 2 3 4 5 6
## C 6 1 6 8 4 11
## H 6 5 6 7 1 10
## W 5 4 1 7 1 12
場所の違いはクラスタにあまり反映されていない。
パッケージkohonenのsom関数を用いて、自己組織化マップを作成する。
res <- som(met1213, grid = somgrid(topo = "hexagonal", xdim = 4, ydim = 4)) # 4×4のグリッドに割り振る
summary(res)
## SOM of size 4x4 with a hexagonal topology and a bubble neighbourhood function.
## The number of data layers is 1.
## Distance measure(s) used: sumofsquares.
## Training data included: 101 objects.
## Mean distance to the closest unit in the map: 23.303.
4×4のグリッドに割り付けられたデータのパターンを確認する(各成分の平均)。
plot(res, type = "codes")
ここに割り付けられたサンプルを示す。
plot(res, type = "mapping")
どれがどのサンプルか分からないので、品種の違いをラベルとして用いる。
var.id <- as.numeric(att1213$Genotype)
品種の違いを数字と色で示す。
plot(res, type = "mapping", labels = var.id, col = var.id)
講義では紹介していないが、tSNEを使って2次元でサンプル間の類似性を表示してみる。
res <- Rtsne(met1213, perplexity = 8)
plot(res$Y, pch = var.id, col = var.id)
同じ品種が2次元上でかたまり(クラスタ)を形成していることがわかる。緑(+)と水色(◇)の品種は、SOMだけでなくtSNEでも同じような混ざって散布される蛍光があることがわかる。
パッケージkernlabのksvm関数を用いて、分類モデルを作成する。
揮発性成分(met1213)から、品種(att1213$Genotype)を予測するモデル。
model <- ksvm(x = met1213, y = att1213$Genotype)
モデルを用いて、訓練データ(met1213)について分類を行う。
genotype.fit <- predict(model, met1213)
分類結果と実際の品種名の関係をtable関数で数え上げる。
table(genotype.fit, att1213$Genotype)
##
## genotype.fit Emerald Endura Farthing Meadowlark Primadonna Scintilla
## Emerald 18 0 0 0 0 0
## Endura 0 18 0 0 0 0
## Farthing 0 0 18 0 0 0
## Meadowlark 0 0 0 18 0 0
## Primadonna 0 0 0 0 17 0
## Scintilla 0 0 0 0 0 12
完璧に分類できているが、これは予測精度を表すものではない。モデル作成に用いられなかったデータで精度を確認する必要がある。
2014年のデータを精度検証用に用いる。
met14 <- data %>% filter(Year == 2014) %>% select(starts_with("X"))
att14 <- data %>% filter(Year == 2014) %>% select(Location, Genotype)
なお、2014年の成分のデータも、12年13年のデータと同じ平均と分散で基準化されなければならない。2012, 13年のデータ(met1213)には、基準化されたときの平均と標準偏差の情報が含まれているので、それを抜き出して、同じ平均と標準偏差で2014年のデータを基準化する。
met14 <- scale(met14, center = attr(met1213, "scaled:center"), scale = attr(met1213, "scaled:scale"))
基準化された2014年のデータについて分類(予測)を行う。
genotype.pred <- predict(model, met14)
予測された品種と、実際の品種間の関係をtable関数で数え上げる。
t <- table(genotype.pred, att14$Genotype)
t
##
## genotype.pred Emerald Endura Farthing Meadowlark Primadonna Scintilla
## Emerald 7 0 0 0 0 0
## Endura 0 6 0 0 0 0
## Farthing 1 0 8 0 1 2
## Meadowlark 1 3 1 9 2 1
## Primadonna 0 0 0 0 5 0
## Scintilla 0 0 0 0 0 3
品種Meadaowlarkは100%Meadowlarkと識別されるが、他の品種も多くMeadowlarkとして分類される。
数え上げられた表tをもとに正解率を計算する。
sum(diag(t)) / sum(t)
## [1] 0.76
58%のサンプルが正確に分類される(42%は分類が誤っている)。