目次

・汎化性能が高いモデルが欲しい

・説明変数が多すぎるとモデルの予測性が下がる

・罰則項付き最小二乗法

・ラグランジュの未定乗数法に帰着して考える

・LASSO回帰/Ridge回帰

・制約の強さλをどう決める？ーk-fold交差検証法

・Rで行う罰則項付き最小二乗法（LASSO回帰）

・おわりに

汎化性能が高いモデルが欲しい

　我々は説明変数をモデルに組み込むことで、その説明変数を使って興味ある現象を予測しようとする。予測の言うからには、今手元にあるデータだけでなく、将来取れるデータでも同じような結果になっていないと意味がないだろう。予測精度の良いモデルとは、どのように構築すればよいのだろう。素朴な感覚としては、今手元のデータで構築したモデルの誤差が小さいものが予測精度がよさそうに思える。ところが、この誤差は説明変数を増やせば増やすほど、言い換えればモデルを複雑にするほど、小さくできる。例え、その説明変数が予測にあまり寄与しなくても、小さくできるのである。もっと厄介なことに、モデルが過剰に複雑だと、新たに得られたデータに関しての予測には全く対応できない状態になる。このような状態は、元のデータだけを過剰に学習した状態、過学習overfittingと呼ばれる。また、傾向として過学習状態では、説明変数の絶対値が大きくなる傾向がある。一方で、説明変数が少なすぎてデータの予測ができない状態を未学習underfittingと呼ぶ。過学習も未学習も予測精度を下げるため、統計モデリングを行う上で避けなければならない。手元のデータだけでなく、将来取れるデータなどでも予測性を担保することを汎化generalizationという。つまり、私たちが最も望んでいるのは、汎化性能の高いモデルなのだ。では、どうやって説明変数は選んでいけばよいのだろう。人間が選べば、そこに恣意的な力が働く可能性も捨てきれない。本稿では、このような変数選択の術として、罰則項付き最小二乗法（正則化最小二乗法）regularized least squares methodについて紹介する。

説明変数が多すぎるとモデルの予測性が下がる

　とはいえ、上記の現象は本当なのか、気になるところだろう。そこでまずは上記の現象の確認を行う。以下のようなデータを解析しよう。

------------------------------------------------------

library(plotn)

x1 <- c(22.3, 18.9, 0.4, 5.5, -1.3, 3.8, 22.1, 12.1, 6.2, 13.9, 6.1, 

        24, 13.2, 18.3, 14.2, -5.4, 9.3, -3.6, -7.5, -8, 7.8, 4.5, 9,

        11.1, -0.5, 1.9, 19.8, 17.9, 14.1, -4.1, -4.7, -4.4, 5.5, 26.8,

        15.6, 9.5, 29.2, 13.7, -5.7, -8.9, -9, 9.3, 12.2, 16.3, 21.4, 

        -4, 0.3, 26.6, 14.7, 2, 13.9, -4.7, 27.8, -3.3, -8.5, 21.5, 20.4)

y <- c(-51.3, -39.6, -7.2, -28.9, -53.4, 11, -74.5, -44.1, 1.6, -2.3, 

       -8.4, -58.8, -41.5, -25.8, -31.9, -0.9, 11.4, -2.3, -1.8, 29.7,

       -15.7, -17.4, -21.8, -26, 16.7, 14.7, -33.2, -47.4, 6.8, 47.6,

       39.3, 18.5, -26.1, -27, -28.9, -11.9, -38.1, 3.1, -5.3, 14.3, 

       -9.9, -7.9, -30.9, 10.8, -46.7, 17.6, 45.1, -62.8, -24.1, -5.7,

       -10, 8.1, -53.6, -7.5, 21.6, -32.2, -32.2)

d <- data.frame(x1 = x1, y = y, group = 1)

head(d)

##     x1     y group

## 1 22.3 -51.3     1

## 2 18.9 -39.6     1

## 3  0.4  -7.2     1

## 4  5.5 -28.9     1

## 5 -1.3 -53.4     1

## 6  3.8  11.0     1

plotn(y ~ x1, data = d) #図1の描画

