4. 予測・回帰：データから連続値を予測する

必要なパッケージ

require(here)
require(pls) # required to use pcr, plsr functions
require(glmnet)
require(ranger)

背景

科学的な研究では、測定・制御が容易な変数を用いて、別の変数の挙動を制御・予測する必要がある場合が多々あります。例えば、肥料や灌漑、日射などの生育条件によって植物の成長を制御する場合が該当します。前者の変数（要因）から後者の変数（応答）を予測・制御するためには、一般的に回帰分析という統計手法が用いられます。また、以下の条件を満たす場合は、重回帰分析（MLR）も適切な方法です。

1.要因数が少ない。 2.要因間に有意な冗長性がない。 3.要因と応答の関係がよく分かっている。

しかし、これら3つの条件のうち1つでも満たされていない場合、MLRモデルは不適切で使い物にならないことがあります。例えば、要因の数が多く、要因と応答の関係がよく分かっておらず、単に予測精度の高いモデルを構築することが目的の場合です。例えば、次のような例が考えられます。

1. 化学者がスペクトルデータから化学物質の化学組成を推定しようとする場合。この場合、因子は数百にも及ぶスペクトルの各波長における測定値であり、隣接する波長間には強い共線性が存在する。
2. 育種家が、ゲノム中に分布する多数のDNA多型のデータから、動物や植物の形質の表現型を予測したい場合。これは、ゲノミック予測と呼ばれる方法である（Meuwissen et al. 2001） この方法を用いると、DNA多型から動物や植物の遺伝的能力を予測することが可能となり（このような選抜方法をゲノミック選抜と呼ぶ）、形質評価に要する時間の短縮や育種の迅速化が期待される。
3. 農学者が、気温や日照時間などの気候条件の時間的変化に対する植物の応答をモデル化したい場合。

本講演では、主成分回帰、部分最小二乗（PLS）回帰、リッジ回帰、LASSOなどの多変量回帰分析について紹介します。これらは、要因の数が多く、要因間に強い共線性がある場合に、予測モデルを構築するための手法です。ここでは、モデルを構築することに重点を置いています。モデルを構築する際に重視するのは「応答を予測すること」であり、必ずしも変数間の潜在的な関係を解明することではないことに注意します。たとえば，主成分回帰，PLS，リッジ回帰は，応答に大きな影響を与えない因子をスクリーニングするのには適していません．しかし，予測が最終目標で，測定する因子の数を制限する現実的な必要性がない場合は，これらの方法は非常に有用なツールとなり得ます．ところで、要因の数がサンプル数よりはるかに大きい場合、MLRモデルを使用することが困難になります。これは “large p small n”問題と呼ばれるものです。一方、上述した4つの手法はいずれもこの問題に対処することができ ます。上記のようなゲノム予測では、ゲノム全体のDNA多型（因子）の数がサンプル数を超える場合がほとんどであり、上記のような多変量回帰分析手法を用いる必要があります。本講義では、実データの解析を通して、多変量回帰手法を解説します。

重回帰分析（MLR）

　本講義では、Gilbertら（2015, PLoS ONE 10: e0138494）で用いられた全てのデータを使用する。このデータを用いて、揮発性成分の量に基づきブルーベリーの香りに対する人間の嗜好性を予測するモデルを構築してみます。

まず、ブルーベリーのデータをロードする。

bb <- read.csv(here("data", "blueberry.csv"), stringsAsFactors = T)   # load data
head(bb) # check the first six lines

##                  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

次に、ロードされたデータから揮発性成分のデータ（20列目から最後の列まで）を説明変数\(\mathbf{x}\)（ファクター）として抽出します。なお、揮発性成分のデータは、成分によって値の幅が大きく異なるため、分析前にすべての成分について平均0、分散1に標準化します。本講義では、香りの嗜好性を目的変数\(y\)（応答）として用います。

x <- as.matrix(bb[, 20:ncol(bb)])  # volatile components
x <- scale(x) # scale the variables
y <- bb$Flavor # preference for flavor

