第6回講義：変数群間の関係のモデル化
岩田洋佳
2022/08/30
必要なパッケージ
require(tidyverse)
require(corrplot)
require(CCA)
require(CCP)
require(pls)
require(corrr)
require(igraph)
require(blockmodels)背景
近年、様々な計測技術の発展に伴い、同一サンプルに対して複数の異なる手法でデータを収集されることは珍しくない。例えば、オミクス技術の発展に伴い、野外で栽培された植物の表現型、ゲノム、遺伝子発現（トランスクリプトーム）、ミネラル（イオノーム）、代謝物（メタボローム）、さらには微生物叢（マイクロバイオーム）といった異なるオミクスデータを統合した解析も可能となってきている。現在、こうしたデータを解析するための手法がさまざまに開発されている。同じサンプルを多次元データとして計測した複数のデータセット間の関係をモデル化できれば、両データセットに共通するメカニズム、あるいは両データセットをつなぐメカニズムを解明したり、一方のデータセットから他方を予測したりことも可能になる。これまでの講義では、複数の変数を要約する方法（第3回講義）と、複数の変数と一つの変数の関係をモデル化する方法（第4回講義）を説明してきた。本講義では、複数変数の集合（変数のグループ）間の関係をモデル化する方法について説明する。具体的には、正準相関分析と多変量 PLS 回帰について説明する。
解析データ
ここでは、Gilbertら（2015, PLoS ONE 10: e0138494）のブルーベリーデータを解析対象データとして使用する。この研究は、複数のブルーベリー品種を用いた多環境（複数の試験場、複数の栽培年）試験（Multi-environmental trials: MET）において、嗜好性などの官能評価結果と、酸度、糖度、多数の揮発成分などの生化学的特性との関係を分析したものである。ここでは、これらの変数を、官能評価の変数群、糖度・酸度の変数群、揮発性成分の変数群に分け、それぞれの関係を解析する。
まず、データをロードする。
bb <- read.csv("blueberry.csv", row.names = 1, stringsAsFactors = T)データを確認する。
head(bb)[, 1:20]## 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官能評価に関する変数群、糖度・酸度に関する変数群、揮発性成分に関する変数群を抽出する。
bb.liking <- bb %>% dplyr::select(Overall.Liking, Texture, Sweetness, Sourness, Flavor)
bb.sugaci <- bb %>% dplyr::select(SSC, TA, pH, Fructose, Glucose, Sucrose)
bb.volcom <- bb %>% dplyr::select(17:67)抽出された変数について、平均0、分散1に標準化する。
bb.liking <- scale(bb.liking)
bb.sugaci <- scale(bb.sugaci)
bb.volcom <- scale(bb.volcom)変数のグループ間の関係を見てみよう。まず、官能評価に関する変数と、糖度や酸度に関する変数との相関を見てみる。
corrplot::corrplot(cor(bb.liking, bb.sugaci))
次に、官能評価に関する変数と揮発性成分に関する変数との相関について見てみる。
corrplot::corrplot(cor(bb.liking, bb.volcom))
糖度と酸度、つまり揮発性成分が官能評価と関係していることがわかる。この関係については、次節以降で分析する。
正準相関分析(CCA)
まず、正準相関分析を行う。正準相関分析は、2つの変数グループの関係を分析するためによく使われる。主成分分析では、1つの変数群（多変量）内の変動をうまくまとめられるような変数の線形結合を求める。例えば、第一主成分では、変数の線形結合から得られるスコアの分散が最大となるように線形結合の重みが決定される。一方、正準相関分析では、2つのグループそれぞれについて線形結合を求め、各グループの線形結合から得られるスコアが2つのグループ間の相関を最大にするように重みが決定される。つまり、主成分分析の目的が1つのグループの変動をまとめることであるのに対し、正準相関分析の目的は、2つのグループ間の関係をまとめる線形結合を抽出することである。
正準相関分析を理解するためには、以下の文献が参考になる。 https://www.jstatsoft.org/article/download/v023i12/212 ここで使われているRパッケージの論文である。
それでは、官能評価に関する変数群（bb.liking）と、糖度・酸度に関する変数群（bb.sugaci）の間に見出された関係を分析してみよう。まず、2つの変数グループの間の相関係数を計算する。
r <- cor(bb.liking, bb.sugaci)
r## SSC TA pH Fructose Glucose
## Overall.Liking 0.51477215 -0.40602941 0.2108530 0.5426572 0.40775563
## Texture 0.14080843 0.09484586 -0.2826011 0.0977051 0.02121231
## Sweetness 0.58506854 -0.52358535 0.3486516 0.6199120 0.48268452
## Sourness -0.09534468 0.75619822 -0.7947494 -0.2106074 -0.22223585
## Flavor 0.53716720 0.08850711 -0.2646968 0.4422369 0.27885512
## Sucrose
## Overall.Liking 0.16619396
## Texture -0.08859778
## Sweetness 0.23099296
## Sourness -0.26724321
## Flavor -0.01842639前図に示すように、Sweetnessと糖度（SSD、Fructoseなど）、Sournessと酸度（TA、pH）には強い相関があることがわかる。一方、Textureは、酸味や糖度とは強い関係がない。
それでは、CCAパッケージに含まれる関数ccを用いて、正準相関分析を行ってみる。
cc.res <- CCA::cc(bb.liking, bb.sugaci)まず、正準変量（2つのグループ間の相関が最大になるように抽出された線形結合）間の相関を見よう。
cc.res$cor## [1] 0.8633334 0.7492393 0.3569748 0.1568234 0.0316673合計5つの値が表示される。これらは順に、1番目、2番目、…、5番目の正準相関係数と呼ばれる。これらは、上述した線形結合によって得られた正準変量間の相関を表している。例えば、第1正準相関係数は、官能評価変数群の第1正準変量と、糖度・酸度変数群の第1正準変量との相関を表している。
では、線形結合の重みを確認してみよう。まず、官能評価変数群の正準変量の重み（係数）を確認する。
cc.res$xcoef## [,1] [,2] [,3] [,4] [,5]
## Overall.Liking 0.37558775 -0.6669887 0.2653372 -0.5098692 2.4725880
## Texture -0.01751421 0.4921982 0.8261774 1.0995902 -0.5928885
## Sweetness -0.64229527 -1.0921066 2.2192481 -1.1718423 -2.4784086
## Sourness 0.78028546 -1.2505500 1.7946523 -1.4371707 -0.3648066
## Flavor 0.06077118 0.2895173 -2.8439756 1.2043985 0.4137175例えば、最初の正準変量については、Overall.liking、Sweetness、Sournessへの係数が大きく、Sweetnessはマイナス、Sournessはプラスのスコアになっていることがわかる。
次に、糖度、酸味の変数群の正準変量に対する係数を確認してみる。
cc.res$ycoef## [,1] [,2] [,3] [,4] [,5]
## SSC -0.09327369 -0.46280816 -1.1923756 0.7729025 -1.0027797
## TA 0.49752954 0.39227248 -1.1337631 -0.8631082 0.2851472
## pH -0.47261924 0.71957396 -1.2110286 -0.5140783 0.3332396
## Fructose -0.05136096 -0.61917357 0.5367650 -0.5297348 1.5742211
## Glucose -0.07895030 0.21512869 0.3562663 -0.6039259 -0.3984859
## Sucrose -0.15873722 -0.03217157 0.1281382 -0.2763046 -0.7128678例えば、最初の正準変数では、酸性度に関係する2つの変数TAとpHに対する係数が大きく、TAは正、pHは負である。
こうして得られた線形結合を用いて、正準スコアを算出する。まず、官能評価変数群の正準スコアを確認する。
head(cc.res$scores$xscores)## [,1] [,2] [,3] [,4] [,5]
## 2012 C1 Emerald -0.4787260 1.06445333 -0.85168474 1.01248085 1.14016543
## 2012 C1 Endura 1.4462456 0.43651453 -1.96914892 0.77525507 -0.41150205
## 2012 C1 Farthing -0.4371230 0.33067509 -1.44325673 0.55560692 0.99855366
## 2012 C1 Meadowlark -0.3726521 1.88890992 0.20151694 2.26847620 0.32682781
## 2012 C1 Primadonna -1.0289363 -0.01410947 -0.08112546 -0.07352395 -0.42155439
## 2012 C1 Scintilla -0.2561780 -0.83410523 1.66838970 -1.30293650 -0.02699176このスコアは、正準変数の係数を用いて算出することができる。
myscores <- bb.liking %*% cc.res$xcoef
head(myscores)## [,1] [,2] [,3] [,4] [,5]
## 2012 C1 Emerald -0.4787260 1.06445333 -0.85168474 1.01248085 1.14016543
## 2012 C1 Endura 1.4462456 0.43651453 -1.96914892 0.77525507 -0.41150205
## 2012 C1 Farthing -0.4371230 0.33067509 -1.44325673 0.55560692 0.99855366
## 2012 C1 Meadowlark -0.3726521 1.88890992 0.20151694 2.26847620 0.32682781
## 2012 C1 Primadonna -1.0289363 -0.01410947 -0.08112546 -0.07352395 -0.42155439
## 2012 C1 Scintilla -0.2561780 -0.83410523 1.66838970 -1.30293650 -0.02699176同様に、酸味と糖分の変数グループの正準スコアを確認してみよう。
head(cc.res$scores$yscores)## [,1] [,2] [,3] [,4] [,5]
## 2012 C1 Emerald 0.4866830 0.31997491 0.05530756 -0.35797933 1.2349727
## 2012 C1 Endura 0.8649419 0.32346688 -0.57731698 -0.33853779 1.1371146
## 2012 C1 Farthing -0.2309199 -0.03476054 0.78433524 0.06185598 1.3019831
## 2012 C1 Meadowlark -0.4099384 1.53712939 0.76787713 0.02675023 1.2792454
## 2012 C1 Primadonna -1.1253494 0.07829178 -0.45846697 -0.53843707 1.8820045
## 2012 C1 Scintilla -0.1761669 -1.69961723 -0.64500092 -0.66245531 0.30283432つのグループ間の正準スコアの相関を計算する。
rho <- cor(cc.res$scores$xscores, cc.res$scores$yscores)
corrplot::corrplot(rho)
この図から、正準変数間の相関係数は、同じ次数の正準変数同士を除いて0であり、変量の次数が増えるに従って係数が最大値から最小値に増加することがよくわかる。同じ次数の正準変数間の相関は相関行列の対角成分であるため、次のように抽出することができる。
diag(rho)## [1] 0.8633334 0.7492393 0.3569748 0.1568234 0.0316673この値は、
cc.res$cor## [1] 0.8633334 0.7492393 0.3569748 0.1568234 0.0316673として既に確認されていたものになる。
散布図を描いて確認してみよう。
op <- par(mfrow = c(2,3))
for(i in 1:5) {
plot(cc.res$scores$xscores[,i], cc.res$scores$yscores[,i],
main = paste0("CC", i), xlab = "CS.X", ylab = "CS.Y")
abline(0,1)
}
par(op)
CCPパッケージのp.asym関数は、正準相関係数に基づいて有意な変量数を決定するために使用できる。
rho <- cc.res$cor
n <- nrow(bb.liking) # number of samples
p <- ncol(bb.liking)
q <- ncol(bb.sugaci)
CCP::p.asym(rho, n, p, q, tstat = "Wilks")## Wilks' Lambda, using F-approximation (Rao's F):
## stat approx df1 df2 p.value
## 1 to 5: 0.09497543 15.01188527 30 562.0000 0.00000000
## 2 to 5: 0.37295663 8.11428763 20 468.5940 0.00000000
## 3 to 5: 0.85025586 1.98118376 12 375.9882 0.02484448
## 4 to 5: 0.97442826 0.62140339 6 286.0000 0.71314027
## 5 to 5: 0.99899718 0.07227538 2 144.0000 0.93030841この検定の結果は、上位3つの変数が有意であることを示しています。
ここで、変数のグループとcanonical scoreの間の相関係数を計算して、変数とcanonical variatesの関係を見てみよう。相関係数は、最初の変数群をX、2番目の変数群をYとすると、以下の4つの組み合わせで計算できる。
- Xの変数（ここでは官能評価）と、Xの正準変数のスコア
- Yの変数（ここでは糖度と酸度）とYの正準変数のスコア
- Xの変数とYの正準変量のスコア
- Yの変数とYの正準変量のスコア
同じ順番の正準変量間の相関が高ければ、1と3、2と4は似たような値をとると予想される。
元の変数と正準変数の相関係数は、compute関数で計算することができる。
op <- par(mfrow = c(2,2))
ccl <- CCA::comput(bb.liking, bb.sugaci, cc.res)
corrplot::corrplot(ccl$corr.X.xscores, title = "Type 1", mar = c(0, 0, 2, 0)) # type 1
corrplot::corrplot(ccl$corr.Y.yscores, title = "Type 2", mar = c(0, 0, 2, 0)) # type 2
corrplot::corrplot(ccl$corr.X.yscores, title = "Type 3", mar = c(0, 0, 2, 0)) # type 3
corrplot::corrplot(ccl$corr.Y.xscores, title = "Type 4", mar = c(0, 0, 2, 0)) # type 4
par(op)タイプ1と3、タイプ2と4は、第1と第2の正準変数が非常に似た値であることに注意する。
また、plt.cc関数を使えば、これらを同時にチェックすることができる。まず、第1、第2の正準変量を確認してみる。
plt.cc(cc.res, var.label = T,
ind.names = paste0(substr(bb$Genotype, 1, 2), substr(bb$Location, 1, 1)))
なお、左の図で使われている相関係数は、タイプ1とタイプ4（つまり、Xの正準変数に対する2つの変数群からの相関）である。
正則化正準相関分析(RCCA)
正準相関分析は、同一グループ内に強い相関を持つ変数がある場合や、変数の数が多い場合、安定した結果が得られない。変数の数がサンプルの数より多い場合もあり、その場合は通常の正準相関分析が適用できない。そのような場合は、正則化正準相関分析（RCCA）を行うことができる。ここでは、変数数が多いvolative変数群とsensory変数群との間の正則化正準相関分析を行うことにする。
正規化パラメータは estim.regul 関数を用いて、1個抜き交差検証により \(\lambda_X, \lambda_Y\) を求めることができる。これらのパラメータは、1個抜き交差検証で変数間の相関が最大になるように選択される。しかし、サンプル数が多い場合、この計算は非常に時間がかかるため、ここでは行わない。
#res.regul <- CCA::estim.regul(bb.liking, bb.volcom, plt = TRUE,
# grid1 = seq(0, 0.2, l=5),
# grid2 = seq(0, 0.2, l=5))上記のコマンドを実行すると、以下のような結果が得られる。詳細な最適化はしていませんが、正則化パラメータはどちらもおおよそ0.05と0.05の値をとればよさそうである。 ! CV_res
では、正則化パラメータの値を用いて、正則化正準相関分析を行ってみよう。
rcc.res <- CCA::rcc(bb.liking, bb.volcom, 0.05, 0.05)
rcc.res$cor## [1] 0.8467243 0.7607413 0.7082493 0.5117604 0.3628217p.asym関数を使ったテストで、有効な正準変数を確認してみよう。
rho <- rcc.res$cor
n <- nrow(bb.liking) # number of samples
p <- ncol(bb.liking)
q <- ncol(bb.volcom)
CCP::p.asym(rho, n, p, q, tstat = "Wilks")## Wilks' Lambda, using F-approximation (Rao's F):
## stat approx df1 df2 p.value
## 1 to 5: 0.0380906 1.7405054 255 478.6591 1.132791e-07
## 2 to 5: 0.1345682 1.2598567 200 385.9101 2.824223e-02
## 3 to 5: 0.3194326 0.9194113 147 291.6630 7.150352e-01
## 4 to 5: 0.6409379 0.5085486 96 196.0000 9.998645e-01
## 5 to 5: 0.8683604 0.3193182 47 99.0000 9.999829e-01最初の2つの変量は有意である。
次のコマンドは、正準スコアと2つの変数のグループの間の相関をプロットする。また、その時のサンプル全体のcanonical scoreをプロットする。
CCA::plt.cc(rcc.res, var.label = T,
ind.names = paste0(substr(bb$Genotype, 1, 2), substr(bb$Location, 1, 1)))
文字が小さく少し確認しにくいが、X1575.95.0、X821.55.6、X470.82.6が大きな正の相関を、X556.24.1が負の相関を持っていることが分かる。
官能評価についても、抽出された正準変量が、官能評価と糖度・酸度との正準相関分析におけるものと非常に類似していることがわかる。
corrplot関数により、揮発性物質と正準変量との関係を詳細に確認することもできる。
corrplot::corrplot(cor(rcc.res$scores$xscores, bb.volcom))
予想通り、X1575.95.0, X821.55.6, X470.82.6
は大きな正の相関を持ち、X556.24.1 は負の相関を持つ。
なお、正準変量は以下のように取り出して可視化できる。
plot(rcc.res$scores$xscores[,1:2],
col = bb$Genotype, pch = as.numeric(bb$Location),
xlab = "Canonical variate 1",
ylab = "Canonical variate 2",
main = "Canonical scores based on X variables")
legend("topright",
legend = levels(bb$Genotype),
col = 1:nlevels(bb$Genotype),
pch = 1)
legend("bottomright",
legend = levels(bb$Location),
pch = 1:nlevels(bb$Location),
col = 1)
また、主成分スコアのように、さらなる分析に利用することも可能です。
y <- rcc.res$scores$xscores[,1]
variety <- bb$Genotype
location <- bb$Location
year <- bb$Year
res <- lm(y ~ variety * location * year)
anova(res)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## variety 5 49.934 9.9869 27.3292 < 2.2e-16 ***
## location 2 25.836 12.9182 35.3509 9.794e-13 ***
## year 1 0.025 0.0255 0.0698 0.79212
## variety:location 9 7.693 0.8548 2.3392 0.01837 *
## variety:year 5 1.340 0.2680 0.7334 0.59979
## location:year 2 0.673 0.3367 0.9215 0.40079
## variety:location:year 9 16.313 1.8125 4.9600 1.198e-05 ***
## Residuals 117 42.755 0.3654
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1第1正準スコアは、品種、場所、品種×場所×年次の交互作用の効果が非常に有意であり、年次の主効果は有意でないことを示している。
また、変数群Yについて得られた正準スコアに基づいても、同様のグラフが描ける。
plot(rcc.res$scores$yscores[,1:2],
col = bb$Genotype, pch = as.numeric(bb$Location),
xlab = "Canonical variate 1",
ylab = "Canonical variate 2",
main = "Canonical scores based on Y variables")
legend("topright",
legend = levels(bb$Genotype),
col = 1:nlevels(bb$Genotype),
pch = 1)
legend("bottomright",
legend = levels(bb$Location),
pch = 1:nlevels(bb$Location),
col = 1)
スコアは似ているが、まったく同じではない。
X変数の正準スコアで行われたように、Y変数の正準スコアでANOVAを実行する。
y <- rcc.res$scores$yscores[,1]
variety <- bb$Genotype
location <- bb$Location
year <- bb$Year
res <- lm(y ~ variety * location * year)
anova(res)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## variety 5 59.865 11.9730 35.6617 < 2.2e-16 ***
## location 2 20.231 10.1153 30.1284 2.790e-11 ***
## year 1 0.077 0.0774 0.2306 0.6320
## variety:location 9 5.796 0.6440 1.9181 0.0558 .
## variety:year 5 1.553 0.3106 0.9251 0.4674
## location:year 2 0.627 0.3134 0.9333 0.3961
## variety:location:year 9 13.661 1.5179 4.5211 4.079e-05 ***
## Residuals 117 39.282 0.3357
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1結果はまたしても似ているが、まったく同じではない。
多変量PLS回帰
第2回目の講義では、説明変数と応答変数の共分散を最大化する潜在変数に基づく重回帰を行う方法であるPLS回帰を紹介した。PLS回帰は、応答変数が多変量である場合にも使用することができる。その場合，PLS 回帰は正準相関分析に似たアプローチになる．違いは，正準相関分析では，変数グループXとYの間の相関を最大化するように潜在変数が抽出される．PLS 回帰では，変数グループ（被説明変数Xと説明変数Y）間の共分散を最大化するように潜在変数が抽出される．そして、被説明変数Xの潜在変数を説明変数Yの潜在変数に回帰させ、これに基づいてXとYの回帰モデルを作成する。
詳しくは https://ieeexplore.ieee.org/document/6553737 を参考にするとよい。
さて、Yが多次元であることを除けば、レクチャー02で説明したのと同じように分析することができる。実行してみよう。
pls.res <- pls::plsr(bb.liking ~ bb.volcom, ncomp = 10)結果を確認する。
summary(pls.res)## Data: X dimension: 151 51
## Y dimension: 151 5
## Fit method: kernelpls
## Number of components considered: 10
## TRAINING: % variance explained
## 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps 7 comps
## X 12.28 22.51 29.50 38.44 45.46 51.26 54.89
## Overall.Liking 30.82 42.98 43.42 47.54 56.26 56.76 57.44
## Texture 21.29 21.98 43.10 43.96 44.30 45.99 53.57
## Sweetness 19.71 34.81 34.92 43.20 52.55 53.46 54.07
## Sourness 0.10 40.31 54.34 62.94 62.97 67.12 67.16
## Flavor 21.46 23.30 30.31 30.91 38.35 47.82 49.18
## 8 comps 9 comps 10 comps
## X 60.48 63.70 66.34
## Overall.Liking 57.45 58.98 60.99
## Texture 55.31 55.35 55.71
## Sweetness 55.12 58.28 60.71
## Sourness 67.96 71.80 72.83
## Flavor 51.77 51.79 52.68交差検証を行って、潜在変数の数を決定する。
pls.res <- pls::plsr(bb.liking ~ bb.volcom, ncomp = 10, validation = "LOO")
plot(RMSEP(pls.res))
潜在変数数が6の場合について、交差検証での予測値と実際の観測値との関係を確認する。
plot(pls.res, ncomp = 6, line = T)
すべての標本について、潜在変数の値をプロットする。散布点を品種によって色分けする。
plot(pls.res, plottype = "scores", comps = 1:6, col = bb$Genotype)
上記で算出した正準変量とPLS回帰の潜在変数とのスコアの相関係数を算出し、可視化する。
r <- cor(rcc.res$scores$xscores, pls.res$scores)
corrplot::corrplot(r)
1番目の正準変量がPLSの2番目の潜在変数に強く関係し、2番目の正準変量がPLSの1番目の潜在変数に強く関係していることがわかる。
潜在変数のスコア計算で使用される線形結合の重みをプロットする。
plot(pls.res, plottype = "loadings", comps = 1:6, lty = 1)
この重みに基づいて標本の得点が計算される。
第1潜在変数と第2潜在変数の重みと、それに基づいて計算されたサンプルのスコアをプロットして表示する。
plot(pls.res, plottype = "biplot", comps = 1:2)
グラフは主成分分析のバイプロットに似ている。
得られた潜在変数をもとに、回帰を行う。潜在変数の数が2、4、6のときに得られた回帰モデルをもとに、説明変数（この場合は揮発性物質）の回帰係数をプロットする。
plot(pls.res, "coef", ncomp = c(2,4,6), legendpos = "bottomright", lty = 1)
潜在変数の数が変われば、当然説明変数の回帰係数も変わる。前節の重みとは異なり、変数群Yの各変数に対して係数が計算されていることに注意しよう。
最後に、潜在変数が 6 個の場合の説明変数（ここでは volatile matter）の回帰係数を確認しよう。
corrplot::corrplot(t(pls.res$coefficients[,,6]), is.corr = F)
Sournessにみられる関係を見ると、各揮発成分と第一の正準変数の関係が似ていることがわかる。
ネットワーク解析
多数の要因間の関係を確認したり、理解したり、あるいは、モデル化するために、ネットワーク解析がよく用いられる。ここでは、相関計数に基づくネットワークの構築と、ネットワーク構造をもとにした要因のグルーピングを行ってみる。
まずは、官能評価変数群と糖度・酸度変数群をまとめてネットワークを描いてみよう。corrrパッケージのcorrelate関数を用いて、要因間の相関係数を計算する。
cor.ls <- corrr::correlate(data.frame(bb.liking, bb.sugaci))## Correlation computed with
## • Method: 'pearson'
## • Missing treated using: 'pairwise.complete.obs'少数点以下、2桁で図示する。
cor.ls %>% corrr::fashion()## term Overall.Liking Texture Sweetness Sourness Flavor SSC TA
## 1 Overall.Liking .52 .89 -.43 .62 .51 -.41
## 2 Texture .52 .35 .27 .67 .14 .09
## 3 Sweetness .89 .35 -.52 .63 .59 -.52
## 4 Sourness -.43 .27 -.52 .27 -.10 .76
## 5 Flavor .62 .67 .63 .27 .54 .09
## 6 SSC .51 .14 .59 -.10 .54 -.08
## 7 TA -.41 .09 -.52 .76 .09 -.08
## 8 pH .21 -.28 .35 -.79 -.26 .01 -.77
## 9 Fructose .54 .10 .62 -.21 .44 .82 -.17
## 10 Glucose .41 .02 .48 -.22 .28 .68 -.10
## 11 Sucrose .17 -.09 .23 -.27 -.02 .17 -.09
## pH Fructose Glucose Sucrose
## 1 .21 .54 .41 .17
## 2 -.28 .10 .02 -.09
## 3 .35 .62 .48 .23
## 4 -.79 -.21 -.22 -.27
## 5 -.26 .44 .28 -.02
## 6 .01 .82 .68 .17
## 7 -.77 -.17 -.10 -.09
## 8 .14 .11 .15
## 9 .14 .81 .28
## 10 .11 .81 .47
## 11 .15 .28 .47相関係数に基づくネットワークを描く。ネットワークにおける要因（ノード）間の繋がり（エッジ）の有無は、min_corというオプションで指定できる。min_corで指定された値よりも「絶対値が」小さいものは繋がりがない（ノード間のエッジがない）とされる。
corrr::network_plot(cor.ls, min_cor = 0.3)
corrr::network_plot(cor.ls, min_cor = 0.5)
相関行列から、ネットワークの要因間のつながり（ノード間の連結、エッジの有無）を表す隣接行列（adjucency matrix）を得る。エッジの有無を決めるしきい値は、ここでは、0.5とする。
rmat <- cor.ls %>% dplyr::select(-term)
amat <- matrix(0, nrow(rmat), ncol(rmat))
rownames(amat) <- colnames(amat) <- colnames(rmat)
amat[abs(rmat) > 0.5] <- 1
head(amat)## Overall.Liking Texture Sweetness Sourness Flavor SSC TA pH
## Overall.Liking 0 1 1 0 1 1 0 0
## Texture 1 0 0 0 1 0 0 0
## Sweetness 1 0 0 1 1 1 1 0
## Sourness 0 0 1 0 0 0 1 1
## Flavor 1 1 1 0 0 1 0 0
## SSC 1 0 1 0 1 0 0 0
## Fructose Glucose Sucrose
## Overall.Liking 1 0 0
## Texture 0 0 0
## Sweetness 1 0 0
## Sourness 0 0 0
## Flavor 0 0 0
## SSC 1 1 0隣接行列から、逆にネットワーク（グラフ構造）を描いてみる。描画には、igraphパッケージを用いる。
g <- igraph::graph_from_adjacency_matrix(amat)
plot(g, vertex.shape="none", vertex.label.color="black")
では、全てのデータ（官能評価、糖度・酸度、揮発性物質）を含んだ相関行列を計算し、それをもとに同じ解析を行ってみよう。
cor.all <- corrr::correlate(data.frame(bb.liking, bb.sugaci, bb.volcom))## Correlation computed with
## • Method: 'pearson'
## • Missing treated using: 'pairwise.complete.obs'ネットワークを図示する。
corrr::network_plot(cor.all, min_cor = 0.3)## Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
corrr::network_plot(cor.all, min_cor = 0.5)## Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
連結行列を計算する。ここでは、連結の有無を決めるためのしきい値を0.4とする。
rmat <- cor.all %>% dplyr::select(- term)
amat <- matrix(0, nrow(rmat), ncol(rmat))
rownames(amat) <- colnames(amat) <- colnames(rmat)
colnames(bb.liking)## [1] "Overall.Liking" "Texture" "Sweetness" "Sourness"
## [5] "Flavor"amat[abs(rmat) > 0.4] <- 1
head(amat)[,1:10]## Overall.Liking Texture Sweetness Sourness Flavor SSC TA pH
## Overall.Liking 0 1 1 1 1 1 1 0
## Texture 1 0 0 0 1 0 0 0
## Sweetness 1 0 0 1 1 1 1 0
## Sourness 1 0 1 0 0 0 1 1
## Flavor 1 1 1 0 0 1 0 0
## SSC 1 0 1 0 1 0 0 0
## Fructose Glucose
## Overall.Liking 1 1
## Texture 0 0
## Sweetness 1 1
## Sourness 0 0
## Flavor 1 0
## SSC 1 1igraphを用いて隣接関係を表すグラフを描く。
g <- igraph::graph_from_adjacency_matrix(amat)
plot(g, vertex.shape="none", vertex.label.color="black")
こうして得られた隣接行列をもとに、要因のグルーピングを行うことができる。ここでは、ネットワーク構造解析に用いられる確率的ブロックモデル（stochastic block model:SBM）を用いる。SBMは、例えば、人間の社会的なつながりから潜在的なグループを推定したり、生態系内での生物間の関連（生態系サービスの有無）をもとに種の属する潜在的なグループを推定したりするのに用いられる。
model <- blockmodels::BM_bernoulli("SBM_sym", amat, verbosity = 0, plotting = "")
model$estimate()ネットワーク構造を説明するのに最適な潜在的なグループ（ブロック）数を確認する。
Q <- which.max(model$ICL)
Q## [1] 3モデルの推定の中では、各要因の属するグループ（ブロック）の推定も行われている。
Z <- model$memberships[[Q]]$Z
head(Z)## [,1] [,2] [,3]
## [1,] 0.001610306 0.9967794 0.001610306
## [2,] 0.001610306 0.9967794 0.001610306
## [3,] 0.001610306 0.9967794 0.001610306
## [4,] 0.001610306 0.9967794 0.001610306
## [5,] 0.001610306 0.9967794 0.001610306
## [6,] 0.001610306 0.9967794 0.001610306行は要因、列がグループ（ブロック）に対応している。
colnames(Z) <- paste("block", 1:3, sep = "")
rownames(Z) <- colnames(amat)
head(Z)## block1 block2 block3
## Overall.Liking 0.001610306 0.9967794 0.001610306
## Texture 0.001610306 0.9967794 0.001610306
## Sweetness 0.001610306 0.9967794 0.001610306
## Sourness 0.001610306 0.9967794 0.001610306
## Flavor 0.001610306 0.9967794 0.001610306
## SSC 0.001610306 0.9967794 0.001610306要因をグルーピングするには、もっとも大きな値が振られているところに割り振ればよい。
gid <- apply(Z, 1, which.max)
names(gid) <- rownames(Z)
gid## Overall.Liking Texture Sweetness Sourness Flavor
## 2 2 2 2 2
## SSC TA pH Fructose Glucose
## 2 2 2 2 2
## Sucrose X590.86.3 X616.25.1 X107.87.9 X110.62.3
## 3 3 1 3 1
## X105.37.3 X623.42.7 X123.51.3 X1576.87.0 X71.41.0
## 3 3 3 1 1
## X1576.95.0 X108.88.3 X556.24.1 X66.25.1 X6728.26.3
## 1 1 2 1 1
## X928.95.0 X111.27.3 X123.92.2 X110.43.0 X111.71.7
## 1 1 3 3 1
## X1191.16.8 X106.70.7 X7785.70.8 X928.68.7 X142.62.1
## 3 3 3 3 1
## X110.93.0 X111.13.7 X123.66.0 X3681.71.8 X142.92.7
## 3 3 3 1 1
## X2497.18.9 X104.76.7 X5989.27.5 X470.82.6 X122.78.1
## 1 3 3 1 3
## X821.55.6 X78.70.6 X124.19.6 X2639.63.6 X53398.83.7
## 3 3 3 3 1
## X112.40.3 X119.36.8 X106.26.3 X141.27.5 X112.12.9
## 3 3 3 3 3
## X629.50.5 X3879.26.3 X689.67.8 X1139.30.6 X5989.33.3
## 3 3 3 3 3
## X43219.68.7 X582.16.1
## 3 3今回は、官能評価と糖度・酸度で概ね１つのブロックが得られ、揮発性物質で別のブロックが構成されている。ただし、揮発性物質にも官能評価や糖度・酸度と同じブロックに割り振られるものがあることに注意する。