------------------------------------------------------

　さあ、このデータ、皆さんはどんな回帰を考えるだろう。1次回帰？、ちょっと曲がってる気もするし2次回帰？　正直、データ数も少ないし、決めるのは難しい盤面だ。けれども、生態学ではこういうデータを解析する必要が出てくるものだ。そこで、とりあえずx1の多項式で表される信じて、4次式までの多項式回帰を試してみよう。

------------------------------------------------------

fit1 <- lm(y ~ 1, data = d) #0次回帰

fit2 <- summary(fit1)

fit2

## 

## Call:

## lm(formula = y ~ 1, data = d)

## 

## Residuals:

##     Min      1Q  Median      3Q     Max 

## -60.867 -18.567   3.633  20.433  61.233 

## 

## Coefficients:

##             Estimate Std. Error t value Pr(>|t|)    

## (Intercept)  -13.633      3.657  -3.728 0.000452 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 27.61 on 56 degrees of freedom

fit3 <- lm(y ~ x1, data = d) #1次回帰

fit4 <- summary(fit3)

fit4

## 

## Call:

## lm(formula = y ~ x1, data = d)

## 

## Residuals:

##     Min      1Q  Median      3Q     Max 

## -57.418 -12.369  -2.204  13.313  43.990 

## 

## Coefficients:

##             Estimate Std. Error t value Pr(>|t|)    

## (Intercept)    1.655      3.230   0.512     0.61    

## x1            -1.817      0.235  -7.733 2.38e-10 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 19.29 on 55 degrees of freedom

## Multiple R-squared:  0.5209, Adjusted R-squared:  0.5122 

## F-statistic: 59.79 on 1 and 55 DF,  p-value: 2.383e-10

fit5 <- lm(y ~ x1 + I(x1^2), data = d) #2次回帰

fit6 <- summary(fit5)

fit6

## 

## Call:

## lm(formula = y ~ x1 + I(x1^2), data = d)

## 

## Residuals:

##     Min      1Q  Median      3Q     Max 

## -57.728 -12.189  -0.869  12.052  43.178 

## 

## Coefficients:

##             Estimate Std. Error t value Pr(>|t|)   

## (Intercept)  2.38015    3.38912   0.702  0.48551   

## x1          -1.52074    0.46615  -3.262  0.00192 **

## I(x1^2)     -0.01705    0.02310  -0.738  0.46369   

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 19.37 on 54 degrees of freedom

## Multiple R-squared:  0.5257, Adjusted R-squared:  0.5081 

## F-statistic: 29.92 on 2 and 54 DF,  p-value: 1.797e-09

fit7 <- lm(y ~ x1 + I(x1^2) + I(x1^3), data = d) #3次回帰

fit8 <- summary(fit7)

fit8

## 

## Call:

## lm(formula = y ~ x1 + I(x1^2) + I(x1^3), data = d)

## 

## Residuals:

##     Min      1Q  Median      3Q     Max 

## -59.594 -11.504  -1.864  11.566  41.107 

## 

## Coefficients:

##              Estimate Std. Error t value Pr(>|t|)   

## (Intercept)  4.429609   4.679656   0.947   0.3482   

## x1          -1.438476   0.486083  -2.959   0.0046 **

## I(x1^2)     -0.060399   0.071703  -0.842   0.4034   

## I(x1^3)      0.001495   0.002339   0.639   0.5256   

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 19.47 on 53 degrees of freedom

## Multiple R-squared:  0.5293, Adjusted R-squared:  0.5026 

## F-statistic: 19.86 on 3 and 53 DF,  p-value: 9.268e-09

fit9 <- lm(y ~ x1 + I(x1^2) + I(x1^3) + I(x1^4), data = d) #4次回帰

fit10 <- summary(fit9)

fit10

## 

## Call:

## lm(formula = y ~ x1 + I(x1^2) + I(x1^3) + I(x1^4), data = d)

## 

## Residuals:

##     Min      1Q  Median      3Q     Max 

## -58.265 -12.548  -1.948  10.393  42.026 

## 

## Coefficients:

##               Estimate Std. Error t value Pr(>|t|)