model <- 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.2377 -1.1035  0.2182  1.1823  4.2505 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  26.090244   0.172676 151.094   <2e-16 ***
## xX590.86.3   -0.004047   0.267382  -0.015   0.9880    
## xX616.25.1    0.365328   1.011393   0.361   0.7186    
## xX107.87.9    0.031380   0.223962   0.140   0.8888    
## xX110.62.3   -0.293967   0.870077  -0.338   0.7361    
## xX105.37.3   -0.234993   0.271721  -0.865   0.3890    
## xX623.42.7   -0.524796   0.365005  -1.438   0.1533    
## xX123.51.3    0.616594   0.291950   2.112   0.0369 *  
## xX1576.87.0   0.014019   0.649720   0.022   0.9828    
## xX71.41.0    -0.051013   0.379477  -0.134   0.8933    
## xX1576.95.0   0.236332   0.899344   0.263   0.7932    
## xX108.88.3    0.384118   0.387203   0.992   0.3233    
## xX556.24.1   -0.112571   0.283833  -0.397   0.6924    
## xX66.25.1     0.137948   0.519976   0.265   0.7913    
## xX6728.26.3   0.086646   0.572576   0.151   0.8800    
## xX928.95.0    0.002723   0.505718   0.005   0.9957    
## xX111.27.3    0.102854   0.543895   0.189   0.8504    
## xX123.92.2   -0.128738   0.353715  -0.364   0.7166    
## xX110.43.0    0.799937   0.359170   2.227   0.0280 *  
## xX111.71.7    0.402944   0.540582   0.745   0.4576    
## xX1191.16.8   0.259713   0.320507   0.810   0.4195    
## xX106.70.7   -0.669875   0.345479  -1.939   0.0550 .  
## xX7785.70.8   0.201538   0.319971   0.630   0.5301    
## xX928.68.7   -0.567805   0.273288  -2.078   0.0400 *  
## xX142.62.1   -0.128043   0.404280  -0.317   0.7521    
## xX110.93.0   -0.504116   0.691014  -0.730   0.4672    
## xX111.13.7   -0.130650   0.282783  -0.462   0.6450    
## xX123.66.0    0.687723   0.382614   1.797   0.0750 .  
## xX3681.71.8   0.379985   0.474027   0.802   0.4245    
## xX142.92.7   -0.587334   0.405532  -1.448   0.1504    
## xX2497.18.9   0.453516   0.553456   0.819   0.4143    
## xX104.76.7    0.315275   0.216514   1.456   0.1482    
## xX5989.27.5   0.430133   0.386421   1.113   0.2681    
## xX470.82.6   -0.445748   0.431374  -1.033   0.3037    
## xX122.78.1   -0.039858   0.253892  -0.157   0.8755    
## xX821.55.6   -0.149547   0.518663  -0.288   0.7736    
## xX78.70.6    -0.583889   0.480431  -1.215   0.2268    
## xX124.19.6   -0.166183   0.330788  -0.502   0.6164    
## xX2639.63.6   0.097056   0.309106   0.314   0.7541    
## xX53398.83.7 -0.228021   0.398761  -0.572   0.5686    
## xX112.40.3    0.135563   0.510710   0.265   0.7912    
## xX119.36.8    0.025524   0.260316   0.098   0.9221    
## xX106.26.3   -0.310367   0.354309  -0.876   0.3829    
## xX141.27.5    0.419287   0.418848   1.001   0.3190    
## xX112.12.9    0.054094   0.444709   0.122   0.9034    
## xX629.50.5    0.176002   0.284271   0.619   0.5371    
## xX3879.26.3  -0.478527   0.373210  -1.282   0.2024    
## xX689.67.8    0.098166   0.910028   0.108   0.9143    
## xX1139.30.6   0.793001   0.452894   1.751   0.0827 .  
## xX5989.33.3   0.467656   0.332604   1.406   0.1625    
## xX43219.68.7 -0.448695   0.399875  -1.122   0.2642    
## xX564.94.3   -0.137871   0.293203  -0.470   0.6391    
## xX582.16.1    0.140409   0.299138   0.469   0.6397    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.211 on 111 degrees of freedom
## Multiple R-squared:  0.5698, Adjusted R-squared:  0.3683 
## F-statistic: 2.827 on 52 and 111 DF,  p-value: 2.361e-06

ここでは、まず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 
## 26.090243902 -0.004046873  0.365328076  0.031379883 -0.293967318 -0.234992938 
##   xX623.42.7   xX123.51.3  xX1576.87.0    xX71.41.0  xX1576.95.0   xX108.88.3 
## -0.524795870  0.616594378  0.014019290 -0.051013199  0.236331701  0.384117706 
##   xX556.24.1    xX66.25.1  xX6728.26.3   xX928.95.0   xX111.27.3   xX123.92.2 
## -0.112570936  0.137947608  0.086646036  0.002722914  0.102854415 -0.128738406 
##   xX110.43.0   xX111.71.7  xX1191.16.8   xX106.70.7  xX7785.70.8   xX928.68.7 
##  0.799937046  0.402944212  0.259712522 -0.669874930  0.201538472 -0.567804624 
##   xX142.62.1   xX110.93.0   xX111.13.7   xX123.66.0  xX3681.71.8   xX142.92.7 
## -0.128043384 -0.504116174 -0.130649804  0.687723045  0.379984959 -0.587334352 
##  xX2497.18.9   xX104.76.7  xX5989.27.5   xX470.82.6   xX122.78.1   xX821.55.6 
##  0.453515967  0.315275498  0.430133434 -0.445747925 -0.039858452 -0.149547434 
##    xX78.70.6   xX124.19.6  xX2639.63.6 xX53398.83.7   xX112.40.3   xX119.36.8 
## -0.583889353 -0.166182826  0.097056112 -0.228020582  0.135562758  0.025524348 
##   xX106.26.3   xX141.27.5   xX112.12.9   xX629.50.5  xX3879.26.3   xX689.67.8 
## -0.310366646  0.419286891  0.054093673  0.176002433 -0.478526899  0.098165572 
##  xX1139.30.6  xX5989.33.3 xX43219.68.7   xX564.94.3   xX582.16.1 
##  0.793001469  0.467656438 -0.448695120 -0.137870547  0.140408566

