はじめに
この資料では、これまでに解説した方法を含めて様々な方法を用いてブルベリーの嗜好性に関するデータの解析を行う。解析には、以下の手法が例として用いられている。
- 分散分析
- 混合線形モデル
- 多変量分散分析
- 主成分分析
- カーネル主成分分析
- t-SNE
- 重回帰
- PCR、PLS回帰
- 正則化回帰（ridge, LASSO）
- random forest回帰
- 正準相関分析
- PLS回帰（多変量版）
必要なパッケージ
require(here)
require(lme4)
require(lmerTest)
require(kernlab)
require(Rtsne)
require(pls)
require(glmnet)
require(ranger)
require(tuneRanger)利用するデータ
データは、Gilbertら（2015, PLoS ONE 10: e0138494） を用いる。この研究は、複数のブルーベリー品種について多環境試験（複数試験地、複数年次に渡る栽培試験）を行い、嗜好性などの官能評価の結果と、酸度、糖含量、多数の揮発性成分などの生化学的特性との関係について解析を行ったものである。 データは、同論文のSupporting Informationに含まれている。
データの読み込み
データは、同論文のSupporting Informationから得られるデータを、csvファイルとして保存したデータを準備した。
bb <- read.csv(here("data", "blueberry.csv"), stringsAsFactors = T)
head(bb)## gID Year Location Harvest Genotype Panelists Overall.Liking
## 1 2012 C1 Emerald 2012 C 1 Emerald 95 23.5
## 2 2012 C1 Endura 2012 C 1 Endura 95 17.5
## 3 2012 C1 Farthing 2012 C 1 Farthing 95 24.2
## 4 2012 C1 Meadowlark 2012 C 1 Meadowlark 95 22.3
## 5 2012 C1 Primadonna 2012 C 1 Primadonna 95 24.0
## 6 2012 C1 Scintilla 2012 C 1 Scintilla 95 24.3
## Texture Sweetness Sourness Flavor SSC TA SSC.TA pH Fructose Glucose
## 1 24.9 22.4 10.7 25.0 10.4 0.5110 20.4 3.43 48.74 49.39
## 2 25.1 19.2 22.6 27.8 10.5 0.6199 16.9 3.36 47.79 46.49
## 3 23.9 23.6 11.8 26.5 10.6 0.3064 34.6 3.55 51.31 51.66
## 4 28.7 21.5 11.2 24.4 9.2 0.2631 34.8 3.78 42.57 45.11
## 5 23.6 25.5 10.4 25.8 11.3 0.2346 48.2 3.95 56.83 56.35
## 6 23.6 24.6 14.4 24.9 12.9 0.4719 27.3 3.40 62.40 62.46
## Sucrose Total.Sugar X590.86.3 X616.25.1 X107.87.9 X110.62.3 X105.37.3
## 1 0.13 98.26 1.05 7.20 3.79 4.79 0
## 2 0.19 94.47 1.40 8.31 3.56 4.32 0
## 3 0.18 103.15 0.82 5.27 21.69 4.53 0
## 4 0.12 87.79 0.85 5.25 5.87 4.20 0
## 5 0.21 113.39 1.82 7.33 6.88 7.47 0
## 6 1.76 126.62 1.76 5.28 0.77 4.83 0
## X623.42.7 X123.51.3 X1576.87.0 X71.41.0 X1576.95.0 X108.88.3 X556.24.1
## 1 0.11 0.33 4.46 0.94 1.92 0.03 0.09
## 2 0.43 0.55 3.73 0.70 2.32 0.08 8.70
## 3 0.13 0.59 3.05 0.87 1.21 0.00 7.65
## 4 0.20 0.24 2.12 1.23 1.52 0.02 43.87
## 5 0.36 0.50 4.93 1.53 1.67 0.16 216.42
## 6 0.04 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 X111.71.7
## 1 1196.56 4345.51 355.57 190.77 0 0.72 0.90
## 2 907.13 5165.84 401.18 254.44 0 1.35 0.53
## 3 1058.29 3146.05 358.01 166.66 0 1.61 0.35
## 4 1380.10 4123.45 390.25 245.92 0 2.62 0.49
## 5 1766.28 4999.20 290.46 226.99 0 0.78 1.02
## 6 1379.83 3488.36 272.86 165.15 0 1.69 0.78
## X1191.16.8 X106.70.7 X7785.70.8 X928.68.7 X142.62.1 X110.93.0 X111.13.7
## 1 0.00 0.08 0.00 0.32 0.42 2.42 0.01
## 2 0.22 0.26 0.05 0.43 0.33 2.54 0.16
## 3 0.28 0.05 0.00 0.07 0.37 1.53 0.00
## 4 0.00 0.09 0.00 0.46 0.33 12.18 0.15
## 5 0.00 0.25 0.00 0.31 0.65 5.21 0.12
## 6 0.21 0.21 0.00 0.10 0.50 1.31 0.00
## X123.66.0 X3681.71.8 X142.92.7 X2497.18.9 X104.76.7 X5989.27.5 X470.82.6
## 1 1.54 0 4.27 0.00 0 0.79 1.87
## 2 1.42 0 10.75 20.05 0 4.94 10.63
## 3 0.04 0 1.92 5.78 0 0.43 3.29
## 4 3.00 0 6.63 0.00 0 4.22 1.26
## 5 1.96 0 7.46 0.00 0 2.51 1.92
## 6 0.00 0 2.10 5.75 0 0.00 0.21
## X122.78.1 X821.55.6 X78.70.6 X124.19.6 X2639.63.6 X53398.83.7 X112.40.3
## 1 0.13 0.50 36.42 0.74 0.21 0.71 0.42
## 2 0.42 4.79 16.18 0.69 0.18 3.38 0.05
## 3 0.07 0.32 10.86 0.26 0.29 1.41 0.01
## 4 0.20 1.62 21.46 0.25 0.18 0.46 0.25
## 5 0.07 0.34 10.41 0.89 0.20 0.69 0.17
## 6 0.00 2.25 2.90 0.72 0.00 0.15 0.00
## X119.36.8 X106.26.3 X141.27.5 X112.12.9 X629.50.5 X3879.26.3 X689.67.8
## 1 1.42 0.70 1.63 0.38 0.11 0.00 2.32
## 2 0.38 0.61 1.94 1.39 0.11 0.14 3.30
## 3 0.11 0.39 0.52 0.11 0.05 0.00 0.66
## 4 0.00 0.18 0.03 0.85 0.09 0.42 14.76
## 5 0.00 0.55 1.67 0.05 0.06 0.00 1.52
## 6 0.21 0.26 0.87 1.67 0.06 0.07 2.14
## X1139.30.6 X5989.33.3 X43219.68.7 X564.94.3 X582.16.1
## 1 1.78 0.14 1.27 0 0.00
## 2 5.46 1.65 0.66 0 0.31
## 3 2.08 0.15 0.20 0 0.22
## 4 4.74 0.81 1.50 0 0.12
## 5 1.65 0.83 0.41 0 0.17
## 6 4.64 0.00 0.00 0 0.00データの内容を確認する。
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" "X590.86.3"
## [21] "X616.25.1" "X107.87.9" "X110.62.3" "X105.37.3"
## [25] "X623.42.7" "X123.51.3" "X1576.87.0" "X71.41.0"
## [29] "X1576.95.0" "X108.88.3" "X556.24.1" "X66.25.1"
## [33] "X6728.26.3" "X928.95.0" "X111.27.3" "X123.92.2"
## [37] "X110.43.0" "X111.71.7" "X1191.16.8" "X106.70.7"
## [41] "X7785.70.8" "X928.68.7" "X142.62.1" "X110.93.0"
## [45] "X111.13.7" "X123.66.0" "X3681.71.8" "X142.92.7"
## [49] "X2497.18.9" "X104.76.7" "X5989.27.5" "X470.82.6"
## [53] "X122.78.1" "X821.55.6" "X78.70.6" "X124.19.6"
## [57] "X2639.63.6" "X53398.83.7" "X112.40.3" "X119.36.8"
## [61] "X106.26.3" "X141.27.5" "X112.12.9" "X629.50.5"
## [65] "X3879.26.3" "X689.67.8" "X1139.30.6" "X5989.33.3"
## [69] "X43219.68.7" "X564.94.3" "X582.16.1"summary(bb[,1:19])## gID Year Location Harvest
## 2012 C1 Emerald : 1 Min. :2012 A:13 Min. :1.000
## 2012 C1 Endura : 1 1st Qu.:2012 C:53 1st Qu.:1.000
## 2012 C1 Farthing : 1 Median :2013 H:53 Median :2.000
## 2012 C1 Meadowlark: 1 Mean :2013 W:45 Mean :1.994
## 2012 C1 Primadonna: 1 3rd Qu.:2014 3rd Qu.:3.000
## 2012 C1 Scintilla : 1 Max. :2014 Max. :3.000
## (Other) :158
## Genotype Panelists Overall.Liking Texture
## Emerald :27 Min. : 72.0 Min. : 6.60 Min. :12.90
## Endura :27 1st Qu.: 88.0 1st Qu.:19.07 1st Qu.:23.00
## Farthing :27 Median : 93.0 Median :22.40 Median :25.00
## Meadowlark:27 Mean : 92.3 Mean :22.28 Mean :24.83
## Primadonna:25 3rd Qu.: 95.0 3rd Qu.:25.50 3rd Qu.:27.12
## Scintilla :18 Max. :109.0 Max. :34.80 Max. :32.10
## (Other) :13
## Sweetness Sourness Flavor SSC
## Min. :12.10 Min. : 5.90 Min. :17.00 Min. : 8.10
## 1st Qu.:20.27 1st Qu.:11.20 1st Qu.:24.38 1st Qu.:10.20
## Median :22.95 Median :14.20 Median :26.00 Median :10.90
## Mean :22.93 Mean :15.56 Mean :26.09 Mean :10.93
## 3rd Qu.:25.52 3rd Qu.:19.00 3rd Qu.:28.00 3rd Qu.:11.62
## Max. :33.50 Max. :36.90 Max. :33.60 Max. :13.90
##
## TA SSC.TA pH Fructose
## Min. :0.0661 Min. : 9.50 Min. :2.820 Min. :28.40
## 1st Qu.:0.2358 1st Qu.: 19.27 1st Qu.:3.225 1st Qu.:44.59
## Median :0.3763 Median : 30.20 Median :3.465 Median :48.63
## Mean :0.4141 Mean : 34.97 Mean :3.472 Mean :48.99
## 3rd Qu.:0.5426 3rd Qu.: 47.15 3rd Qu.:3.680 3rd Qu.:53.93
## Max. :0.9849 Max. :145.20 Max. :4.410 Max. :69.86
##
## Glucose Sucrose Total.Sugar
## Min. :31.48 Min. :0.04 Min. : 60.66
## 1st Qu.:45.02 1st Qu.:0.30 1st Qu.: 90.37
## Median :48.96 Median :1.25 Median : 99.56
## Mean :49.88 Mean :1.64 Mean :100.54
## 3rd Qu.:54.97 3rd Qu.:2.38 3rd Qu.:109.54
## Max. :76.01 Max. :5.69 Max. :140.36
## 年次、収穫時期、パネリストのIDが数値データとなっているので、factor（因子）に変換する。
bb$Year <- as.factor(bb$Year)
bb$Harvest <- as.factor(bb$Harvest)
bb$Panelists <- as.factor(bb$Panelists)
summary(bb[,1:19])## gID Year Location Harvest Genotype
## 2012 C1 Emerald : 1 2012:51 A:13 1:55 Emerald :27
## 2012 C1 Endura : 1 2013:50 C:53 2:55 Endura :27
## 2012 C1 Farthing : 1 2014:63 H:53 3:54 Farthing :27
## 2012 C1 Meadowlark: 1 W:45 Meadowlark:27
## 2012 C1 Primadonna: 1 Primadonna:25
## 2012 C1 Scintilla : 1 Scintilla :18
## (Other) :158 (Other) :13
## Panelists Overall.Liking Texture Sweetness Sourness
## 95 :29 Min. : 6.60 Min. :12.90 Min. :12.10 Min. : 5.90
## 90 :28 1st Qu.:19.07 1st Qu.:23.00 1st Qu.:20.27 1st Qu.:11.20
## 93 :17 Median :22.40 Median :25.00 Median :22.95 Median :14.20
## 85 :11 Mean :22.28 Mean :24.83 Mean :22.93 Mean :15.56
## 101 :11 3rd Qu.:25.50 3rd Qu.:27.12 3rd Qu.:25.52 3rd Qu.:19.00
## 108 :10 Max. :34.80 Max. :32.10 Max. :33.50 Max. :36.90
## (Other):58
## Flavor SSC TA SSC.TA
## Min. :17.00 Min. : 8.10 Min. :0.0661 Min. : 9.50
## 1st Qu.:24.38 1st Qu.:10.20 1st Qu.:0.2358 1st Qu.: 19.27
## Median :26.00 Median :10.90 Median :0.3763 Median : 30.20
## Mean :26.09 Mean :10.93 Mean :0.4141 Mean : 34.97
## 3rd Qu.:28.00 3rd Qu.:11.62 3rd Qu.:0.5426 3rd Qu.: 47.15
## Max. :33.60 Max. :13.90 Max. :0.9849 Max. :145.20
##
## pH Fructose Glucose Sucrose
## Min. :2.820 Min. :28.40 Min. :31.48 Min. :0.04
## 1st Qu.:3.225 1st Qu.:44.59 1st Qu.:45.02 1st Qu.:0.30
## Median :3.465 Median :48.63 Median :48.96 Median :1.25
## Mean :3.472 Mean :48.99 Mean :49.88 Mean :1.64
## 3rd Qu.:3.680 3rd Qu.:53.93 3rd Qu.:54.97 3rd Qu.:2.38
## Max. :4.410 Max. :69.86 Max. :76.01 Max. :5.69
##
## Total.Sugar
## Min. : 60.66
## 1st Qu.: 90.37
## Median : 99.56
## Mean :100.54
## 3rd Qu.:109.54
## Max. :140.36
## 3年間4試験地で3回収穫されて計測されたデータであることが分かる。
なお、後半（20列目以降）のデータは、ガスクロマトグラフで計測された揮発性成分の計測値であり、それぞれの具体的な物質名は、以下のファイルに保存されている。
vc <- read.csv(here("data", "volatileComp.csv"))
head(vc)## CAS_id Name Formula
## 1 590-86-3 3-Methyl Butanal C5 H10 O
## 2 616-25-1 1-Penten-3-ol C5 H10 O
## 3 107-87-9 2-Pentanone C5 H10 O
## 4 110-62-3 Pentanal C5 H10 O
## 5 105-37-3 Propanoic acid, ethyl ester C5 H10 O2
## 6 623-42-7 Butanoic acid, methyl ester C5 H10 O2どの品種が、どの環境（年次、試験地）で栽培試験されているかを集計表を作って確認する。
table(paste(bb$Year, bb$Location), bb$Genotype)##
## Chickadee Emerald Endura Farthing FL01-25 FL06-435 FL06-510 FL06-571
## 2012 C 0 3 3 3 0 0 0 0
## 2012 H 0 3 3 3 0 0 0 0
## 2012 W 0 3 3 3 0 0 0 0
## 2013 C 0 3 3 3 0 0 0 0
## 2013 H 0 3 3 3 0 0 0 0
## 2013 W 0 3 3 3 0 0 0 0
## 2014 A 1 0 0 0 1 1 1 1
## 2014 C 0 3 3 3 0 0 0 0
## 2014 H 0 3 3 3 0 0 0 0
## 2014 W 0 3 3 3 0 0 0 0
##
## FL07-290 FL08-022 FL10-107 FL10-186 Flicker Kestrel Meadowlark
## 2012 C 0 0 0 0 0 0 3
## 2012 H 0 0 0 0 0 0 3
## 2012 W 0 0 0 0 0 0 3
## 2013 C 0 0 0 0 0 0 3
## 2013 H 0 0 0 0 0 0 3
## 2013 W 0 0 0 0 0 0 3
## 2014 A 1 1 1 1 1 1 0
## 2014 C 0 0 0 0 0 0 3
## 2014 H 0 0 0 0 0 0 3
## 2014 W 0 0 0 0 0 0 3
##
## Primadonna Scintilla Vireo Windsor
## 2012 C 3 3 0 0
## 2012 H 3 3 0 0
## 2012 W 3 0 0 0
## 2013 C 3 3 0 0
## 2013 H 2 3 0 0
## 2013 W 3 0 0 0
## 2014 A 0 0 1 1
## 2014 C 2 3 0 0
## 2014 H 3 3 0 0
## 2014 W 3 0 0 06品種が複数の環境で栽培されており、残りの13品種は、1環境でのみ栽培されていることが分かる。 ここでは、まず、全試験地で栽培された5品種のデータにしぼって解析をすることとする。
selector <- bb$Genotype %in% c("Emerald", "Endura", "Farthing", "Meadowlark", "Primadonna")
bb5 <- bb[selector, ]遺伝子型（G）、環境（E）、遺伝子型と環境の交互作用（GxE）の効果
嗜好性のスコアや、生化学物質の計測値に関して遺伝子型や環境がどのように影響しているのかを確認してみる。
y <- bb5$SSC
df <- data.frame(y, G = bb5$Genotype, E = factor(paste0(bb5$Year, bb5$Location)))
model <- lm(y ~ G * E, data = df)
anova(model)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## G 4 47.417 11.8542 25.3136 5.716e-14 ***
## E 8 24.939 3.1174 6.6569 8.802e-07 ***
## G:E 32 31.265 0.9770 2.0863 0.003704 **
## Residuals 88 41.210 0.4683
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1可溶性固形物含量(soluble solid content: SSC) では、遺伝子型と交互作用が1%水準で有意である。
df$y <- bb5$TA
model <- lm(y ~ G * E, data = df)
anova(model)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## G 4 1.2714 0.31785 17.6858 1.098e-10 ***
## E 8 2.3740 0.29675 16.5119 1.113e-14 ***
## G:E 32 1.0043 0.03138 1.7462 0.02164 *
## Residuals 88 1.5815 0.01797
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1滴定酸度（titratable acidity: TA）では、SSCに比べると交互作用項の有意性が少し低い。
嗜好性についても検定してみる。
df$y <- bb5$Overall.Liking
model <- lm(y ~ G * E, data = df)
anova(model)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## G 4 516.38 129.094 15.1620 1.819e-09 ***
## E 8 245.42 30.678 3.6031 0.001156 **
## G:E 32 611.91 19.122 2.2459 0.001575 **
## Residuals 88 749.26 8.514
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1交互作用が有意である。なお、パネリストの効果を変量効果として入れた、混合モデル（母数効果と変量効果を含んだモデル）を適用してみる。
df2 <- data.frame(df, P = bb5$Panelists)
model <- lmer(y ~ G * E + (1 | P), data = df2)
anova(model)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## G 517.12 129.280 4 76.015 17.1829 4.455e-10 ***
## E 141.14 17.643 8 46.115 2.3449 0.0331928 *
## G:E 611.54 19.111 32 75.977 2.5400 0.0004746 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1環境による効果が小さくなった。ただし、混合モデルにおける母数効果の検定は様々な議論があり、解釈には注意が必要である。以降、パネリストの違いは、考慮せずに解析をおこなうこととする。
mlogp <- NULL
col.id <- 7:71
for(i in col.id) {
df$y <- bb5[, i]
model <- lm(y ~ G * E, data = df)
mlogp <- cbind(mlogp, -log(anova(model)$Pr[1:3]))
}
colnames(mlogp) <- colnames(bb5)[col.id]
op <- par(mar = c(10, 5, 2, 2))
barplot(mlogp, beside = T, las = 2, legend = c("G", "E", "G:E"), ylab = "-log(p)")
par(op)変数によって、G,E,GxEの効果の有意性が様々であることが分かる。
多変量分散分析
複数の変数をまとめて多変量分散分析（MANOVA）を行って見る。
colnames(bb5)## [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" "X590.86.3"
## [21] "X616.25.1" "X107.87.9" "X110.62.3" "X105.37.3"
## [25] "X623.42.7" "X123.51.3" "X1576.87.0" "X71.41.0"
## [29] "X1576.95.0" "X108.88.3" "X556.24.1" "X66.25.1"
## [33] "X6728.26.3" "X928.95.0" "X111.27.3" "X123.92.2"
## [37] "X110.43.0" "X111.71.7" "X1191.16.8" "X106.70.7"
## [41] "X7785.70.8" "X928.68.7" "X142.62.1" "X110.93.0"
## [45] "X111.13.7" "X123.66.0" "X3681.71.8" "X142.92.7"
## [49] "X2497.18.9" "X104.76.7" "X5989.27.5" "X470.82.6"
## [53] "X122.78.1" "X821.55.6" "X78.70.6" "X124.19.6"
## [57] "X2639.63.6" "X53398.83.7" "X112.40.3" "X119.36.8"
## [61] "X106.26.3" "X141.27.5" "X112.12.9" "X629.50.5"
## [65] "X3879.26.3" "X689.67.8" "X1139.30.6" "X5989.33.3"
## [69] "X43219.68.7" "X564.94.3" "X582.16.1"変数群のIDリストを作っておく。
liking.id <- 7:11
sugaci.id <- c(12, 13, 15:18)
volati.id <- 20:71まずは、嗜好性の変数群についてMANOVAを行なってみる。
var.id <- liking.id
x <- as.matrix(bb5[, var.id])
res <- manova(x ~ bb5$Genotype * bb5$Location * bb5$Year)
summary(res)## Df Pillai approx F num Df den Df Pr(>F)
## bb5$Genotype 4 1.50466 10.4919 20 348 < 2.2e-16
## bb5$Location 2 0.62133 7.6614 10 170 4.671e-10
## bb5$Year 2 0.31731 3.2057 10 170 0.0008423
## bb5$Genotype:bb5$Location 8 0.69364 1.7718 40 440 0.0033372
## bb5$Genotype:bb5$Year 8 0.69926 1.7885 40 440 0.0028934
## bb5$Location:bb5$Year 4 0.29093 1.3648 20 348 0.1366548
## bb5$Genotype:bb5$Location:bb5$Year 16 1.32291 1.9787 80 440 7.933e-06
## Residuals 88
##
## bb5$Genotype ***
## bb5$Location ***
## bb5$Year ***
## bb5$Genotype:bb5$Location **
## bb5$Genotype:bb5$Year **
## bb5$Location:bb5$Year
## bb5$Genotype:bb5$Location:bb5$Year ***
## Residuals
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1試験地と年次の交互作用を除いて、高度に有意である。
var.id <- sugaci.id
x <- as.matrix(bb5[, var.id])
res <- manova(x ~ bb5$Genotype * bb5$Location * bb5$Year)
summary(res)## Df Pillai approx F num Df den Df Pr(>F)
## bb5$Genotype 4 1.74449 11.0859 24 344 < 2.2e-16
## bb5$Location 2 0.96387 13.0237 12 168 < 2.2e-16
## bb5$Year 2 1.30266 26.1527 12 168 < 2.2e-16
## bb5$Genotype:bb5$Location 8 0.85922 1.8385 48 528 0.0007774
## bb5$Genotype:bb5$Year 8 0.76863 1.6162 48 528 0.0068833
## bb5$Location:bb5$Year 4 1.09507 5.4033 24 344 1.404e-13
## bb5$Genotype:bb5$Location:bb5$Year 16 1.18622 1.3553 96 528 0.0206837
## Residuals 88
##
## bb5$Genotype ***
## bb5$Location ***
## bb5$Year ***
## bb5$Genotype:bb5$Location ***
## bb5$Genotype:bb5$Year **
## bb5$Location:bb5$Year ***
## bb5$Genotype:bb5$Location:bb5$Year *
## Residuals
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1糖度・酸度については、全ての効果が有意である。嗜好性に比べて、年次効果が特に大きい。
var.id <- volati.id
x <- as.matrix(bb5[, var.id])
res <- manova(x ~ bb5$Genotype * bb5$Location * bb5$Year)
summary(res)## Df Pillai approx F num Df den Df Pr(>F)
## bb5$Genotype 4 3.9197 37.556 208 160 < 2.2e-16
## bb5$Location 2 1.7988 6.532 104 76 1.348e-15
## bb5$Year 2 1.9425 24.689 104 76 < 2.2e-16
## bb5$Genotype:bb5$Location 8 6.0346 2.598 416 352 < 2.2e-16
## bb5$Genotype:bb5$Year 8 6.1303 2.774 416 352 < 2.2e-16
## bb5$Location:bb5$Year 4 3.1676 2.927 208 160 2.706e-12
## bb5$Genotype:bb5$Location:bb5$Year 16 9.3054 1.390 832 832 1.085e-06
## Residuals 88
##
## bb5$Genotype ***
## bb5$Location ***
## bb5$Year ***
## bb5$Genotype:bb5$Location ***
## bb5$Genotype:bb5$Year ***
## bb5$Location:bb5$Year ***
## bb5$Genotype:bb5$Location:bb5$Year ***
## Residuals
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1揮発性成分については、いずれの効果も高度に有意である。遺伝子型の効果は、3つの変数群の中で最も大きい。 3つの変数群の中で、最も遺伝的な支配が大きな変数群である。このことは、分散分析の結果の棒グラフからも見て取れる。
多次元データの視覚化
多次元データの視覚化のために、主成分分析をおこなってみる。ここでは、対象を生化学成分とする。
col.id <- c(sugaci.id, volati.id)
x <- scale(bb5[, col.id])res <- prcomp(x)plot(res$x, col = bb5$Genotype, pch = as.numeric(bb5$Location), asp = 1)
品種間の差のほうが、栽培地間の差に比べて顕著である。
各主成分と元の変数の関係は、biplot関数により視覚化できる。
biplot(res)
次に、カーネル主成分分析を行なって、非線形の主成分分析を行ってみる。
res <- kpca(x, kernel = "rbfdot", kpar = list(sigma = 0.05), features = 2)
plot(rotated(res), col = bb5$Genotype, pch = as.numeric(bb5$Location), asp = 1)
結果は、強くハイパーパラメータであるsigmaに依存する。sigmaを小さくすると結果は線形の主成分分析に近づく。
res <- kpca(x, kernel = "rbfdot", kpar = list(sigma = 0.001), features = 2)
plot(rotated(res), col = bb5$Genotype, pch = as.numeric(bb5$Location), asp = 1)
様々な、カーネルを指定することができる。
res <- kpca(x, kernel = "laplacedot", kpar = list(sigma = 0.05), features = 2)
plot(rotated(res), col = bb5$Genotype, pch = as.numeric(bb5$Location), asp = 1)
次に、最近、多次元データの視覚化法としてよく用いられるようになったt-SNEを適用してみる。
res <- Rtsne(x, perplexity = 30)
plot(res$Y, col = bb5$Genotype, pch = as.numeric(bb5$Location), asp = 1)
t-SNEでは近い関係がより近く、遠い関係がより遠く表されるため、まとまりが明瞭になる。
なお、t-SNEはハイパーパラメータperplexityに結果が依存する。
res <- Rtsne(x, perplexity = 5)
plot(res$Y, col = bb5$Genotype, pch = as.numeric(bb5$Location), asp = 1)
嗜好性と生化学物質の計測値の関係
まずは、揮発性成分を除いた生化学物質と、嗜好性の各項目の関係を見てみる。
var.id <- sugaci.id
x <- scale(as.matrix(bb5[, var.id]))
y <- bb5$Overall.Liking
model <- lm(y ~ x)
summary(model)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.6897 -2.0125 0.0379 2.2110 6.1675
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.8128 0.2754 79.202 < 2e-16 ***
## xSSC 0.7109 0.4720 1.506 0.1345
## xTA -2.1561 0.4443 -4.853 3.52e-06 ***
## xpH -0.8204 0.4513 -1.818 0.0714 .
## xFructose 1.5294 0.5982 2.557 0.0118 *
## xGlucose -0.6217 0.5027 -1.237 0.2184
## xSucrose 0.2924 0.3224 0.907 0.3662
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.176 on 126 degrees of freedom
## Multiple R-squared: 0.4013, Adjusted R-squared: 0.3728
## F-statistic: 14.07 on 6 and 126 DF, p-value: 3.242e-12TAと果糖の効果が有意である。
揮発性物質との関係をモデル化してみる。
var.id <- volati.id
x <- scale(as.matrix(bb5[, var.id]))
model <- lm(y ~ x)
summary(model)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.2288 -1.3666 0.1629 1.5708 6.3325
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.81278 0.26915 81.044 <2e-16 ***
## xX590.86.3 0.29672 0.50491 0.588 0.5584
## xX616.25.1 0.91515 1.20919 0.757 0.4514
## xX107.87.9 -0.05995 0.35526 -0.169 0.8664
## xX110.62.3 0.12609 1.51790 0.083 0.9340
## xX105.37.3 0.20497 0.45604 0.449 0.6543
## xX623.42.7 -1.06166 0.56659 -1.874 0.0646 .
## xX123.51.3 0.15649 0.50861 0.308 0.7591
## xX1576.87.0 0.32023 1.21568 0.263 0.7929
## xX71.41.0 -0.28374 0.69514 -0.408 0.6842
## xX1576.95.0 -0.83877 1.15487 -0.726 0.4698
## xX108.88.3 0.03352 0.60961 0.055 0.9563
## xX556.24.1 -0.68992 0.55108 -1.252 0.2142
## xX66.25.1 0.94831 0.98935 0.959 0.3407
## xX6728.26.3 0.43513 1.26722 0.343 0.7322
## xX928.95.0 -0.17040 0.86585 -0.197 0.8445
## xX111.27.3 0.70247 1.01390 0.693 0.4904
## xX123.92.2 -0.17341 0.69435 -0.250 0.8034
## xX110.43.0 0.53180 0.63439 0.838 0.4044
## xX111.71.7 -0.79001 0.80056 -0.987 0.3267
## xX1191.16.8 -0.03961 0.55886 -0.071 0.9437
## xX106.70.7 -0.26939 0.58623 -0.460 0.6471
## xX7785.70.8 0.53328 0.52080 1.024 0.3089
## xX928.68.7 -0.35116 0.47190 -0.744 0.4590
## xX142.62.1 -0.26962 0.80312 -0.336 0.7380
## xX110.93.0 -0.55934 1.22160 -0.458 0.6483
## xX111.13.7 -0.15298 0.48245 -0.317 0.7520
## xX123.66.0 0.07187 0.46017 0.156 0.8763
## xX3681.71.8 0.19271 0.70452 0.274 0.7852
## xX142.92.7 0.15113 0.60533 0.250 0.8035
## xX2497.18.9 0.08936 0.95737 0.093 0.9259
## xX104.76.7 0.07634 0.36604 0.209 0.8353
## xX5989.27.5 1.61541 0.64210 2.516 0.0139 *
## xX470.82.6 -1.55693 0.83216 -1.871 0.0650 .
## xX122.78.1 -0.72708 0.43746 -1.662 0.1004
## xX821.55.6 -0.41033 1.28308 -0.320 0.7500
## xX78.70.6 -0.84245 0.61073 -1.379 0.1716
## xX124.19.6 -0.63313 0.41115 -1.540 0.1275
## xX2639.63.6 -0.02487 0.50357 -0.049 0.9607
## xX53398.83.7 0.64723 0.69537 0.931 0.3548
## xX112.40.3 -0.29379 0.53593 -0.548 0.5851
## xX119.36.8 -0.73574 0.46410 -1.585 0.1168
## xX106.26.3 -0.11644 0.58489 -0.199 0.8427
## xX141.27.5 1.61525 0.69454 2.326 0.0226 *
## xX112.12.9 -1.63931 1.14669 -1.430 0.1567
## xX629.50.5 0.30676 0.51773 0.593 0.5552
## xX3879.26.3 0.74214 0.64683 1.147 0.2547
## xX689.67.8 -0.34328 1.67331 -0.205 0.8380
## xX1139.30.6 0.15300 1.18997 0.129 0.8980
## xX5989.33.3 -0.25775 0.66956 -0.385 0.7013
## xX43219.68.7 0.44276 0.72634 0.610 0.5439
## xX564.94.3 0.51015 0.49403 1.033 0.3049
## xX582.16.1 0.12918 0.49959 0.259 0.7966
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.104 on 80 degrees of freedom
## Multiple R-squared: 0.6369, Adjusted R-squared: 0.4009
## F-statistic: 2.699 on 52 and 80 DF, p-value: 3.121e-05モデルの説明力（\(R^2\)）は先のモデルより高くなったが、変数の効果は不明瞭である。
揮発性成分とそれ以外の両者でモデルを作成してみる。
var.id <- c(sugaci.id, volati.id)
x <- scale(as.matrix(bb5[, var.id]))
model <- lm(y ~ x)
summary(model)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.493 -1.398 -0.009 1.410 5.096
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.81278 0.23002 94.830 <2e-16 ***
## xSSC 1.26672 0.67596 1.874 0.0649 .
## xTA -1.89893 0.76387 -2.486 0.0152 *
## xpH -0.38238 0.56555 -0.676 0.5011
## xFructose 0.87995 0.70211 1.253 0.2140
## xGlucose -0.40115 0.65716 -0.610 0.5435
## xSucrose 0.12370 0.44335 0.279 0.7810
## xX590.86.3 0.22200 0.44889 0.495 0.6224
## xX616.25.1 1.48883 1.12455 1.324 0.1896
## xX107.87.9 0.18173 0.31096 0.584 0.5607
## xX110.62.3 -0.27564 1.31870 -0.209 0.8350
## xX105.37.3 0.05080 0.39269 0.129 0.8974
## xX623.42.7 -0.87221 0.49368 -1.767 0.0814 .
## xX123.51.3 0.50988 0.44626 1.143 0.2569
## xX1576.87.0 -0.44296 1.08022 -0.410 0.6829
## xX71.41.0 -0.13612 0.59964 -0.227 0.8211
## xX1576.95.0 -1.46656 1.01929 -1.439 0.1544
## xX108.88.3 -0.23717 0.53394 -0.444 0.6582
## xX556.24.1 -0.74131 0.48778 -1.520 0.1328
## xX66.25.1 0.34822 0.95559 0.364 0.7166
## xX6728.26.3 0.95230 1.09936 0.866 0.3892
## xX928.95.0 1.28869 0.81365 1.584 0.1175
## xX111.27.3 -1.46283 0.98185 -1.490 0.1405
## xX123.92.2 -0.05199 0.64471 -0.081 0.9359
## xX110.43.0 0.09477 0.59218 0.160 0.8733
## xX111.71.7 0.58006 0.83139 0.698 0.4876
## xX1191.16.8 -0.06604 0.49720 -0.133 0.8947
## xX106.70.7 0.22417 0.52647 0.426 0.6715
## xX7785.70.8 0.08719 0.47438 0.184 0.8547
## xX928.68.7 -0.47239 0.42631 -1.108 0.2714
## xX142.62.1 0.44113 0.70922 0.622 0.5359
## xX110.93.0 -1.10859 1.12055 -0.989 0.3257
## xX111.13.7 0.30468 0.42370 0.719 0.4743
## xX123.66.0 -0.15350 0.43365 -0.354 0.7244
## xX3681.71.8 -0.13605 0.62994 -0.216 0.8296
## xX142.92.7 0.41879 0.53068 0.789 0.4325
## xX2497.18.9 -0.74157 0.85783 -0.864 0.3901
## xX104.76.7 0.27586 0.34572 0.798 0.4275
## xX5989.27.5 1.31042 0.60509 2.166 0.0336 *
## xX470.82.6 -1.17422 0.73829 -1.590 0.1160
## xX122.78.1 -0.59947 0.39854 -1.504 0.1368
## xX821.55.6 -0.88036 1.14504 -0.769 0.4444
## xX78.70.6 -0.23771 0.61748 -0.385 0.7014
## xX124.19.6 -0.51440 0.37057 -1.388 0.1693
## xX2639.63.6 0.07223 0.43743 0.165 0.8693
## xX53398.83.7 1.19620 0.62573 1.912 0.0598 .
## xX112.40.3 -0.57383 0.48254 -1.189 0.2382
## xX119.36.8 -0.90976 0.41648 -2.184 0.0321 *
## xX106.26.3 -0.51737 0.51055 -1.013 0.3142
## xX141.27.5 1.32295 0.62410 2.120 0.0374 *
## xX112.12.9 -0.38672 1.08158 -0.358 0.7217
## xX629.50.5 0.49319 0.46416 1.063 0.2914
## xX3879.26.3 0.33096 0.57419 0.576 0.5661
## xX689.67.8 0.83367 1.51358 0.551 0.5834
## xX1139.30.6 0.80390 1.05584 0.761 0.4489
## xX5989.33.3 -0.72329 0.60209 -1.201 0.2335
## xX43219.68.7 -0.04831 0.69435 -0.070 0.9447
## xX564.94.3 0.33147 0.46119 0.719 0.4746
## xX582.16.1 0.52561 0.50862 1.033 0.3048
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.653 on 74 degrees of freedom
## Multiple R-squared: 0.7547, Adjusted R-squared: 0.5625
## F-statistic: 3.926 on 58 and 74 DF, p-value: 2.545e-08モデルの説明力（\(R^2\)）はさらに上がる。
モデルを当てはめた値と、観察値間の関係を図示してみる。
y.fit <- fitted(model)
range <- range(c(y, y.fit))
plot(y, y.fit, xlim = range, ylim = range)
abline(0, 1)
よく当てはまっているように見える。
cor(y, y.fit)## [1] 0.8687432両者の相関も高い（なお、この値の2乗がモデルの説明力\(R^2\)となる）。
しかし、モデルは柔軟なほど（パラメータが多いほど）当てはまりがよくなるため、特に変数の数が大きく異る場合は当てはまりの良さでは評価できない。
特に、得られたモデルを予測に利用したい場合は、予測値と観察値の一致度を見る必要がある。
ここでは、まず、10分割交差検証を用いて評価をしてみる。
n.fold <- 10
id <- 1:length(y) %% 10 + 1 # 1〜10の番号を割り振る
head(id, 20)## [1] 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 1id.rand <- sample(id) # 無作為に並べ直し
head(id.rand, 20)## [1] 4 5 2 8 4 2 10 7 2 4 1 9 5 8 9 5 10 2 4 3df <- data.frame(y, x)
y.pred <- rep(NA, length(y))
for(i in 1:n.fold) {
test <- id.rand == i
df.train <- df[!test, ]
df.test <- df[test, ]
model <- lm(y ~ ., data = df.train)
y.pred[test] <- predict(model, df.test)
}
range <- range(c(y, y.pred))
plot(y, y.pred, xlim = range, ylim = range)
abline(0, 1)
cor(y, y.pred)## [1] 0.4991842予測値と観察値の一致度は、モデルを当てはめた値と観察値の一致度に比べて下がることが分かる。 未知のデータに対して予測を行っているので、当然と言えば当然である。
予測モデルを構築する場合は、予測値と観察値の一致度を指標として、モデルの改良をする。
なお、上の例では、データを無作為に10分割して予測を行ったが、実用した場合の精度を評価するために、2年分のデータをもとにモデルをつくり、3年目の予測を行って精度評価をおこなってみる。
df <- data.frame(y, x)
df.train <- df[bb5$Year != "2014", ]
df.test <- df[bb5$Year == "2014", ]2012年, 2013年のデータをモデルの訓練に、2014年のデータをモデルの精度評価に用いる。
model <- lm(y ~ ., data = df.train)
y.pred <- predict(model, df.test)## Warning in predict.lm(model, df.test): prediction from a rank-deficient fit may
## be misleadingエラーがでる。
apply(df.train, 2, var)## y SSC TA pH Fructose Glucose
## 13.07539837 1.06478711 1.06609644 0.97621850 1.19333914 1.34879180
## Sucrose X590.86.3 X616.25.1 X107.87.9 X110.62.3 X105.37.3
## 1.11149175 0.80602318 0.53259247 1.35503218 0.31792756 0.35705280
## X623.42.7 X123.51.3 X1576.87.0 X71.41.0 X1576.95.0 X108.88.3
## 0.78963122 0.44784902 0.79483382 0.62568368 0.84964473 0.06382981
## X556.24.1 X66.25.1 X6728.26.3 X928.95.0 X111.27.3 X123.92.2
## 1.26560355 0.74070165 0.87601415 0.35355445 0.37015974 0.45841576
## X110.43.0 X111.71.7 X1191.16.8 X106.70.7 X7785.70.8 X928.68.7
## 0.60041088 0.81575617 0.79989700 0.46499731 0.23936101 1.31441753
## X142.62.1 X110.93.0 X111.13.7 X123.66.0 X3681.71.8 X142.92.7
## 0.08733661 1.10550354 1.11051897 1.27946324 0.10800142 0.48824474
## X2497.18.9 X104.76.7 X5989.27.5 X470.82.6 X122.78.1 X821.55.6
## 0.31458459 1.41056016 0.79965607 0.80038359 0.24170470 0.84890075
## X78.70.6 X124.19.6 X2639.63.6 X53398.83.7 X112.40.3 X119.36.8
## 0.85084834 0.48866726 1.09871404 0.61986529 0.86014613 1.13770053
## X106.26.3 X141.27.5 X112.12.9 X629.50.5 X3879.26.3 X689.67.8
## 0.75866431 1.01272900 0.89078859 0.90395080 1.34341149 1.08493739
## X1139.30.6 X5989.33.3 X43219.68.7 X564.94.3 X582.16.1
## 0.97504654 0.82861406 1.05275007 0.00000000 0.94961736どうやら、1つの変数が訓練データ内に変動が無いようである。除くこととする。
df.train <- df.train[, colnames(df) != "X564.94.3"]
df.test <- df.test[, colnames(df) != "X564.94.3"]model <- lm(y ~ ., data = df.train)
y.pred <- predict(model, df.test)エラーは出ない。
range <- range(c(df.test$y, y.pred))
plot(df.test$y, y.pred, xlim = range, ylim = range)
abline(0, 1)
cor(df.test$y, y.pred)## [1] 0.3836946予測精度はさらに下がることが分かる。 では、ここから別の方法を適用して、予測精度を改良していこう。
多変量回帰の利用
まず、主成分回帰分析（principal component regression: PCR）を行ってみる。PCRでは、xの主成分分析を行っておき、その主成分得点とyの回帰分析を行う。主成分分析を行うことにより、回帰分析に用いられる説明変数の次元を下げるとともに、多重共線性を取り除く。
n.comp <- 30
model <- pcr(y ~ ., n.comp, data = df.train)
summary(model)## Data: X dimension: 89 57
## Y dimension: 89 1
## Fit method: svdpc
## Number of components considered: 30
## TRAINING: % variance explained
## 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps 7 comps 8 comps
## X 19.027 32.671 42.55 51.14 57.44 62.60 67.47 71.00
## y 4.424 8.145 15.88 18.80 19.92 33.52 33.66 37.65
## 9 comps 10 comps 11 comps 12 comps 13 comps 14 comps 15 comps
## X 74.25 76.75 78.99 81.03 82.90 84.60 86.18
## y 37.70 37.95 39.67 40.37 40.96 41.34 41.46
## 16 comps 17 comps 18 comps 19 comps 20 comps 21 comps 22 comps
## X 87.50 88.75 89.78 90.76 91.64 92.47 93.26
## y 41.52 43.61 45.91 45.97 46.97 47.06 51.82
## 23 comps 24 comps 25 comps 26 comps 27 comps 28 comps 29 comps
## X 93.96 94.61 95.20 95.71 96.16 96.58 96.91
## y 55.59 56.39 58.47 58.50 58.77 62.04 63.52
## 30 comps
## X 97.23
## y 65.40主成分数が増えるとともにx、yの変動の多くの部分が説明されることがわかる。では、主成分数はどのように決めると良いだろうか？ここでは、やはり交差検証を用いて主成分数を決めてみる。あらかじめ少し多めの主成分数を設定しておき、交差検証をオプションをつけて改めてpcr関数を実行してみる。
model <- pcr(y ~ ., n.comp + 20, data = df.train, validation = "LOO")
plot(model, "validation")
描かれるグラフは、予測におけるRMSE（Root Mean Squared Error of prediction: RMSEP）である。グラフから主成分数は30程度が良いことがわかる。なお、赤い破線で示されているのは歪み（bias）補正済みのRMSEPである。同じく30主成分程度で良い予測精度が最大になることがわかる。
最もRMSEPが小さくなる主成分数は、以下のようにすると結果から取り出すことができる。
n.comp <- which.min(model$validation$PRESS) # 予測残差平方和が最小になる主成分数
n.comp## [1] 8予測に最適と判断された主成分数でPCRを行い、yの予測値と観察値の関係をみてみる。
y.pred <- predict(model, newdata = df.test, ncomp = n.comp)range <- range(c(df.test$y, y.pred))
plot(df.test$y, y.pred, xlim = range, ylim = range)
abline(0, 1)
cor(df.test$y, y.pred)## [1] 0.7321363重回帰モデルより、高い精度のモデルができる。
次に、PLS回帰分析を用いて回帰モデルを構築してみる。
n.comp <- 30
model <- plsr(y ~ ., n.comp, data = df.train, validation = "LOO")
plot(model, "validation")
PLSでは、PCRに比べ少ない数の潜在変数で予測できる。
n.comp <- which.min(model$validation$PRESS)
n.comp## [1] 2y.pred <- predict(model, newdata = df.test, ncomp = n.comp)range <- range(c(df.test$y, y.pred))
plot(df.test$y, y.pred, xlim = range, ylim = range)
abline(0, 1)
cor(df.test$y, y.pred)## [1] 0.729506正則化回帰の利用
次に、リッジ（ridge）回帰を用いて予測モデルを構築してみる。リッジ回帰は、glmnetという関数を用いて実行できる。オプションのalphaを0にしたときがリッジ回帰、1にするとLASSO、0-1の間の値はelastic netとよばれる手法に対応する。
y <- df.train$y
x <- as.matrix(df.train[, -1])
model <- glmnet(x, y, alpha = 0) # リッジ回帰モデル得られた回帰係数の図示をする。
plot(model, "lambda", label = T)
リッジ回帰では、講義でも説明したように、回帰係数の大きさに対してペナルティをかけることで訓練データに対する過度なあてはめ（over-fitting）を防ぎ、未知データに対する予測精度を向上させる。回帰係数の大きさに対するペナルティの大きさは、\(\lambda\)というパラメータによってコントロールされている。上の図は、\(\lambda\)が大きくなるに従って、回帰係数の大きさに対するペナルティが強くなっていき、最後は0に縮合（shrink）してしまう様子が示されている。ペナルティの強さを決める方法には様々あるが、ここでは、PCRやPLSの成分数を決めたのと同様に、交差検証を用いる。
関数cv.glmnetにより、交差検証で\(\lambda\)の最適化ができる。
model <- cv.glmnet(x, y, alpha = 0, nfolds = length(y), grouped = F)
plot(model)
RMSEPを最小にする\(\lambda\)を用いた場合の\(y\)の予測値を求めて、精度評価をする。
x.test <- as.matrix(df.test[, -1])
y.pred <- predict(model, newx = x.test, s = "lambda.min")range <- range(c(df.test$y, y.pred))
plot(df.test$y, y.pred, xlim = range, ylim = range)
abline(0, 1)
cor(df.test$y, y.pred)## lambda.min
## [1,] 0.6693491PCRやPLSほどは精度が出ていないことが分かる。
次に、LASSOを用いたモデル構築を行ってみる。
model <- glmnet(x, y, alpha = 1) # LASSOモデル
plot(model, "lambda", label = T)
リッジ回帰と同じく、LASSOも回帰係数の大きさにペナルティをかける。リッジ回帰が回帰係数の2乗和に対してペナルティをかけるのに対し、LASSOでは絶対値の和に対してペナルティをかける。このことが、\(\lambda\)が大きくなったときの係数の収縮のしかたが異なることにつながっている。
関数cv.glmnetを用いて交差検証で\(\lambda\)を最適化する。
model <- cv.glmnet(x, y, alpha = 1, nfolds = length(y), grouped = F)
plot(model)
最適化された\(\lambda\)で予測値を計算する。
y.pred <- predict(model, newx = x.test, s = "lambda.min")range <- range(c(df.test$y, y.pred))
plot(df.test$y, y.pred, xlim = range, ylim = range)
abline(0, 1)
cor(df.test$y, y.pred)## lambda.min
## [1,] 0.7336718リッジ回帰に比べて、高い精度のモデルが得られた。
なお、いずれの回帰モデルも
\[ y_i = b_0 + \sum_{j=1}^J b_j x_{ij} + e_i \] という式で表され、回帰係数\(b_j\)の値だけが異なる。そこで、上述した5つのモデルの回帰係数\(b_j\)を比較してみる。
# 重回帰
model <- lm(y ~ ., data = df.train)
b.mlr <- coef(model)[-1]
# PCR
model <- pcr(y ~ ., data = df.train, validation = "LOO")
b.pcr <- coef(model, ncomp = which.min(model$validation$PRESS))[,,1]
# PLS
model <- plsr(y ~ ., data = df.train, validation = "LOO")
b.pls <- coef(model, ncomp = which.min(model$validation$PRESS))[,,1]
# ridge
model <- cv.glmnet(x = as.matrix(df.train[, -1]), y = df.train$y, alpha = 0,
nfolds = length(y), grouped = F)
b.ridge <- coef(model, s = "lambda.min")[-1,1]
# lasso
model <- cv.glmnet(x = as.matrix(df.train[, -1]), y = df.train$y, alpha = 1,
nfolds = length(y), grouped = F)
b.lasso <- coef(model, s = "lambda.min")[-1,1]b.all <- data.frame(b.mlr, b.pcr, b.pls, b.ridge, b.lasso)
pairs(b.all)
PCR, PLS, リッジ回帰で、同様の回帰係数が得られていることが分かる。一方、重回帰と他のモデルの回帰係数は大きく異る。LASSOも他のモデルと異なるが、重要な変数はある程度一致している。
比較的精度が高かったPCR, PLS, LASSOについて回帰係数を棒グラフを用いて図示する。
op <- par(mar = c(8, 4, 2, 2))
barplot(t(b.all[, c("b.pcr", "b.pls", "b.lasso")]), beside = T, las = 2, legend = c("PCR", "PLS", "LASSO"))
par(op)棒グラフを用いて回帰係数\(b_i\)を視覚的に表すことにより、嗜好性に強くかかわる変数と、そうでなない変数がある事がわかる。また、嗜好性への影響の正負もわかる。
機械学習法の利用
最後に、機械学習法を用いて予測モデルを作ってみる。
まずは、様々なデータに汎用的に用いられているrandom forestを適用する。
まずは、モデルを構築する。
model <- ranger(y ~ ., data = df.train)
model## Ranger result
##
## Call:
## ranger(y ~ ., data = df.train)
##
## Type: Regression
## Number of trees: 500
## Sample size: 89
## Number of independent variables: 57
## Mtry: 7
## Target node size: 5
## Variable importance mode: none
## Splitrule: variance
## OOB prediction error (MSE): 9.607521
## R squared (OOB): 0.2652215予測をする。
y.pred <- predict(model, data = df.test)$predictionsrange <- range(c(df.test$y, y.pred))
plot(df.test$y, y.pred, xlim = range, ylim = range)
abline(0, 1)
cor(df.test$y, y.pred)## [1] 0.634338予測精度は多変量回帰に比べて低い。しかし、機械学習では、ハイパーパラメータを最適化することが重要である。
そこで、tuneRanger関数を用いて、3つのハイパーパラメータ（mtry, min.node.size, sample.fraction）を最適化する。
task <- makeRegrTask(data = df.train, target = "y")
res <- tuneRanger(task, show.info = F)
res$recommended.pars## mtry min.node.size sample.fraction mse exec.time
## 1 21 10 0.6808397 8.831716 0.1132ハイパーパラメータを最適化されたモデルで予測する。
pred <- predict(res$model, newdata = df.test)
y.pred <- pred$data$responserange <- range(c(df.test$y, y.pred))
plot(df.test$y, y.pred, xlim = range, ylim = range)
abline(0, 1)
cor(df.test$y, y.pred)## [1] 0.5850388残念ながら、予測精度の大きな改良はみられない。
なお、random forestでも変数の重要度を評価できる。
pars <- res$recommended.pars
model <- ranger(y ~ ., data = df.train, importance = "permutation",
mtry = pars$mtry,
min.node.size = pars$min.node.size,
sample.fraction = pars$sample.fraction)
importance(model)## SSC TA pH Fructose Glucose
## 8.686109e-01 1.207896e+00 3.176977e-01 1.184989e+00 1.064287e+00
## Sucrose X590.86.3 X616.25.1 X107.87.9 X110.62.3
## 1.978088e-01 4.549886e-02 -9.638693e-03 -9.431292e-03 4.113836e-02
## X105.37.3 X623.42.7 X123.51.3 X1576.87.0 X71.41.0
## -2.627353e-02 7.588712e-02 -1.881558e-02 8.314545e-02 9.124974e-02
## X1576.95.0 X108.88.3 X556.24.1 X66.25.1 X6728.26.3
## 2.894042e-02 1.060768e-02 7.276029e-02 1.583107e-01 6.438577e-03
## X928.95.0 X111.27.3 X123.92.2 X110.43.0 X111.71.7
## 4.397677e-03 2.691044e-02 8.358215e-03 1.014712e-02 9.768453e-03
## X1191.16.8 X106.70.7 X7785.70.8 X928.68.7 X142.62.1
## -8.993501e-04 -8.031288e-02 -1.676667e-04 3.613982e-02 9.652285e-04
## X110.93.0 X111.13.7 X123.66.0 X3681.71.8 X142.92.7
## -2.688403e-02 -1.523315e-02 -4.407579e-02 -1.020546e-02 -3.979032e-02
## X2497.18.9 X104.76.7 X5989.27.5 X470.82.6 X122.78.1
## 2.274053e-01 -9.740107e-05 -5.590511e-02 1.878263e-03 1.006695e-02
## X821.55.6 X78.70.6 X124.19.6 X2639.63.6 X53398.83.7
## 8.094032e-02 1.160430e-01 9.043543e-03 3.556612e-02 6.005868e-02
## X112.40.3 X119.36.8 X106.26.3 X141.27.5 X112.12.9
## -4.942633e-02 2.932127e-02 1.040219e-02 8.251648e-03 3.762604e-01
## X629.50.5 X3879.26.3 X689.67.8 X1139.30.6 X5989.33.3
## -1.096063e-02 3.031417e-02 5.133754e-02 -2.404518e-02 2.147746e-02
## X43219.68.7 X582.16.1
## 3.125251e-02 4.827291e-02棒グラフを用いて表示してみる。
op <- par(mar = c(8, 4, 2, 2))
barplot(importance(model), las = 2)
par(op)揮発性物質の中では、
vc[vc$CAS_id == "112-12-9", ]## CAS_id Name Formula
## 44 112-12-9 2-Undecanone C11 H22 Oの重要性が高そうである。ただし、回帰係数の符号からすると、同物質は多いほど嗜好性が低くなるのかもしれない。
多変量と多変量の関係
多次元の嗜好性スコアと多次元の生化学物質間の関係を解析する。
ここでは、正準相関分析を用いる。
var.id <- c(sugaci.id, volati.id)
x <- scale(bb5[, var.id])
var.id <- liking.id
y <- scale(bb5[, var.id])res <- cancor(x, y)x.score <- x %*% res$xcoef
y.score <- y %*% res$ycoefplot(x.score, col = bb5$Genotype, pch = as.numeric(bb5$Location), asp = 1,
xlab = "CC1", ylab = "CC2", main = "Canonical components for X")
plot(y.score, col = bb5$Genotype, pch = as.numeric(bb5$Location), asp = 1,
xlab = "CC1", ylab = "CC2", main = "Canonical components for Y")
op <- par(mfrow = c(1, 2))
for(i in 1:2) {
plot(x.score[,i], y.score[,i], col = bb5$Genotype, pch = as.numeric(bb5$Location), asp = 1,
xlab = "CC_X", ylab = "CC_Y", main = paste("Canonical component", i))
}
par(op)barplot(res$cor, names.arg = paste0("CC", 1:length(res$cor)),
xlab = "Canonical Component", ylab = "Correlation coefficient")
par(mfrow = c(1, 2))
plot(res$xcoef, type = "n", xlab = "CC1", ylab = "CC2")
text(res$xcoef, rownames(res$xcoef))
abline(h = 0, v = 0)
plot(res$ycoef, type = "n", xlab = "CC1", ylab = "CC2")
text(res$ycoef, rownames(res$ycoef))
abline(h = 0, v = 0)
多変量から多変量を予測する
最後に変数群Xから変数群Yを予測するモデルを構築する。ここでは、PLS回帰を用いる。 PLS回帰では、XとYの共分散が最大になるような1次結合を求める。 なお、XとYの相関を最大にするような1次結合を求めるのが、先述した正準相関分析である。
まずは、訓練用と評価用のデータを作成する。
x.train <- x[bb5$Year != "2014",]
y.train <- y[bb5$Year != "2014",]
x.test <- x[bb5$Year == "2014",]
y.test <- y[bb5$Year == "2014",]次に、PLS回帰を行う。
n.comp <- 30
model <- plsr(y.train ~ x.train, ncomp = n.comp, validation = "LOO")
plot(model, "validation")
おおよそ、5つの潜在変数で説明できることがわかる。
n.comp <- apply(model$validation$PRESS, 1, which.min)
n.comp## Overall.Liking Texture Sweetness Sourness Flavor
## 5 4 4 13 5ここでは、酸度に対する嗜好性の予測精度は落ちるが、他の変数を考慮して5を採択する。
n.comp <- 5
y.pred <- predict(model, newdata = x.test)[,,n.comp]1変数ずつ図示して精度を確認する。
op <- par(mfrow = c(2,3))
for(i in 1:ncol(y.test)) {
range <- range(c(y.test[,i], y.pred[,i]))
plot(y.test[,i], y.pred[,i], xlim = range, ylim = range, main = colnames(y.test)[i])
abline(0, 1)
}
par(op)
予測値と観察値の相関係数を計算する。
cor(y.pred, y.test)## Overall.Liking Texture Sweetness Sourness Flavor
## Overall.Liking 0.69579345 0.1958729 0.715440409 -0.5381014 0.52986314
## Texture 0.28641555 0.5228885 0.151629293 0.2303116 0.59659109
## Sweetness 0.69238682 0.0548687 0.763353901 -0.6628817 0.45149431
## Sourness -0.51351036 0.1137000 -0.618661333 0.8561681 0.01554278
## Flavor 0.06206554 0.1978386 -0.001608782 0.4211799 0.50174012対角の数値が精度である。甘さや酸っぱさに対する予測精度は高い。 一方、食感、香りに対する予測精度は甘さや酸っぱさに比べると低い。 また、全体的な嗜好性の予測値（および観察値）が、甘さの観察値（および予測値）と強い関係を持っており、 甘さが嗜好性の重要な因子であることが分かる。
最後に変量Yのに対する回帰係数を図示する。モデルは、 \[ {\bf y}_i = {\bf b}_0 + {\bf x}_i' {\bf B}_k + {\bf e}_i \] となる。\({\bf y}_i\)はi番目の標本の嗜好性の観察値のベクトル（5次元）、\({\bf x}_i'\)はi番目の標本の生化学物質の観察値のベクトル｛58次元｝である。\({\bf B}_k\)は、回帰係数を要素とする行列｛58 x 5｝であり、各行がXの各変数に対応し、各列がYの各変数に対応する。
op <- par(mar = c(8, 4, 2, 2), mfrow = c(2, 1))
bmat <- coef(model, ncomp = n.comp)[,,1]
range <- range(bmat)
colors <- c("red", "orange", "blue", "yellow", "green")
barplot(t(bmat[1:29,]), beside = T, las = 2, col = colors, ylim = range, legend = colnames(bmat))
barplot(t(bmat[30:58,]), beside = T, las = 2, col = colors, ylim = range)
par(op)556.24.1、110.43.0、582.16.1は、いずれもパネリストの感性スコアに関連が強そうに見える。
vc[vc$CAS_id %in% c("556-24-1", "110-43-0", "582-16-1"),]## CAS_id Name Formula
## 12 556-24-1 Methyl isovalerate C6 H12 O2
## 18 110-43-0 2-Heptanone C7 H14 O
## 53 582-16-1 2,7-Dimethyl Naphthalene C12 H12ただし、予測用のモデルの回帰係数の解釈は、変数間に存在する相関関係により難しい。 解釈のためのモデリングと、予測のためのモデリングは、目指すものが異なることに十分注意する。