## (Intercept)  3.4219934  5.1516507   0.664    0.509

## x1          -1.1492949  0.7717053  -1.489    0.142

## I(x1^2)     -0.0339957  0.0904585  -0.376    0.709

## I(x1^3)     -0.0028603  0.0092863  -0.308    0.759

## I(x1^4)      0.0001129  0.0002329   0.485    0.630

## 

## Residual standard error: 19.62 on 52 degrees of freedom

## Multiple R-squared:  0.5314, Adjusted R-squared:  0.4954 

## F-statistic: 14.74 on 4 and 52 DF,  p-value: 4.089e-08

------------------------------------------------------

さて、負の相関があるのはほぼ確実なので、1次回帰の結果は納得だろう。問題は2次以上の項を考慮するかである。さて、当てはまりの良さの指標であるR二乗値を確認すると、モデルが複雑にするほど大きな値を示していることがわかるだろう。これらの関係を図示してみよう。

------------------------------------------------------

xx <- seq(min(d$x1), max(d$x1), length = 200)

coefs1 <- coef(fit2)[,1]

yy1 <- rep(coefs1[1], length(xx))

coefs2 <- coef(fit4)[,1]

yy2 <- coefs2[1] + coefs2[2]*xx

coefs3 <- coef(fit6)[,1]

yy3 <- coefs3[1] + coefs3[2]*xx + coefs3[3]*(xx)^2

coefs4 <- coef(fit8)[,1]

yy4 <- coefs4[1] + coefs4[2]*xx + coefs4[3]*(xx)^2 + coefs4[4]*(xx)^3

coefs5 <- coef(fit10)[,1]

yy5 <- coefs5[1] + coefs5[2]*xx + coefs5[3]*(xx)^2 + coefs5[4]*(xx)^3 + coefs5[5]*(xx)^4

plotn(y ~ x1, data = d, legend = T, leg.lab = 0:4, pch.leg = NA, lty.leg = 1, 

      pt.col.leg = c("darkgrey","red","blue","darkgreen","purple")) #図2の描画

overdraw(points(xx, yy1, col = "darkgrey", type = "l"),

         points(xx, yy2, col = "red", type = "l"),

         points(xx, yy3, col = "blue", type = "l"),

         points(xx, yy4, col = "darkgreen", type = "l"),

         points(xx, yy5, col = "purple", type = "l"))

------------------------------------------------------

上の図はどう映るだろうか。確かに3次や4次はデータに合わせて、曲線をフィットしているので、データの平均値を通っている気もする。けれども一方で、本当に真の構造はこんなにうねうねしてるのかという気持ちにもなる。次に、考えられるシンプルな変数選択は、有意な効果があればそれを選択する方法である。ところが、1次の項ですら、モデルによって有意になったりならなかったりだ。果たして、この中から予測性のいいモデルを選べるだろうか。また、確認しておくと、係数の絶対値の和はモデルが複雑になるほど、基本的に大きな値を示す。

------------------------------------------------------

#係数の絶対値の和

sum(abs(coefs1))

## [1] 13.63333

sum(abs(coefs2))

## [1] 3.472766

sum(abs(coefs3))

## [1] 3.917942

sum(abs(coefs4))

## [1] 5.92998

sum(abs(coefs5))

## [1] 4.608257

------------------------------------------------------

　そこで、始めに得たデータと同じ条件で、新たに次のようなデータが取れたとする。このデータも図示してみよう。

------------------------------------------------------

x1_2 <- c(-4.7, 16.8, 28.1, 9.7, 25.5, -4.6, 1.3, 13.7, 5.4, 12.7, 

          12.3, 25.6, 17.8, 24.7, -8.3, 22.4, -5.5, 18.9, 3.3, 26.5,

          27.4, 8.1, 4.5, 5, 19.3, 20.5, 14.7, -8.9, 29.7, 27, -3.8,

          11.7, 0.8, 3.2, 26.1, 8.9, 25.3, -7.7, 0.7, -8.4, 28.9, 24.1,

          -6.4, 15.8, 13.1, 16.8, 8.2, 28.3, 22.9, 1.1, 22.3, 23.6, 

          -7.9, 4.7, -6, 3.8, 10.1)