barplot(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.7548588

rmse <- sqrt(sum((y - y.fit)^2) / length(y)) # root mean square error square root
rmse # display of RMSE

## [1] 1.819253

上記の計算を一度に行うための関数を自作してみましょう。

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.7548588 1.8192526

では、得られた重回帰モデルの予測精度を評価しましょう。ここでいう予測とは、モデル構築に用いたデータ（「訓練データ」あるいは「学習データ」とよばれる）とは異なる未知のデータにモデルを当てはめ、要因から応答を求めることです。予測の精度を評価するためには、未知データに対して予測した応答と実際に測定した応答との「答え合わせ」をする必要があります。そのためによく使われるのが交差検証（cross-validation）という手法です。交差検証とは、すべてのデータを訓練データと検証データに分け、前者をモデル構築に、後者を予測精度の評価に使用する方法です。この例では、データを5分割し、そのうち4分割を訓練データ、残り1分割を検証データとし、5分割交差検証を用いています。5分割交差検証では、計算された予測値の精度は分割の無作為性に依存します。実際には、計算された予測精度の平均と分散は、この処理を複数回繰り返すことで算出されます（講義では1回のみ行う）。分割数とサンプル数が等しい場合、1つのサンプルを取り除いたデータが訓練データ、除いたサンプルが検証データとなり、これを1個抜き交差検証と呼 びます。この場合、予測精度は常に同じになるので、交差検証を繰り返す必要はありません。

では、交差検証の準備をします。全データを無作為（random）に5群に分割する。まずは、1〜5のIDをほぼ同じ数だけ生成します。

n.sample <- length(y) # find length of y (number of samples)
n.sample # display the number of samples

## [1] 164

n.fold <- 5 # 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 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3
##  [38] 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5
##  [75] 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2
## [112] 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4
## [149] 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5

無作為にサンプルを5つのグループに分類するために、得られたIDを並び替えます。

id.rand <- sample(id) # Randomize the classes (non-reciprocal extraction of the same number)
id.rand # display randomized classes

##   [1] 5 1 4 1 2 3 5 1 2 3 3 4 2 4 2 4 3 5 3 5 3 1 3 4 2 5 1 5 4 3 3 5 4 3 2 2 2
##  [38] 2 1 4 3 5 3 2 3 2 2 1 2 4 5 5 1 2 2 1 4 5 5 2 5 3 2 1 1 2 2 1 5 1 4 3 4 1
##  [75] 1 1 5 3 1 5 1 1 1 2 2 3 2 3 5 1 1 5 5 1 3 4 4 4 5 2 3 5 5 4 5 4 5 1 2 3 5
## [112] 4 4 1 5 5 4 1 3 1 4 3 4 4 2 2 1 4 2 5 2 3 3 3 3 4 1 4 3 1 3 4 4 4 2 2 4 4
## [149] 5 3 2 5 5 3 1 5 3 2 3 5 1 2 4 4

なお、Rでは乱数は擬似的に生成されます。したがって、乱数の種（seed）を指定することで、一定の規則に従って同じ値の乱数を連続的に生成することができます。ここでは、解析結果を再現できるように、乱数種を設定しています。

set.seed(1234)

では、無作為に分割された5群のうち、4群を訓練データ、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 5
    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
}
head(y.pred.mlr)

## [1] 24.22421 26.18458 25.67558 24.66115 23.79049 29.59364

y.pred.mlrというオブジェクトには、重回帰モデルによる予測値が保存されています。

これを、実際の観察値yと比較します。

my.plot(y, y.pred.mlr, main = "MLR", 
    xlab = "Observed", ylab = "Predicted") # Display results using your own functions

## [1] 0.3923739 3.1473499

予測値と観測値の相関は0.40に減少し、RMSEは2.94に増加しました。予測値と観測値の散布図でも、推定値と観測値の関係と比較して、関係の強さが低下していることがわかる。特に、予測値では一部のサンプルが過大評価されています。

