必要なパッケージ
require(here)
require(corrplot)
require(CCA)
require(CCP)背景
近年では、さまざまな計測技術の発展により、同一サンプルに対して複数の異なる手法でデータを収集することが珍しくなくなってきています。例えば、オミクス技術の進展により、野外で栽培された植物について、表現型、ゲノム、遺伝子発現（トランスクリプトーム）、ミネラル（イオノーム）、代謝物（メタボローム）、さらには微生物叢（マイクロバイオーム）といった異なるオミクスデータを統合した解析が可能となっています。現在、このようなデータを解析するためのさまざまな手法が開発されています。
同じサンプルに対して多次元データとして計測された複数のデータセット間の関係をモデル化することができれば、両データセットに共通するメカニズムや、両データセットをつなぐメカニズムを解明することができるようになります。また、一方のデータセットから他方を予測することも可能になります。
これまでの講義では、複数の変数を要約する方法（第3回講義）や、複数の変数と一つの変数の関係をモデル化する方法（第4回講義）について説明してきました。本講義では、複数の変数集合（変数のグループ）間の関係をモデル化する方法について説明します。具体的には、正準相関分析と多変量PLS回帰について取り上げます。
解析データ
ここでは、Gilbertら（2015年, PLoS ONE 10: e0138494）のブルーベリーデータを解析対象として使用します。今回の解析では、これらの変数を官能評価の変数群、糖度・酸度の変数群、そして揮発性成分の変数群に分け、それぞれの関係を解析していきます。
まず、データをロードします。
bb <- read.csv(here("data", "blueberry.csv"), stringsAsFactors = T)データを確認します。
head(bb)[, 1:20]## 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
## 1 0.13 98.26 1.05
## 2 0.19 94.47 1.40
## 3 0.18 103.15 0.82
## 4 0.12 87.79 0.85
## 5 0.21 113.39 1.82
## 6 1.76 126.62 1.76官能評価に関する変数群、糖度・酸度に関する変数群、揮発性成分に関する変数群を抽出します。
attach(bb)
bb.liking <- data.frame(Overall.Liking, Texture, Sweetness, Sourness, Flavor)
bb.sugaci <- data.frame(SSC, TA, pH, Fructose, Glucose, Sucrose)
bb.volcom <- bb[,20:ncol(bb)]
detach(bb)抽出された変数について、平均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つのグループ間の関係を要約するための線形結合を抽出することにあります。このように、正準相関分析は、複数の変数を持つ2つのデータセット間の関連性を探るための強力な手法です。
それでは、官能評価に関する変数群（bb.liking）と、糖度・酸度に関する変数群（bb.sugaci）の間に見出された関係を分析してみましょう。まず、2つの変数グループの間の相関係数を計算します。
r <- cor(bb.liking, bb.sugaci)
r## SSC TA pH Fructose Glucose
## Overall.Liking 0.45787543 -0.42996614 0.2673230 0.4532099 0.365503820
## Texture 0.16088639 0.07924312 -0.2609702 0.1018886 0.000214648
## Sweetness 0.54773193 -0.54330888 0.3927478 0.5454639 0.440375612
## Sourness -0.08716536 0.71489370 -0.7632588 -0.1694723 -0.211233829
## Flavor 0.52783170 0.08369449 -0.2615606 0.4125374 0.243228580
## Sucrose
## Overall.Liking 0.17225758
## Texture -0.03432326
## Sweetness 0.22450654
## Sourness -0.22705848
## Flavor 0.02695600前図に示すように、Sweetnessと糖度（SSD、Fructoseなど）、Sournessと酸度（TA、pH）には強い相関があることがわかります。一方、Textureは、酸味や糖度とは強い関係がありません。
それでは、CCAパッケージに含まれる関数ccを用いて、正準相関分析を行ってみます。
cc.res <- CCA::cc(bb.liking, bb.sugaci)まず、正準変量（2つのグループ間の相関が最大になるように抽出された線形結合）間の相関を見ましょう。
cc.res$cor## [1] 0.83888450 0.71212691 0.31000515 0.19194542 0.03741143合計5つの値が表示される。これらは順に、1番目、2番目、…、5番目の正準相関係数と呼ばれる。これらは、上述した線形結合によって得られた正準変量間の相関を表しています。例えば、第1正準相関係数は、官能評価変数群の第1正準変量と、糖度・酸度変数群の第1正準変量との相関を表しています。
では、線形結合の重みを確認してみましょう。まず、官能評価変数群の正準変量の重み（係数）を確認します。
cc.res$xcoef## [,1] [,2] [,3] [,4] [,5]
## Overall.Liking 0.7209142 -0.2832845 0.4088056 -0.8860023 -2.3362418
## Texture -0.1467923 0.4017240 0.7838556 1.2201445 0.3925561
## Sweetness -0.3631756 -1.2102858 2.1963875 -0.7309514 2.7661220
## Sourness 1.2463888 -1.0267501 2.0114684 -1.3976928 0.7555027
## Flavor -0.2481775 -0.0127644 -2.7488954 0.8790375 -0.8459130例えば、最初の正準変量については、Overall.liking、Sweetness、Sournessへの係数が大きく、Sweetnessはマイナス、Sournessはプラスのスコアになっていることがわかります。
次に、糖度、酸味の変数群の正準変量に対する係数を確認してみます。
cc.res$ycoef## [,1] [,2] [,3] [,4] [,5]
## SSC -0.1136845 -0.627903270 -1.05764552 0.58418620 -0.9627302
## TA 0.3744957 0.478617789 -1.27566520 -0.68730888 0.4511134
## pH -0.6452807 0.619622423 -1.23019073 -0.50100270 0.2730009
## Fructose 0.1810166 -0.438035255 0.37155000 0.02555497 1.9709074
## Glucose -0.1815024 0.152493034 0.51130154 -1.29524185 -1.1728889
## Sucrose -0.1066091 -0.001958023 -0.07040604 0.37152976 0.4542422例えば、最初の正準変数では、酸性度に関係する2つの変数TAとpHに対する係数が大きく、TAは正、pHは負です。
こうして得られた線形結合を用いて、正準スコアを算出します。まず、官能評価変数群の正準スコアを確認します。
head(cc.res$scores$xscores)## [,1] [,2] [,3] [,4] [,5]
## [1,] -0.71968984 0.9840584 -0.82173198 0.7392187 -1.30192277
## [2,] 1.00540193 0.2965460 -1.83825689 0.5692594 -0.05491962
## [3,] -0.58716417 0.2250007 -1.35938013 0.2250832 -1.15759988
## [4,] -0.80998501 1.7061058 0.15213993 2.1612782 -0.72123114
## [5,] -1.03500483 -0.1772706 -0.09349533 -0.1059804 0.38145968
## [6,] 0.04584174 -0.6042149 1.67421659 -1.2355707 0.34747700このスコアは、正準変数の係数を用いて算出することができます。
myscores <- bb.liking %*% cc.res$xcoef
head(myscores)## [,1] [,2] [,3] [,4] [,5]
## [1,] -0.71968984 0.9840584 -0.82173198 0.7392187 -1.30192277
## [2,] 1.00540193 0.2965460 -1.83825689 0.5692594 -0.05491962
## [3,] -0.58716417 0.2250007 -1.35938013 0.2250832 -1.15759988
## [4,] -0.80998501 1.7061058 0.15213993 2.1612782 -0.72123114
## [5,] -1.03500483 -0.1772706 -0.09349533 -0.1059804 0.38145968
## [6,] 0.04584174 -0.6042149 1.67421659 -1.2355707 0.34747700同様に、酸味と糖分の変数グループの正準スコアを確認してみましょう。
head(cc.res$scores$yscores)## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.42088277 0.462718039 0.1286006 -0.8230133 0.22118241
## [2,] 0.78190097 0.516706496 -0.6133432 -0.5073797 0.51440606
## [3,] -0.19179031 -0.001662266 0.9964985 -0.5986725 0.06189234
## [4,] -0.62527635 1.556670219 0.9081019 -0.5027458 0.06165177
## [5,] -1.15515313 -0.063643179 -0.1659345 -1.3651238 0.38199778
## [6,] 0.05656579 -1.723777136 -0.5101177 -1.0725300 -0.005069052つのグループ間の正準スコアの相関を計算します。
rho <- cor(cc.res$scores$xscores, cc.res$scores$yscores)
corrplot::corrplot(rho)
この図から、正準変数間の相関係数は、同じ次数の正準変数同士を除いて0であり、変量の次数が増えるに従って係数が最大値から最小値に増加することがよくわかります。同じ次数の正準変数間の相関は相関行列の対角成分であるため、次のように抽出することができます。
diag(rho)## [1] 0.83888450 0.71212691 0.31000515 0.19194542 0.03741143この値は、
cc.res$cor## [1] 0.83888450 0.71212691 0.31000515 0.19194542 0.03741143として既に確認されていたものになります。
散布図を描いて確認してみましょう。
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.1269511 13.8210057 30 614.0000 0.00000000
## 2 to 5: 0.4284939 7.4489275 20 511.7101 0.00000000
## 3 to 5: 0.8693760 1.8580740 12 410.3830 0.03780229
## 4 to 5: 0.9618089 1.0223468 6 312.0000 0.41061274
## 5 to 5: 0.9986004 0.1100238 2 157.0000 0.89588184この検定の結果は、上位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.8129140 0.7395093 0.6819648 0.5252012 0.3504838p.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.05222106 1.6751521 260 538.5161 3.471140e-07
## 2 to 5: 0.15396686 1.2718012 204 433.8445 2.055055e-02
## 3 to 5: 0.33978817 0.9469252 150 327.6431 6.450875e-01
## 4 to 5: 0.63520828 0.5717887 98 220.0000 9.990309e-01
## 5 to 5: 0.87716112 0.3238458 48 111.0000 9.999851e-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について得られた正準スコアに基づいても、同様のグラフが描けます。
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 18 79.453 4.4141 15.3885 < 2.2e-16 ***
## location 2 15.639 7.8196 27.2611 1.910e-10 ***
## year 1 0.558 0.5579 1.9449 0.16578
## variety:location 9 8.296 0.9218 3.2135 0.00163 **
## variety:year 5 0.946 0.1892 0.6597 0.65475
## location:year 2 1.313 0.6564 2.2884 0.10595
## variety:location:year 9 13.017 1.4464 5.0424 9.527e-06 ***
## Residuals 117 33.561 0.2868
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1結果はまたしても似ているが、まったく同じではありません。