y_2 <- c(0.1, -11.6, -76, -24.6, -29.4, 21.6, 22.2, -78.4, -11.8, -17.5, 

         -21.9, -47.2, -44.3, -26.9, 15.6, -36.1, 12, -25.1, -19.3, -2.2,

         -46.7, -11.7, 3.5, 0.6, -35.1, -24.4, -32.2, 19, -80.3, -24, 9.5,

         -22.6, -20.9, -15.8, -36.3, 1.5, 0.4, -9.2, 17.5, -3.1, -76.9, 

         -52.4, 23.4, -64.2, 1.9, -30.7, -20.6, -46.3, -5.4, 0, -63.6, 

         -30, 3.6, -27.8, 50.6, -20.7, 14.7)

d2 <- data.frame(x1 = x1_2, y = y_2, group = 2)

d <- rbind(d,d2)

plotn(y ~ x1, data = d, mode = "m", group = "group", pch = c(16, 2), 

      col.dot = c("black","red"),

      legend = T, leg.lab = 0:4, pch.leg = NA, lty.leg = 1, 

      pt.col.leg = c("darkgrey","red","blue","darkgreen","purple")) #図3の描画

overdraw(points(xx, yy1, col = "darkgrey", type = "l"),

         points(xx, yy2, col = "red", type = "l"),

         points(xx, yy3, col = "blue", type = "l"),

         points(xx, yy4, col = "darkgreen", type = "l"),

         points(xx, yy5, col = "purple", type = "l"))

------------------------------------------------------

どうだろう。なんだか、両端部分は結構怪しいんじゃないか、そんな気がしてくるだろう。では、予測精度の指標として、モデルから予測される値と実際のデータの偏差の二乗和を確認してみよう。group 1のデータで行った多項式回帰の係数を使って以下のように計算する。

------------------------------------------------------

#データ1のyの推定値

xe <- d$x1[d$group == 1]

ye1_1 <- coefs1[1]

ye1_2 <- coefs2[1] + coefs2[2]*xe

ye1_3 <- coefs3[1] + coefs3[2]*xe + coefs3[3]*(xe)^2

ye1_4 <- coefs4[1] + coefs4[2]*xe + coefs4[3]*(xe)^2 + coefs4[4]*(xe)^3

ye1_5 <- coefs5[1] + coefs5[2]*xe + coefs5[3]*(xe)^2 + coefs5[4]*(xe)^3 + coefs5[5]*(xe)^4 

error1 <- c(sum((d$y[d$group == 1] - ye1_1)^2),

            sum((d$y[d$group == 1] - ye1_2)^2),

            sum((d$y[d$group == 1] - ye1_3)^2),

            sum((d$y[d$group == 1] - ye1_4)^2),

            sum((d$y[d$group == 1] - ye1_5)^2))

error1

## [1] 42695.57 20456.67 20252.38 20097.53 20007.10

#データ2のyの推定値

xe <- d$x1[d$group == 2]

ye2_1 <- coefs1[1]

ye2_2 <- coefs2[1] + coefs2[2]*xe

ye2_3 <- coefs3[1] + coefs3[2]*xe + coefs3[3]*(xe)^2

ye2_4 <- coefs4[1] + coefs4[2]*xe + coefs4[3]*(xe)^2 + coefs4[4]*(xe)^3

ye2_5 <- coefs5[1] + coefs5[2]*xe + coefs5[3]*(xe)^2 + coefs5[4]*(xe)^3 + coefs5[5]*(xe)^4 

error2 <- c(sum((d$y[d$group == 2] - ye2_1)^2),

            sum((d$y[d$group == 2] - ye2_2)^2),

            sum((d$y[d$group == 2] - ye2_3)^2),

            sum((d$y[d$group == 2] - ye2_4)^2),

            sum((d$y[d$group == 2] - ye2_5)^2))

error2

## [1] 46232.83 21896.20 22444.60 22598.95 23676.39