主成分回帰と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: 164 52 
##  Y dimension: 164 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.931    28.54    36.31    43.77    48.89    53.48    57.55    61.25
## y    8.671    10.97    14.56    14.69    14.79    17.76    17.90    18.50
##    9 comps  10 comps  11 comps  12 comps  13 comps  14 comps  15 comps
## X    64.62     67.52     70.10     72.52     74.68     76.75     78.78
## y    21.39     21.50     22.53     27.62     27.75     28.21     36.60
##    16 comps  17 comps  18 comps  19 comps  20 comps  21 comps  22 comps
## X     80.53     82.21     83.71     85.13     86.42     87.63     88.72
## y     36.99     37.76     37.82     40.66     41.08     43.76     43.79
##    23 comps  24 comps  25 comps  26 comps  27 comps  28 comps  29 comps
## X     89.67     90.61     91.51     92.32     93.03     93.70     94.29
## y     46.30     46.31     47.39     47.94     48.11     52.03     53.24
##    30 comps
## X     94.86
## y     53.34

plot(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 1205.423 1178.194 1165.724 1178.029 1169.838 1161.752 1166.358 1174.038
##    9 comps 10 comps 11 comps 12 comps 13 comps 14 comps 15 comps 16 comps
## y 1148.275 1195.473 1174.394 1101.371 1099.716 1097.264 1024.154 1036.499
##   17 comps 18 comps 19 comps 20 comps 21 comps 22 comps 23 comps 24 comps
## y 1056.478 1061.962 1031.782 981.5121 980.5477  970.331   972.35 998.3975
##   25 comps 26 comps 27 comps 28 comps 29 comps 30 comps 31 comps 32 comps
## y 986.4258 954.3053 954.1367 922.8505 898.8521 909.2711 939.7675 914.1738
##   33 comps 34 comps 35 comps 36 comps 37 comps 38 comps 39 comps 40 comps
## y 945.5253 960.1939  935.602 991.1556 996.6053 1011.457 1031.418 1055.206
##   41 comps 42 comps 43 comps 44 comps 45 comps 46 comps 47 comps 48 comps
## y 1114.641 1120.502 1144.087 1144.238 1159.989 1249.578 1262.848 1293.064
##   49 comps 50 comps
## y 1347.087 1400.373

n.comp <- which.min(model$validation$PRESS) # Number of principal components with minimum predicted residual sum of squares
n.comp

## [1] 29

では、予測に最適と判断された主成分数で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.7296593 1.8967048

MLRの結果と比較すると、適合度が若干悪くなり、推定値と予測値の相関係数が減少し、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.5386126 2.4131717

PCR を使用することで、MLR を使用する場合と比較して予測精度が向上していることがわかります。予測値と観測値の相関は 0.47 であり、RMSEP は 2.52 です。

次に、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: 164 52 
##  Y dimension: 164 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.91    22.09    29.55    36.72    41.22    44.00    46.95    51.58
## y    24.48    41.98    49.41    52.31    54.39    55.65    56.11    56.34
##    9 comps  10 comps
## X    55.59     58.27
## y    56.54     56.70

PCRと同様に成分（潜在変数）の数が増えると\(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] 4

PCRの場合と比べて必要な成分数が大幅に削減され、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.7232877 1.9153886

PCR と同様に、推定値と観測値の関係は MLR よりも弱いことがわかります。また、推定値と観測値の相関は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.5274512 2.4186659

PCRと同様に、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.350493

model$lambda.min # the lambda at that time

## [1] 1.006219

log(model$lambda.min) # its logarithm

## [1] 0.006199488

ここで、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.7304473 1.9439990

観測値と推定値の関係は、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.train, y.train, 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.5272258 2.3591620

予測値と観測値の関係は、推定値と観測値の関係と同様に強く、リッジ回帰で得られたモデルは高い予測性を持っていることがわかります。

最後に、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.7012468 2.0290938

推定値と観測値の関係は、リッジ回帰よりは弱いが、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.train, y.train, 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.4462608 2.5224116

PLSやPCRよりは精度が高いが、Ridge回帰より少し精度が低いことがわかります。LASSOでは、目的変数と強い関連を持つ変数のみでモデルを構築するため（弱い変数は回帰係数が0に縮まる）、非ゼロの係数を持つ変数の数が減り、その後の測定の効率は良くなる場合があります。

手法間の類似度

それでは、それぞれの手法で計算された予測値が、どのような関係になっているかを見てみましょう。

y.pred <- data.frame(y.pred.mlr, y.pred.pcr, 
    y.pred.pls, y.pred.rr, y.pred.lasso) # 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