第4回講義：多変量に基づく一連続変量の予測
岩田洋佳
2022/08/30
必要なパッケージ
require(pls) # required to use pcr, plsr functions
require(glmnet)
require(ranger)背景
科学的な研究では、測定・制御が容易な変数を用いて、別の変数の挙動を制御・予測する必要がある場合が多々ある。例えば、肥料や灌漑、日射などの生育条件によって植物の成長を制御する場合がそうである。前者の変数（要因）から後者の変数（応答）を予測・制御するためには、一般的に回帰分析という統計手法が用いられる。また、以下の条件を満たす場合は、重回帰分析（MLR）も適切な方法である。
1.要因数が少ない 2.要因間に有意な冗長性がない。 3.要因と応答の関係がよく分かっている。
しかし、これら3つの条件のうち1つでも満たされていない場合、MLRモデルは不適切で使い物にならないことがある。例えば、要因の数が多く、要因と応答の関係がよく分かっておらず、単に予測精度の高いモデルを構築することが目的の場合である。例えば、次のような例が考えられる。
- 化学者がスペクトルデータから化学物質の化学組成を推定しようとする場合。この場合、因子は数百にも及ぶスペクトルの各波長における測定値であり、隣接する波長間には強い共線性が存在する。
- 育種家が、ゲノム中に分布する多数のDNA多型のデータから、動物や植物の形質の表現型を予測したい場合。これは、ゲノミック予測と呼ばれる方法である（Meuwissen et al. 2001） この方法を用いると、DNA多型から動物や植物の遺伝的能力を予測することが可能となり（このような選抜方法をゲノミック選抜と呼ぶ）、形質評価に要する時間の短縮や育種の迅速化が期待される。
- 農学者が、気温や日照時間などの気候条件の時間的変化に対する植物の応答をモデル化したい場合。
本講演では、主成分回帰、部分最小二乗（PLS）回帰、リッジ回帰、LASSOなどの多変量回帰分析について紹介する。これらは、要因の数が多く、要因間に強い共線性がある場合に、予測モデルを構築するための手法である。ここでは、モデルを構築することに重点を置いている。モデルを構築する際に重視するのは「応答を予測すること」であり、必ずしも変数間の潜在的な関係を解明することではないことに注意する。たとえば，主成分回帰，PLS，リッジ回帰は，応答に大きな影響を与えない因子をスクリーニングするのには適していない．しかし，予測が最終目標で，測定する因子の数を制限する現実的な必要性がない場合は，これらの方法は非常に有用なツールとなり得る．ところで、要因の数がサンプル数よりはるかに大きい場合、MLRモデルを使用することが困難になる。これは “large p small n”問題と呼ばれるものである。一方、上述した4つの手法はいずれもこの問題に対処することができる。上記のようなゲノム予測では、ゲノム全体のDNA多型（因子）の数がサンプル数を超える場合がほとんどであり、上記のような多変量回帰分析手法を用いる必要がある。本講義では、実データの解析を通して、多変量回帰手法を解説する。
重回帰分析（MLR）
本講義では、Gilbertら（2015, PLoS ONE 10: e0138494）のブルーベリーデータを使用する。この研究は、複数のブルーベリー品種を用いた多環境（複数の試験場、複数の栽培年）試験（Multi-environmental trials: MET)において、口当たりなどの官能評価結果と酸度、糖度、多数の揮発成分などの生化学的特性の関係を分析したものである。本講義では、揮発性成分のデータをもとに、ブルーベリーの香りに対する人間の嗜好性を予測するモデルを構築してみる。
まず、ブルーベリーのデータをロードする。
bb <- read.csv("berry.csv", row.names = 1) # load data
head(bb) # check the first six lines## 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次に、ロードされたデータから揮発性成分のデータ（17列目から最後の列まで）を説明変数\(\mathbf{x}\)（ファクター）として抽出する。なお、揮発性成分のデータは、成分によって値の幅が大きく異なるため、分析前にすべての成分について平均0、分散1に標準化する。本講義では、香りの嗜好性を目的変数\(y\)（応答）として用いる。
x <- as.matrix(bb[, 17:ncol(bb)]) # volatile components
x <- scale(x) # scale the variables
y <- bb$Flavor # preference for flavormodel <- lm(y ~ x) # perform multiple regression
summary(model) # see the results of multiple regression##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.7220 -1.1667 0.1598 1.1640 3.9593
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.942384 0.180931 143.383 <2e-16 ***
## xX590.86.3 -0.133870 0.279036 -0.480 0.6325
## xX616.25.1 0.357279 0.811736 0.440 0.6608
## xX107.87.9 -0.021579 0.238015 -0.091 0.9279
## xX110.62.3 0.408335 0.925839 0.441 0.6601
## xX105.37.3 -0.229327 0.277529 -0.826 0.4106
## xX623.42.7 -0.621681 0.386712 -1.608 0.1111
## xX123.51.3 0.662267 0.318438 2.080 0.0401 *
## xX1576.87.0 0.034999 0.790557 0.044 0.9648
## xX71.41.0 -0.277929 0.445528 -0.624 0.5342
## xX1576.95.0 0.400438 0.773677 0.518 0.6059
## xX108.88.3 0.456688 0.409286 1.116 0.2672
## xX556.24.1 -0.094972 0.331164 -0.287 0.7749
## xX66.25.1 -0.008882 0.578791 -0.015 0.9878
## xX6728.26.3 -0.570372 0.761802 -0.749 0.4558
## xX928.95.0 -0.003721 0.540290 -0.007 0.9945
## xX111.27.3 -0.011258 0.663026 -0.017 0.9865
## xX123.92.2 -0.296504 0.383502 -0.773 0.4413
## xX110.43.0 0.503588 0.380999 1.322 0.1893
## xX111.71.7 0.576120 0.494707 1.165 0.2470
## xX1191.16.8 0.363231 0.360385 1.008 0.3160
## xX106.70.7 -0.708879 0.370033 -1.916 0.0583 .
## xX7785.70.8 0.171545 0.311386 0.551 0.5829
## xX928.68.7 -0.560343 0.300606 -1.864 0.0653 .
## xX142.62.1 -0.370256 0.495498 -0.747 0.4567
## xX110.93.0 -0.317338 0.780165 -0.407 0.6851
## xX111.13.7 -0.104884 0.314194 -0.334 0.7392
## xX123.66.0 0.552973 0.307068 1.801 0.0748 .
## xX3681.71.8 -0.101549 0.470853 -0.216 0.8297
## xX142.92.7 -0.214427 0.396072 -0.541 0.5895
## xX2497.18.9 0.725049 0.610392 1.188 0.2377
## xX104.76.7 0.239253 0.238582 1.003 0.3184
## xX5989.27.5 0.840441 0.427130 1.968 0.0519 .
## xX470.82.6 -0.908442 0.475073 -1.912 0.0587 .
## xX122.78.1 -0.049844 0.259562 -0.192 0.8481
## xX821.55.6 0.490430 0.685154 0.716 0.4758
## xX78.70.6 -0.341317 0.423661 -0.806 0.4224
## xX124.19.6 -0.206936 0.278124 -0.744 0.4586
## xX2639.63.6 0.118704 0.315942 0.376 0.7079
## xX53398.83.7 0.075766 0.449456 0.169 0.8665
## xX112.40.3 -0.119190 0.357398 -0.333 0.7395
## xX119.36.8 -0.146940 0.297663 -0.494 0.6227
## xX106.26.3 -0.383523 0.383905 -0.999 0.3202
## xX141.27.5 0.538366 0.449463 1.198 0.2339
## xX112.12.9 -0.236566 0.504087 -0.469 0.6399
## xX629.50.5 -0.042776 0.330439 -0.129 0.8973
## xX3879.26.3 -0.313615 0.402751 -0.779 0.4380
## xX689.67.8 -0.181052 0.939449 -0.193 0.8476
## xX1139.30.6 0.770425 0.536167 1.437 0.1539
## xX5989.33.3 0.246958 0.416238 0.593 0.5543
## xX43219.68.7 -0.154830 0.441315 -0.351 0.7265
## xX582.16.1 0.068757 0.314524 0.219 0.8274
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.223 on 99 degrees of freedom
## Multiple R-squared: 0.5753, Adjusted R-squared: 0.3565
## F-statistic: 2.629 on 51 and 99 DF, p-value: 1.954e-05ここでは、まずMLRで回帰モデルを作成する。
MLRの結果、いくつかの揮発性成分がブルーベリーの香りに対する人の嗜好と有意な関係を持っていることがわかる。また、Multiple R-squared（R2）の値から、人の嗜好のばらつきの57％が同モデルで説明されることがわかる。
次に、推定された回帰係数を見てみよう。
beta.mlr <- coef(model) # extract regression coefficients
beta.mlr # display regression coefficients## (Intercept) xX590.86.3 xX616.25.1 xX107.87.9 xX110.62.3 xX105.37.3
## 25.942384106 -0.133869862 0.357279136 -0.021578980 0.408334958 -0.229327160
## xX623.42.7 xX123.51.3 xX1576.87.0 xX71.41.0 xX1576.95.0 xX108.88.3
## -0.621680679 0.662267282 0.034998576 -0.277928559 0.400438316 0.456688047
## xX556.24.1 xX66.25.1 xX6728.26.3 xX928.95.0 xX111.27.3 xX123.92.2
## -0.094972146 -0.008881508 -0.570371785 -0.003720558 -0.011257983 -0.296503750
## xX110.43.0 xX111.71.7 xX1191.16.8 xX106.70.7 xX7785.70.8 xX928.68.7
## 0.503588127 0.576120002 0.363231061 -0.708879266 0.171544715 -0.560342838
## xX142.62.1 xX110.93.0 xX111.13.7 xX123.66.0 xX3681.71.8 xX142.92.7
## -0.370256118 -0.317337813 -0.104884417 0.552973103 -0.101548866 -0.214426858
## xX2497.18.9 xX104.76.7 xX5989.27.5 xX470.82.6 xX122.78.1 xX821.55.6
## 0.725048765 0.239253316 0.840441396 -0.908441926 -0.049844291 0.490429592
## xX78.70.6 xX124.19.6 xX2639.63.6 xX53398.83.7 xX112.40.3 xX119.36.8
## -0.341316778 -0.206936118 0.118703555 0.075765616 -0.119189709 -0.146939873
## xX106.26.3 xX141.27.5 xX112.12.9 xX629.50.5 xX3879.26.3 xX689.67.8
## -0.383523171 0.538366018 -0.236566205 -0.042776438 -0.313615042 -0.181051979
## xX1139.30.6 xX5989.33.3 xX43219.68.7 xX582.16.1
## 0.770425432 0.246957877 -0.154830018 0.068757195barplot(beta.mlr[-1]) # plot regression coefficients without intercept
いくつかの揮発性成分が、香りの嗜好性にプラスまたはマイナスの影響を及ぼしているようだ。しかし、ここでは多数の要素を含む重回帰を行っているため、得られた結果の解釈には注意が必要である。
ここで、モデルに当てはめて得られた嗜好性の推定値と、実際の観測値との関係を見るために、グラフを描いてみる。
y.fit <- fitted(model) # "estimate" of y (fitted model)
range <- range(c(y, y.fit)) # find the range of observed and estimated values
plot(y, y.fit, xlim = range, ylim = range) # scatter plot of observed y and estimated y.fit
abline(0, 1) # Draw a line y = x on the scatter plot (0 is the intercept, 1 is the slope)
得られた重回帰モデルは、データによく適合しているように見える。観測値と推定値の相関と二乗平均平方根誤差(RMSE)を計算してみよう。
r <- cor(y, y.fit) # Correlation between observed y and estimated y.fit
r # display correlation## [1] 0.7584757rmse <- sqrt(sum((y - y.fit)^2) / length(y)) # root mean square error square root
rmse # display of RMSE## [1] 1.800236上記の計算を一度に行うための関数を自作してみよう。
my.plot <- function(y1, y2, ...) {
range <- range(c(y1, y2))
plot(y1, y2, xlim = range, ylim = range, ...)
abline(0, 1)
r <- cor(y1, y2)
rmse <- sqrt(sum((y1 - y2)^2) / length(y))
print(c(r, rmse))
}
# use the function
my.plot(y, y.fit, xlab = "Oberved", ylab = "Expected", main = "MLR") 
## [1] 0.7584757 1.8002360では、得られた重回帰モデルの予測精度を評価しよう。ここでいう予測とは、モデル構築に用いたデータ（「訓練データ」あるいは「学習データ」とよばれる）とは異なる未知のデータにモデルを当てはめ、要因から応答を求めることである。予測の精度を評価するためには、未知データに対して予測した応答と実際に測定した応答との「答え合わせ」をする必要がある。そのためによく使われるのが交差検証（cross-validation）という手法である。交差検証とは、すべてのデータを訓練データと検証データに分け、前者をモデル構築に、後者を予測精度の評価に使用する方法である。この例では、データを10分割し、そのうち9分割を訓練データ、残り1分割を検証データとし、10分割交差検証を用いている。10分割交差検証では、計算された予測値の精度は分割の無作為性に依存する。実際には、計算された予測精度の平均と分散は、この処理を複数回繰り返すことで算出される（講義では1回のみ行う）。分割数とサンプル数が等しい場合、1つのサンプルを取り除いたデータが訓練データ、除いたサンプルが検証データとなり、これを1個抜き交差検証と呼ぶ。この場合、予測精度は常に同じになるので、交差検証を繰り返す必要はない。
では、交差検証の準備をしよう。全データを無作為（random）に10群に分割する。
n.sample <- length(y) # find length of y (number of samples)
n.sample # display the number of samples## [1] 151n.fold <- 10 # perform 10-part cross-validation
id <- seq(n.sample) %% n.fold + 1# divide by 10 and compute remainder + 1
id # classified into classes 1-10## [1] 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6
## [26] 7 8 9 10 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 1
## [51] 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6
## [76] 7 8 9 10 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 1
## [101] 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6
## [126] 7 8 9 10 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 1
## [151] 2id.rand <- sample(id) # Randomize the classes (non-reciprocal extraction of the same number)
id.rand # display randomized classes## [1] 9 9 9 1 9 10 4 9 8 6 1 5 5 8 3 3 7 1 6 5 9 9 3 3 8
## [26] 6 5 8 9 2 1 5 5 3 7 1 5 8 10 8 6 4 4 4 3 3 2 7 6 8
## [51] 8 8 1 10 1 8 5 9 7 7 6 4 2 6 3 9 10 7 7 10 8 1 1 7 7
## [76] 6 2 5 6 5 2 5 1 2 4 4 4 10 3 2 3 7 9 1 2 10 3 7 7 8
## [101] 4 2 2 1 6 6 2 9 2 9 1 3 9 7 9 10 2 6 4 1 8 6 1 2 5
## [126] 10 8 5 5 10 10 10 10 4 6 2 7 10 3 4 4 8 3 3 6 4 5 10 4 2
## [151] 7乱数は擬似的に生成されるため、乱数の種（seed）を指定することで、一定の規則に従って同じ値の乱数を連続的に生成することができる。ここでは、解析結果を再現できるように、乱数種を設定しています。
set.seed(1234)では、無作為に分割された10群のうち、9群を訓練データ、1群を検証データとして交差検証を実行してみよう。
data <- data.frame(y, x) # prepare data
y.pred.mlr <- rep(NA, n.sample) # container for predictions
for(i in 1:n.fold) { # Repeat for i from 1 to 10
test <- id.rand == i # i-th class becomes validation data
train <- !test # all other classes become teacher data
data.train <- data[train, ] # data for teaching
data.test <- data[test, ] # data for validation
model <- lm(y ~ . , data = data.train) # multiple regression model
y.pred.mlr[test] <- predict(model, data.test) # prediction on data for validation
}
my.plot(y, y.pred.mlr, main = "MLR",
xlab = "Observed", ylab = "Predicted") # Display results using your own functions
## [1] 0.425476 2.936393予測値と観測値の相関は0.42に減少し、RMSEは2.99に増加した。予測値と観測値の散布図でも、推定値と観測値の関係と比較して、関係の強さが低下していることがわかる。特に、予測値では一部のサンプルが過大評価されている。
主成分回帰とPLS（Partial Least Square）回帰
次に、\(\mathbf{X}\)の主成分分析を行い、その主成分スコアと\(\mathbf{y}\)の回帰分析を行う主成分回帰分析(PCR)を行ってみよう。PCRでは、\(\mathbf{X}\)の主成分分析を行い、その主成分得点と\(\mathbf{y}\)の回帰分析を行う。主成分分析により、回帰分析で用いる説明変数の次元が下がり、多重共線性が取り除かれる。
data <- data.frame(y, x) # prepare data
n.comp <- 30 # number of principal components
model <- pls::pcr(y ~ . , n.comp, data = data) # Principal component regression model
summary(model)## Data: X dimension: 151 51
## Y dimension: 151 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 16.174 30.10 38.96 47.16 52.55 57.58 61.96 65.73
## y 9.761 10.72 11.75 13.44 19.18 19.93 20.26 22.22
## 9 comps 10 comps 11 comps 12 comps 13 comps 14 comps 15 comps
## X 68.55 71.15 73.48 75.65 77.65 79.49 81.13
## y 26.55 27.05 29.25 32.98 33.45 33.53 34.25
## 16 comps 17 comps 18 comps 19 comps 20 comps 21 comps 22 comps
## X 82.69 84.04 85.34 86.46 87.53 88.56 89.48
## y 35.95 36.13 41.31 41.86 41.87 43.31 46.03
## 23 comps 24 comps 25 comps 26 comps 27 comps 28 comps 29 comps
## X 90.36 91.22 92.02 92.79 93.51 94.20 94.78
## y 46.43 46.71 46.86 47.72 50.62 51.64 51.95
## 30 comps
## X 95.32
## y 53.04plot(model)
主成分の数が増えれば増えるほど、xとyの変動のより大きな部分が説明されることがわかる。では、主成分の数はどのように決めればよいだろうか。
ここでは、交差検証を用いて主成分数を決定する。あらかじめ少し大きめの主成分数を設定し、交差検証のオプション（validation = “CV”）を選択して再度pcr関数を実行してみよう。
model <- pcr(y ~ ., n.comp + 20, data = data, validation = "CV") # cross-validation
plot(model, "validation", legendpos = "topright") # visualize the result of cross-validation
描かれているグラフは、予測値の二乗平均平方根誤差（RMSEP）である。このグラフから、主成分数は30程度にすることが望ましいことがわかる。赤い破線は、バイアス補正を行った場合のRMSEPである。30個程度の主成分で予測精度が最大になることがわかります。
この結果から、RMSEPが最も低くなる主成分数を以下のように抽出することができる。
model$validation$PRESS # Display the sum of squared predicted residuals## 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps 7 comps 8 comps
## y 1084.686 1066.604 1063.959 1082.195 995.7825 1002.537 1009.145 1021.932
## 9 comps 10 comps 11 comps 12 comps 13 comps 14 comps 15 comps 16 comps
## y 973.2604 965.2623 965.9067 939.7228 940.5229 943.0533 940.0754 941.2773
## 17 comps 18 comps 19 comps 20 comps 21 comps 22 comps 23 comps 24 comps
## y 930.1743 895.8641 883.1367 877.5594 870.978 844.0637 837.7323 850.4012
## 25 comps 26 comps 27 comps 28 comps 29 comps 30 comps 31 comps 32 comps
## y 859.8153 840.3016 840.4764 804.8797 816.8977 833.9352 841.6213 823.4685
## 33 comps 34 comps 35 comps 36 comps 37 comps 38 comps 39 comps 40 comps
## y 837.3731 841.4446 883.2579 889.9753 903.3275 913.1755 909.2198 922.1753
## 41 comps 42 comps 43 comps 44 comps 45 comps 46 comps 47 comps 48 comps
## y 907.4984 909.0544 967.2052 994.8937 1025.705 1044.274 1167.551 1182.03
## 49 comps 50 comps
## y 1211.992 1248.006n.comp <- which.min(model$validation$PRESS) # Number of principal components with minimum predicted residual sum of squares
n.comp## [1] 28では、予測に最適と判断された主成分数でPCRを行い、yの推定値と観測値の関係を見てみよう。
model <- pcr(y ~ . , n.comp, data = data) # Principal Component Regression Model
y.fit <- predict(model)[,,n.comp] # compute the estimate
my.plot(y, y.fit, main = "PCR", xlab = "Observed", ylab = "Fitted") # plot the results
## [1] 0.7186424 1.9208920MLRの結果と比較すると、適合度が若干悪くなり、推定値と予測値の相関係数が減少し、RMSEが増加していることがわかる。
次にPCRで得られた回帰係数をMLRで得られた回帰係数と比較してみよう。
beta.pcr <- as.vector(coef(model, n.comp)) # Extract regression coefficients
plot(beta.mlr[-1], beta.pcr) # compare regression coefficients
abline(0, 1) # y = x-linear
abline(v = 0, lty = 2) # x = 0, dashed line
abline(h = 0, lty = 2) # y = 0, dashed line
係数の推定値には差があるものの、互いに関係があることがわかる。
次に、交差検証による予測精度の評価を行ってみよう。なお、主成分の数は学習データごとに決定する必要がある。全データを用いて主成分数を決定すると、モデル構築の重要な要素である主成分数の決定に未知データの情報が利用されることになる。
y.pred.pcr <- rep(NA, n.sample) # Prepare predictions
for(i in 1:n.fold) { # Repeat for i from 1 to 10
test <- id.rand == i # i-th class becomes validation data
train <- !test # all other classes become teacher data
data.train <- data[train, ] # data for teaching
data.test <- data[test, ] # data for validation
n.comp <- 50 # first set to 50
model <- pls::pcr(y ~ . , n.comp, data = data.train, validation = "CV") # PC regression model
n.comp <- which.min(model$validation$PRESS) # Number of principal components with highest prediction accuracy
y.pred.pcr[test] <- predict(model, data.test)[,,n.comp] # prediction to validation data
}
my.plot(y, y.pred.pcr, main = "PCR",
xlab = "Observed", ylab = "Predicted") # Display results using your own functions
## [1] 0.5435401 2.3853600PCR を使用することで、MLR を使用する場合と比較して予測精度が向上していることがわかる。予測値と観測値の相関は 0.59 であり、RMSEP は 2.29 である。
次に、PLS回帰分析を用いて回帰モデルを構築してみる。
n.comp <- 10 # number of latent variables
model <- pls::plsr(y ~ . , n.comp, data = data) # PLS regression model
summary(model) # display results## Data: X dimension: 151 51
## Y dimension: 151 1
## 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 8 comps
## X 14.44 20.93 27.81 36.91 44.13 47.78 53.25 57.96
## y 25.38 42.21 49.18 52.18 54.04 55.52 56.02 56.44
## 9 comps 10 comps
## X 60.64 63.24
## y 56.83 57.07PCRと同様に成分（潜在変数）の数が増えると\(x, y\)の説明量が増えるが、PLS回帰ではPCRと比較して特に\(y\)の説明量の増加が顕著である。
ここで、PCRの時と同様に交差検証で回帰に使用する成分数を決定してみよう。
model <- pls::plsr(y ~ . , n.comp + 10, data = data, validation = "CV") # Cross-validation
plot(model, "validation", legendpos = "topright") # Graphical representation of cross-validation results
n.comp <- which.min(model$validation$PRESS) # Number of principal components with minimum predicted residual sum of squares
n.comp## [1] 5PCRの場合と比べて必要な成分数が大幅に削減され、4種類程度で十分であることがわかる。
さて、最適化された成分数をもとに、再度PLS回帰を行い、推定値と実測値の関係を調べてみよう。
model <- pls::plsr(y ~ . , n.comp, data = data) # PLS regression model
y.fit <- predict(model)[,,n.comp] # estimation
my.plot(y, y.fit, main = "PLS", xlab = "Observed", ylab = "Fitted") # plot the results
## [1] 0.7351414 1.8726446PCR と同様に、推定値と観測値の関係は MLR よりも弱いことがわかる。また、推定値と観測値の相関、RMSEはPCRとPLSで同程度であることがわかる。
次に、PLSで得られた回帰係数を他の手法で得られた回帰係数と比較してみよう。
beta.pls <- coef(model, n.comp) # extract regression coefficients
par(mfrow = c(1, 2))
plot(beta.mlr[-1], beta.pls) # compare regression coefficients
abline(0, 1); abline(v = 0, lty = 2); abline(h = 0, lty = 2)
plot(beta.pcr, beta.pls) # compare regression coefficients
abline(0, 1); abline(v = 0, lty = 2); abline(h = 0, lty = 2)
par(mfrow = c(1, 1))PCRと同様、回帰係数はMLRと似ているが、回帰係数の関係はPLSとPCRの方が高いことがわかる。
次に交差検証による予測精度の評価をしてみよう。
y.pred.pls <- rep(NA, n.sample)
for(i in 1:n.fold) { # Repeat for i from 1 to 10
test <- id.rand == i # i-th class becomes validation data
train <- !test # all other classes become teacher data
data.train <- data[train, ] # data for teaching
data.test <- data[test, ] # data for validation
n.comp <- 10 # first set to 10
model <- pls::plsr(y ~ . , n.comp, data = data.train, validation = "CV") # PLS regression model
n.comp <- which.min(model$validation$PRESS) # Number of latent variables with highest prediction accuracy
y.pred.pls[test] <- predict(model, data.test)[,,n.comp] # Prediction on data for validation
}
my.plot(y, y.pred.pls, main = "PLS",
xlab = "Observed", ylab = "Predicted") # Display results using your own functions
## [1] 0.5430825 2.3724459PCRと同様に、PLSとPCRの予測精度に大きな差はない。
リッジ回帰とLASSO
次に、リッジ回帰を用いて予測モデルを構築してみよう。リッジ回帰はglmnetという関数で実行できる。alphaというオプションを0にするとリッジ回帰が、1にするとLASSOが実行され、0と1の間の値がelastic netという手法に対応する。
model <- glmnet::glmnet(x, y, alpha = 0) # Ridge regression model
plot(model, "lambda") # Graphical representation of regression coefficients
リッジ回帰は、回帰係数の大きさにペナルティをかけることで、訓練データへの過適合（overfit）を防ぎ、未知データに対する予測精度を向上させるものである。回帰係数の大きさに対するペナルティの大きさは、\(\lambda\)というパラメータで制御されます。上の図を見ると、\(\lambda\)が大きくなるにつれて、回帰係数の大きさに対するペナルティが強くなり、最終的には0に縮小することがわかります。では、ペナルティの強さはどのように決めればよいのだろうか。様々な方法があるが、ここではPCRやPLSで成分数を決定したのと同様に交差検証を用いて決定する。
cv.glmnetという関数を使ってクロスバリデーションを行い、\(\lambda\)の値を最適化しましょう。
model <- glmnet::cv.glmnet(x, y, alpha = 0) # Ridge regression with cross-validation
plot(model) # plot the results
上図の縦軸はRMSEPで、\(log(\lambda) = 0\)付近で最小となることがわかる。この値は次のようにして求めることができる。
min(model$cvm) # Minimum value of RMSEP## [1] 5.388062model$lambda.min # the lambda at that time## [1] 1.294419log(model$lambda.min) # its logarithm## [1] 0.2580623ここで、RMSEPを最小化する\(\lambda\)を用いた\(y\)の推定値を求め、観測値と比較してみよう。
y.fit <- predict(model, newx = x,
s = "lambda.min") # estimate in lambda that minimizes RMSEP
my.plot(y, y.fit, main = "Ridge",
xlab = "Observed", ylab = "Fitted") # Figure.
## [1] 0.7216016 1.9743858観測値と推定値の関係は、PCRやPLSで得られたものと同様であることがわかる。
では、回帰係数を取り出して、他の手法で推定した値と比較してみよう。
beta.rr <- coef(model,
s = "lambda.min") # estimate in lambda that minimizes RMSEP
par(mfrow = c(1, 2)) # Display in 1 row and 2 columns
plot(beta.mlr[-1], beta.rr[-1]) # compare regression coefficients
abline(0, 1); abline(v = 0, lty = 2); abline(h = 0, lty = 2)
plot(beta.pls, beta.rr[-1]) # Comparison of regression coefficients
abline(0, 1); abline(v = 0, lty = 2); abline(h = 0, lty = 2)
par(mfrow = c(1, 1)) # return one row and one columnリッジ回帰の回帰係数推定値は、PLSの回帰係数推定値と非常によく似ていることがわかる。しかし、リッジ回帰では、回帰係数の値が一般に小さくなる傾向があることも分かる。これは、係数の大きさにペナルティを課している結果による。
それでは、交差検証による予測精度の評価をしてみよう。
y.pred.rr <- rep(NA, n.sample)
for(i in 1:n.fold) { # Repeat for i from 1 to 10
test <- id.rand == i # i-th class becomes validation data
train <- !test # all other classes become teacher data
x.train <- x[train, ] # x's for teaching
y.train <- y[train] # of y for teacher
x.test <- x[test, ] # of x for validation
model <- glmnet::cv.glmnet(x, y, alpha = 0) # PLS regression model
y.pred.rr[test] <- predict(model,
newx = x.test, s = "lambda.min") # Prediction on data for validation
}
my.plot(y, y.pred.rr, main = "RR",
xlab = "Observed", ylab = "Predicted") # Display results using your own function
## [1] 0.7189795 1.9788908予測値と観測値の関係は、推定値と観測値の関係と同様に強く、リッジ回帰で得られたモデルは高い予測性を持っていることがわかる。
最後に、LASSOを用いたモデル構築を行ってみよう。
model <- glmnet::glmnet(x, y, alpha = 1) # LASSO regression model
plot(model, "lambda") # Graphical representation of regression coefficients
リッジ回帰のように、LASSOは回帰係数の大きさにペナルティを課します。Ridge regressionが回帰係数の二乗和にペナルティをかけるのに対して、LASSOは絶対値の和にペナルティをかける。このことから、\(\lambda\)が大きくなると係数の収縮（shrinkage）程度が変わってきます。
ここで、交差検証により\(\lambda\)を最適化しよう。
model <- glmnet::cv.glmnet(x, y, alpha = 1) # LASSO with cross-validation
plot(model) # Graphical representation of results
\(log(\lambda)\)に対してRMSEPは-2付近で最小になることがわかる。
では、最適化された\(\lambda\)の下での推定値を計算し、観測値と比較してみよう。
y.fit <- predict(model, newx = x,
s = "lambda.min") # estimate for a given number of principal components
my.plot(y, y.fit, main = "LASSO",
xlab = "Observed", ylab = "Fitted") # Figure.
## [1] 0.7145421 1.9778361推定値と観測値の関係は、リッジ回帰よりは弱いが、PLSやPCRよりは強いことがわかる。
それでは、回帰係数を計算し、他の手法と比較してみよう。
beta.lasso <- coef(model, s = "lambda.min")
par(mfrow = c(2, 2))
plot(beta.mlr[-1], beta.lasso[-1]) # compare regression coefficients
abline(0, 1); abline(v = 0, lty = 2); abline(h = 0, lty = 2)
plot(beta.pls, beta.lasso[-1]) # Comparison of regression coefficients
abline(0, 1); abline(v = 0, lty = 2); abline(h = 0, lty = 2)
plot(beta.pcr, beta.lasso[-1]) # Comparison of regression coefficients
abline(0, 1); abline(v = 0, lty = 2); abline(h = 0, lty = 2)
plot(beta.rr[-1], beta.lasso[-1]) # Comparison of regression coefficients
abline(0, 1); abline(v = 0, lty = 2); abline(h = 0, lty = 2)
par(mfrow = c(1, 1))LASSOの特徴は、他の手法では係数の小さな推定値がほぼゼロになることが多いことである。LASSOで用いられる係数の絶対値の総和に対するペナルティは、小さな変動しか説明しない変数については係数をゼロにする効果がある。この特徴により、重要な変数を選択することができる。逆に、推定された係数が大きい変数については、リッジ回帰よりもゼロに縮まる傾向が小さる。LASSOとPLS、リッジ回帰とPLS、LASSOとリッジ回帰の比較グラフで確認してみよう。
LASSOを用いた場合の交差検証の予測精度を確認してみよう。
y.pred.lasso <- rep(NA, n.sample)
for(i in 1:n.fold) { # Repeat for i from 1 to 10
test <- id.rand == i # i-th class becomes validation data
train <- !test # all other classes become teacher data
x.train <- x[train, ] # x's for teaching
y.train <- y[train] # of y for teacher
x.test <- x[test, ] # of x for validation
model <- glmnet::cv.glmnet(x, y, alpha = 1) # LASSO
y.pred.lasso[test] <- predict(model, newx = x.test,
s = "lambda.min") # predict to data for verification
}
my.plot(y, y.pred.lasso, main = "LASSO",
xlab = "Observed", ylab = "Predicted") # Display results using your own functions
## [1] 0.7161339 1.9669668PLSやPCRよりは精度が高いが、Ridge回帰より少し精度が低いことがわかる。LASSOでは、目的変数と強い関連を持つ変数のみでモデルを構築するため（弱い変数は回帰係数が0に縮まる）、非ゼロの係数を持つ変数の数が減り、その後の測定の効率は良くなる場合がある。
Random Forest
最後に、機械学習法のひとつであるrandom forest を用いて予測モデルを作ってみる。ranger パッケージの関数を用いてモデルを構築する。
model <- ranger::ranger(y ~ ., data = data)
model## Ranger result
##
## Call:
## ranger::ranger(y ~ ., data = data)
##
## Type: Regression
## Number of trees: 500
## Sample size: 151
## Number of independent variables: 51
## Mtry: 7
## Target node size: 5
## Variable importance mode: none
## Splitrule: variance
## OOB prediction error (MSE): 5.675528
## R squared (OOB): 0.2611456交差検証による予測精度の評価をしてみよう。
y.pred.rf <- rep(NA, n.sample)
for(i in 1:n.fold) { # Repeat for i from 1 to 10
test <- id.rand == i # i-th class becomes validation data
train <- !test # all other classes become teacher data
data.train <- data[train, ] # data for teaching
data.test <- data[test, ] # data for validation
model <- ranger::ranger(y ~ ., data = data.train) # PLS regression model
y.pred.rf[test] <- predict(model, data = data.test)$predictions # Prediction on data for validation
}
my.plot(y, y.pred.rf, main = "RF",
xlab = "Observed", ylab = "Predicted") # Display results using your own function
## [1] 0.5115306 2.3805157なお、random forest でも変数の重要度を評価できる。
model <- ranger::ranger(y ~ ., data = data, importance = "permutation")
ranger::importance(model)## X590.86.3 X616.25.1 X107.87.9 X110.62.3 X105.37.3 X623.42.7
## 0.027927949 0.211127336 -0.036062272 0.090233757 -0.002197971 0.086482752
## X123.51.3 X1576.87.0 X71.41.0 X1576.95.0 X108.88.3 X556.24.1
## 0.119426448 0.073188284 0.038221241 0.106157688 0.252131202 0.157476366
## X66.25.1 X6728.26.3 X928.95.0 X111.27.3 X123.92.2 X110.43.0
## 0.139705112 0.059610892 0.064656548 0.031847172 0.070032719 0.265524681
## X111.71.7 X1191.16.8 X106.70.7 X7785.70.8 X928.68.7 X142.62.1
## 0.080155792 -0.007479923 0.134606485 -0.009756599 0.088871864 0.059275582
## X110.93.0 X111.13.7 X123.66.0 X3681.71.8 X142.92.7 X2497.18.9
## 0.207380845 -0.034167713 0.084817091 0.074161611 0.001402483 0.523009683
## X104.76.7 X5989.27.5 X470.82.6 X122.78.1 X821.55.6 X78.70.6
## -0.005098659 0.160789481 -0.014367640 0.012080367 0.136680486 0.423448561
## X124.19.6 X2639.63.6 X53398.83.7 X112.40.3 X119.36.8 X106.26.3
## 0.016605949 -0.002300743 0.089225724 0.295479562 0.063920175 -0.017302208
## X141.27.5 X112.12.9 X629.50.5 X3879.26.3 X689.67.8 X1139.30.6
## 0.079203702 0.124939332 0.050752931 0.130772865 0.136850388 0.299807668
## X5989.33.3 X43219.68.7 X582.16.1
## 0.033007239 0.052742512 0.092103391棒グラフを用いて表示してみる。
op <- par(mar = c(8, 4, 2, 2))
barplot(ranger::importance(model), las = 2)
par(op)手法間の類似度
それでは、それぞれの手法で計算された予測値が、どのような関係になっているかを見てみよう。
y.pred <- data.frame(y.pred.mlr, y.pred.pcr,
y.pred.pls, y.pred.rr, y.pred.lasso, y.pred.rf) # bundle by row
plot(y.pred) # pairwise scatter between methods
PCR-PLS、Ridge regression-LASSOの関係は強いようである。それらと比較すると、重回帰の予測値は異なっている。特に、2つのサンプルは大きな予測値を持っていることがわかる。
最後に、手法間の距離を計算し、クラスター分析で視覚化してみよう。
d <- dist(t(y.pred)) # Calculate the distance between methods based on the predictions
tre <- hclust(d) # cluster analysis
plot(tre) # plot the figure