barplotn(rbind(error1, error2), beside = T, ylab = "error", names = 0:4,

         legend = T, leg.lab = c("group1", "group2"), leg.sp = 4) #図4の描画

------------------------------------------------------

group 1のデータの誤差は次数が大きくなるほど小さくなる。これはR二乗値からも確認できている。ところが、group 2のデータに対する誤差は次数が上がるほど大きくなることがわかる。つまり、高次項の多項式モデルは、group 1を過剰に学習しており、group 2のデータの予測はうまく行ってないのである。この状態がまさに過学習である。

　実は、これらのデータは以下のような設定で、1次式を基盤として作成している。新たにデータを生成して、その時の様子を確認しよう。

------------------------------------------------------

b1 <- 5

b2 <- -2

n <- 57

x1_3 <- runif(n, -10, 30)

y_3 <- rnorm(n = n, mean = b1 + b2*x1_3, sd = 20)

d3 <- data.frame(x1 = x1_3, y = y_3, group = 3)

d <- rbind(d,d3)

plotn(y ~ x1, data = d, mode = "m", group = "group", pch = c(16, 2, 4), 

      col.dot = c("black","red", "red"),

      legend = T, leg.lab = 0:4, pch.leg = NA, lty.leg = 1, 

      pt.col.leg = c("darkgrey","red","blue","darkgreen","purple")) #図5の描画

overdraw(points(xx, yy1, col = "darkgrey", type = "l"),

         points(xx, yy2, col = "red", type = "l"),

         points(xx, yy3, col = "blue", type = "l"),

         points(xx, yy4, col = "darkgreen", type = "l"),

         points(xx, yy5, col = "purple", type = "l"))

------------------------------------------------------

いよいよ怪しい気分だろう。誤差を確認しよう。

------------------------------------------------------

#データ3のyの推定値

xe <- d$x1[d$group == 3]

ye3_1 <- coefs1[1]

ye3_2 <- coefs2[1] + coefs2[2]*xe

ye3_3 <- coefs3[1] + coefs3[2]*xe + coefs3[3]*(xe)^2

ye3_4 <- coefs4[1] + coefs4[2]*xe + coefs4[3]*(xe)^2 + coefs4[4]*(xe)^3

ye3_5 <- coefs5[1] + coefs5[2]*xe + coefs5[3]*(xe)^2 + coefs5[4]*(xe)^3 + coefs5[5]*(xe)^4 

error3 <- c(sum((d$y[d$group == 3] - ye3_1)^2),

            sum((d$y[d$group == 3] - ye3_2)^2),

            sum((d$y[d$group == 3] - ye3_3)^2),

            sum((d$y[d$group == 3] - ye3_4)^2),

            sum((d$y[d$group == 3] - ye3_5)^2))

error3

## [1] 48109.64 18615.40 19748.39 20001.61 20587.06

barplotn(rbind(error1, error2, error3), beside = T, ylab = "error", names = 0:4,

         legend = T, leg.lab = c("group1", "group2", "group3"), leg.sp = 4) #図6の描画

------------------------------------------------------

すると、1次回帰が最もマシな成績をおさめていることがわかる。こんな感じで、過剰に複雑なモデルは必ずしもデータ一般に対して予測性は高いわけではないのである。

罰則項付き最小二乗法

　そこで汎化性能を高めるために考えられた手法が、罰則項付き最小二乗法regularized least squares methodである。正則化最小二乗法とも呼ばれ、英語の直訳はこちらである。ただ、罰則項付き最小二乗法のほうが、解析の様子が非常にイメージしやすい。罰則項を課すことを正則化するregularization、と表現することもある。最小二乗法は、モデルとデータの偏差二乗和を最小にするパラメータを選択するのであった。ところがこれだけだと、モデルをどんどん複雑にすればより小さな偏差平方和を得ることができるため、あたかも過剰に複雑なモデルが予測性が良いかのようになってしまう。そこで、以下の式を最小にするようなパラメータを探索する。

式の前半は今までの最小二乗法でも登場した偏差平方和そのものである。f(xi)と一般の形で表示しているが、わかりにくければb1 + b2 x1,iという形を代入して想像してもらってもよい。後半は罰則項（ペナルティー項、正則化項）と呼ばれる部分で、まさに式の形は罰則項付き最小二乗法である。罰則項は係数の和をλだけ補正して足している。つまり、もし偏差平方和部分で絶対値の大きな係数が多く推定されてしまうと、その分だけRは大きな値を示す、ことになる。ということは、式全体として最小化されるということは、極端に大きな係数が出ないように、場合によっては係数の値を0に制御しつつも、誤差を最小にするようなパラメータを推定していることになる。λは罰則項による制御の強さを表しており、大きくなればその分だけ罰則が強くなるので、より強く大きな係数が出ないような制御が働くことになる。逆にλが小さい、もしくは0であれば、それは通常の最小二乗法そのものである。λがこの解析の肝であるが、この値は人間が決めるハイパーパラメータである。どうやって決めるかは、あとで紹介することとして、上記の式を最小にする、とはどういうことかをもう少し掘り下げよう。

ラグランジュの未定乗数法に帰着して考える

　上記の罰則項について、以下の制約を満たすようにパラメータを求めることを考える。

この後に登場するのは上記の式についてp = 1 or 2の時なので、これらの時にg(b1……bm)が示す領域を考える。簡単のため2パラメータで考えよう。以下の図でグレーで示される領域がg(b1, b2) ≦ 0である。

このようなパラメータに制約があるときの最小値を求めるためには、ラグランジュの未定乗数法を活用するとよいのであった。そこでラグランジュ乗数λLを導入して、ラグランジュ関数Lを以下のように考える。

では、ラグランジュの未定乗数法の手続きにのっとると、

ここでベクトル∇h（ナブラh）と∇gは関数hや関数gの勾配ベクトルであり、勾配ベクトルはその関数の値が上昇する方向に向いているのであった。これを図示すると以下のようになる。

制約条件gはその内側はg < 0であり、外側はg > 0であるから、∇gは領域の外を向くようなベクトルである。一方、hはもともと偏差平方和の式だったので、だいたいは下に凸のおわん型の関数である。上図のようにhの最小値がg < 0の領域の外に位置する時、∇h = ∇gを満たす点における関数hの勾配ベクトルは∇gとは逆方向を向いていることになる。つまり、式に立ち戻ればλL < 0であることがわかる。そこで、λ'L = - λLと置きなおせば、ラグランジュ関数は罰則項付き二乗和の式を使って、以下のように書きなおせる。

Rを最小にするパラメータの条件はRをパラメータで偏微分し、偏微分係数が0になる点だから、上記の式で考えると、

以上から、Rの最小値を求める問題は、制約条件（罰則項）g ≦ 0の下での偏差平方和の最小値を求める問題に帰着できるということである。λ = λ'Lμであり、λ'Lはラグランジュの未定乗数法により求まる値だから、λの設定は、μの設定に依存していることになる。μと制約条件gの関係を図で確認すると、以下のようになる。

制約条件gの「円」の半径はμの逆数に依存しているから、μが小さいほど大きな円、μが大きいほど小さな円の制約条件となる。μとλは比例しているから、λが大きいときは円が小さくなる。逆に、λが小さいほど円が大きくなり、λ = 0とは無限大の半径を持つ円、すなわち、制約がない最小二乗法となる。当然、円の半径が大きいほど、偏差平方和の関数の真の最小値が制約条件g < 0に含まれる可能性が高くなるので、制約なしの最小二乗法と一致する。これらの点から、λ（μ）がまさに制約の強さである。

LASSO回帰/Ridge回帰

　では、罰則項の次数によって、パラメータ推定はどのように変わるのか、概要を見てみよう。特にp = 1の時をLASSO（ラッソ、Least Absolute Shrinkage and Selection Operator）回帰、p = 2の時をRidge（リッジ）回帰と呼ぶ。

上図の左がLASSO回帰、右がRidge回帰である。LASSO回帰では制約条件が四角形となっているため、鋭い頂点が存在する。そのため、通常の最小二乗法の最適パラメータと比較して、パラメータ全体が小さな値をとる、だけでなく効果の弱いパラメータが0になりやすくなっている。一方、Ridge回帰は効果の弱いパラメータを0にするだけの力はなく、全体的にパラメータを小さくするだけである。全体的にパラメータを小さくすることで、モデルは極端なブレを生み出すことが少なくなり、汎化性能が高まるのである。また、詳細には立ち入らないが、以下のようなLASSO回帰とRidge回帰をα:(1-α)の割合で混ぜた罰則を科すことで、これら2つの回帰のいいとこどりを目指すのが、エラスティック・ネットElastic netと呼ばれる方法である。

制約の強さλをどう決める？ーk-fold交差検証法

　制約の強さλは、k-fold交差検証法 k-fold cross-validationと呼ばれる方法で決定する。交差検証法をクロス・ヴァリデーションと呼ぶこともある。上記の例では、データの生成過程がわかっていたので、新たにデータを生成してゆき、未知のデータでの当てはまりが悪くなることを示した。しかし、実際は1つのデータから予測性の高いモデルを構築せねばならない。そこで、手元のデータをk分割することを考える。分割したデータに個別にD1、D2……Dkと名前を付けることとして、1回目はD1を除いたデータを学習用として回帰を行い、パラメータを推定する。その後、推定したパラメータを使ったモデルとD1検証用データの偏差平方和（平均二乗誤差）を求める。次に、D2を除いたデータでパラメータを推定し、……とこれをDkまで繰り返す。例えば、k = 5であれば以下のようになる。最終検証用データに関しては後述する。

するとk個分の平均二乗誤差が得られる。平均的に小さな平均二乗誤差を得るには、学習用データに合わせるのではなくパラメータ推定に使わなかった検証用データとの平均二乗誤差が小さくならなければならない。このように、手元のデータを疑似的に既知データ（学習用データ）と未知データ（検証用データ）に分けることで、汎化を達成するのである。この手続きを様々なλで実行し、平均的に小さな平均二乗誤差が得られたλをハイパーパラメータとして用いる。

　なお、今回はそこまでの検証をしないが、実はこれだけでは汎化は不十分である。k個分割したとは言えども、すべてのデータを最終的には学習に使ってしまっているので、このk-foldの中で過学習しているかもしれないからである。したがって、厳密な評価には、学習に一度も使わないデータを最終検証用としてより分けておいて、出来上がったモデルと最終検証用データとの誤差をモデルの性能として評価するべきである。

Rで行う罰則項付き最小二乗法（LASSO回帰）

　では、Rで罰則項付き最小二乗法を実行してみよう。特に、効果の弱い項を削除するモデル選択に興味があるので、LASSO回帰を行う。glmnetパッケージが必要なので、あらかじめインストールしておこう。解析のデータは上記の例で出したものとする。glmnet関数の引数、alphaによって実行する回帰が異なり、alpha = 1ならLASSO回帰、0 < alpha < 1ならエラスティック・ネット、alpha = 0ならRidge回帰である（上記のエラスティック・ネットの式参照）。glmnet関数をそのまま実行すれば、様々なλでの解析結果を返してくれる。

------------------------------------------------------

library(glmnet)

x1 <- c(22.3, 18.9, 0.4, 5.5, -1.3, 3.8, 22.1, 12.1, 6.2, 13.9, 6.1, 

        24, 13.2, 18.3, 14.2, -5.4, 9.3, -3.6, -7.5, -8, 7.8, 4.5, 9,

        11.1, -0.5, 1.9, 19.8, 17.9, 14.1, -4.1, -4.7, -4.4, 5.5, 26.8,

        15.6, 9.5, 29.2, 13.7, -5.7, -8.9, -9, 9.3, 12.2, 16.3, 21.4, 

        -4, 0.3, 26.6, 14.7, 2, 13.9, -4.7, 27.8, -3.3, -8.5, 21.5, 20.4)

y <- c(-51.3, -39.6, -7.2, -28.9, -53.4, 11, -74.5, -44.1, 1.6, -2.3, 

       -8.4, -58.8, -41.5, -25.8, -31.9, -0.9, 11.4, -2.3, -1.8, 29.7,

       -15.7, -17.4, -21.8, -26, 16.7, 14.7, -33.2, -47.4, 6.8, 47.6,

       39.3, 18.5, -26.1, -27, -28.9, -11.9, -38.1, 3.1, -5.3, 14.3, 

       -9.9, -7.9, -30.9, 10.8, -46.7, 17.6, 45.1, -62.8, -24.1, -5.7,

       -10, 8.1, -53.6, -7.5, 21.6, -32.2, -32.2)

d <- data.frame(x1 = x1)

d$x1_2 <- x1^2

d$x1_3 <- x1^3

d$x1_4 <- x1^4

d$y <- y

fit_lasso1 <- glmnet(x = as.matrix(d[,1:4]), y = d$y, family = "gaussian", alpha = 1)

#alpha = 1とすることでLASSO回帰を実行

plot(fit_lasso1, xvar="lambda", label=TRUE)#図7の描画

------------------------------------------------------

では、クロス・バリデーションを行い、適切なλを選ぶことにする。cv.glmnet関数を使う。引数nfoldsを変更することでクロス・バリデーションの回数を変更できる。デフォルトは10なので10-folds クロス・バリデーションである。

------------------------------------------------------

fit_lasso_cv1 <- cv.glmnet(x = as.matrix(d[,1:4]), y = d$y, family = "gaussian", alpha = 1)

plot(fit_lasso_cv1)#図8の描画

------------------------------------------------------

平均的に最小の平均二乗誤差を与えるλは以下のように取り出せる。またその時の各係数は以下のようになる。.で示される項は0であることを示す。この解析によって、3次および4次の項を除去することができた。

------------------------------------------------------

fit_lasso_cv1$lambda.min

## [1] 0.4780327

coef(fit_lasso_cv1, s="lambda.min")

## 5 x 1 sparse Matrix of class "dgCMatrix"

##                      s1

## (Intercept)  1.96095034

## x1          -1.49877346

## x1_2        -0.01580933

## x1_3         .         

## x1_4         .

------------------------------------------------------

これらの係数は制約なしの場合と比べて、その絶対値の和が小さくなっていることがわかる。

------------------------------------------------------

coefs3

## (Intercept)          x1     I(x1^2) 

##  2.38015159 -1.52074032 -0.01705034

coefs5

##   (Intercept)            x1       I(x1^2)       I(x1^3)       I(x1^4) 

##  3.4219933835 -1.1492949112 -0.0339956721 -0.0028602510  0.0001129023

------------------------------------------------------

最後に今回推定されたパラメータを使って回帰曲線を描くと以下のとおりである。

------------------------------------------------------

xx <- seq(min(d$x1), max(d$x1), length = 200)

coefs_lasso1 <- coef(fit_lasso_cv1, s="lambda.min")[,1]

yy <- coefs_lasso1[1] + coefs_lasso1[2]*xx + 

  coefs_lasso1[3]*(xx)^2 + coefs_lasso1[4]*(xx)^3 + coefs_lasso1[5]*(xx)^4

plotn(y ~ x1, data = d, pch.leg = NA, lty.leg = 1)#図9の描画

overdraw(points(xx, yy, type = "l"))

------------------------------------------------------

以上でLASSO回帰を使うことで、汎化性能が高いと思われるモデルを構築できた。さて、重要な注意点を述べる。お気づきかと思うが、LASSO回帰は、あくまで未知のデータに対してもある程度、予測精度を担保したモデルを構築する手段であり、真のモデルを推定するものではない。上記がいい例で、データの生成過程は1次式だが、モデル選択で選ばれたのは2次式である。得られたモデルが真のモデルであると信じることは、とんでもない間違いを産む可能性があるので、十分注意しよう。

おわりに

　今回紹介した解析は私が知る限り、あまり生態学では用いられていないように思う。どちらかといえば、これらの解析法は昨今話題の機械学習の分野でよく用いられている。これらの解析の理解で、上記のことを知っておくことに損はないだろう。また、罰則項という発想は今後、生態学界隈のモデル選択でより頻繁に使われるAICでも登場するので、心の隅に留めておいてほしい。