線形モデルにおけるパラメータ選択

はじめに

線形モデルは
$$
Y = \beta_0 + \beta_1X_1 + \beta_2X_2 + \cdots + \beta_pX_p + \epsilon
$$
のような関係を持つ。$Y$は応答変数、$X_i$は予測変数である。予測変数の数を$p$、観測データ数を$n$とする。簡単のためこの記事では最小二乗法によってこのモデルを当てはめる。
この記事では以下のような問題を解決するアプローチを紹介する。

1 予測変数に対して観測データが小さいとき、適切な予測ができない。

予測変数の数と観測データ数の大小は以下のような振る舞いの違いを見せる。

• $n\gg p$のとき求められるモデルの分散は小さくなる。テストデータ(学習に使わなかった観測データ)での性能も良い。
• $n > p$のとき過学習を起こしやすい。テストデータでの性能は悪い。
• $n < p$のとき係数が一意に定まらない(連立一次方程式を解く時にパラメータの数だけ式がないと解けないのと同様)。

2 予測変数のうち応答変数との相関が小さいものが存在する。
相関の無いものの係数は0にすることで計算量、モデルの複雑さが減少する。

部分集合選択

予測変数のうち応答変数と関係があると考えられる部分集合を特定し、関係している変数のみを使ってモデルを当てはめる。

最良部分集合選択

予測変数の考えられる全ての組み合わせてモデルに当てはめ、最良の組み合わせを特定する。全ての組み合わせは$2^p$個である。

アルゴリズム

• ヌルモデルを用意する
• $i=1\sim p$まで繰り返す
  • 予測変数がi個あるモデルを${}_pC_i$通り全て作成する
  • ${}_pC_i$通りからRSS(残差平方和$\sum(y_i-\hat{y_i})^2$)が最小になるものまたは$R^2$(決定係数)が最大となるモデルを抜き出す。
• $1\sim p$個の選び出したモデルとヌルモデルのうちベストなものを交差検証やAIC、BICを用いて選ぶ。選ばれた組み合わせが最良の組み合わせとされる。

詳細

この方法は単純でわかりやすいアプローチではあるが、$p$が大きければ大きいほど計算量が大きくなり効率の面でとても悪い。

ステップワイズ法

最良部分集合選択では$p$が大きければ計算量が大きくなることがわかった。ステップワイズ法はこのことを解決するものである。この方法では変数がないモデルから始めて逐次最良の変数を追加する手法である。全ての変数から生成されたモデルがあり、そこから減らしていくものや、そのどちらも取り入れたハイブリットなものもある。

アルゴリズム

• ヌルモデルを用意する
• $i=1\sim p$まで繰り返す
  • 使われていない予測変数を一つ加えた$p-i$通りのモデル作成する。
  • $p-1$通りからRSS(残差平方和)が最小になるものまたは$R^2$(決定係数)が最大となるモデルを抜き出す
• $1\sim p$個抜き出した中でベストなモデルとヌルモデルのうちベストなものを交差検証やAIC、BICを用いて選ぶ。選ばれた組み合わせが最良とされる。

詳細

最良部分集合選択では$2^p$通りを考えていたのに対して、このアプローチでは$1 + \frac{p(p+1)}{2}$通りのみで済む。$p=10$程度であったとしても$1024$と$56$の違いなのでこれによって計算量が大幅に減ったことがわかる。選んできたものがその個数での最適とは限らないので、これによって生成されたモデルは最良部分集合選択によって生成されたモデルと異なることに注意が必要である。変数を増やしていく他に減らしていく方法もあるが、減らしていく方法は$p$が大きい時にいくつかだけ減らすなどして活用することで、計算量を減らして利用できる。この方法は$n<p$の条件では使えないことに注意する必要がある。

評価について

同じ濃度の部分集合ではRSSが最小、$R^2$が最大となるものを最適なモデルとしたが、異なる濃度間での最適を求めるときは異なる評価法を用いた。RSSが最小、$R^2$が最大となるものを選んだ場合は全ての予測変数で生成したモデルが選ばれるからである。このモデルは訓練データでは良い値を求めることができるが、テストデータでも良い値を得られる保証はない。モデルを生成するときはテスト誤差が小さくしたいので、過学習によるバイアスを考慮した調整または、直接的にテスト誤差を推定する方法を行う。

過学習によるバイアス考慮

マローズのCp

$C_p$は予測変数の種類を考慮したMSEでモデルの評価をする指標であり、値が小さいほど評価が高いモデルとされる。これは以下のような数式で計算される。
$$C_p = \frac{1}{n}(RSS+2p\hat{\sigma}^2) = MSE + \frac{2p\hat{\sigma}^2}{n}$$
$\hat{\sigma}^2$は応答変数の測定値に対する誤差の分散の推定値である。数式からわかるように予測変数の数$p$が大きければ第二項が大きくなり、罰則項としての役割を果たしている。

AIC

最尤法によるモデルの当てはめに対して定義される。AICは以下のような数式となる。
$$
AIC = \frac{1}{n\hat{\sigma}^2}(RSS+2d\hat{\sigma}^2)
$$
このようにAICは$C_p$に比例したものである。つまりAICも値が小さいほどモデルの評価が高くなる。

BIC

ベイズ統計によるモデルの当てはめに対して定義される。BICは以下のような数式となる。
$$
BIC = \frac{1}{n\hat{\sigma}^2}(RSS+\log(n)d\hat{\sigma}^2)
$$
このようにBICはAICの2を$\log(n)$に書き換えたものである。これは観測データ数が多いほど罰則項の値を大きくする役割を果たしている。$\log(7) > 2$であるので、AICよりも罰則項の値はほとんどの場合大きくなる。BICも同様に値が小さいほどモデルの評価が高くなる。

自由度調整済み決定係数

通常の決定係数$R^2$は
$$
R^2 = 1-\frac{RSS}{TSS}
$$
と定義される。TSSは総変動である($\sum(y_i-\bar y)^2$)。RSSはモデルの予測変数を加えるほど値が小さくなり、$R^2$は大きくなる。この影響を小さくするために変更した自由度調整済み決定係数は以下のようになる。
$$
R^2 = 1-\frac{\frac{RSS}{n-p-1}}{\frac{TSS}{n-1}}
$$
$p$が0の時は通常のものと同じである。通常の$R^2$では$p$が大きくなるほど第二項の分子が小さくなるので、$(n-p-1)$で割ることによって是正している。この項のおかげで、良い変数を追加したときは小さくなり、悪い変数を追加したときは大きくなる。分母は通常の$R^2$と値が大きく異なるのを防ぐために$n-1$で割った。自由度調整済み決定係数は先に紹介した方法と異なり値が大きいほどモデルの評価が高くなる。

直接テスト誤差を推定

ホールドアウト検証

観測データを訓練用データとテストようデータにランダムに二分し、訓練データでモデルの作成を行い、テストデータでモデルの評価を行う手法である。これによって訓練には使っていない真新しいデータでモデルの評価を行うことができる。

交差検証

ホールドアウト検証を改修した手法であり、データの選択による推定値のばらつき、観測データ数が少ないことによりモデルをうまく生成できないという問題を解決したものである。LOOCVなど交差検証にも様々な手法があるが、k分割交差検証について説明する。

k分割交差検証

観測データをk個に分割し、そのうちの一個をテストデータに、$(k-1)$個を訓練データにする。これをk回繰り返し評価する。i回目の評価は$MSE_i$とすると全体の評価は以下のように行う(MSEである必要はない)。
$$
CV_{(k)} = \frac{1}{k}\sum_{i=1}^k MSE_i
$$
$k=n$のときをLOOCVと言う。kが大きいときは計算量が大きくなってしまうので$k=5,10$がよく使われる。

縮小推定

求められる係数が通常よりも縮小される(正規化)推定である。縮小されることにより分散が小さくなる。
最小二乗法では係数を
$$
RSS = \sum_{i=1}^n(y_i - \beta_0 - \sum_{j=1}^p\beta_j x_{ij})^2
$$
が最小となるようなものとして推定する。これから紹介するアプローチは不要とされる予測変数の係数が0に近い形で出るように改修したものである。

リッジ回帰

リッジ回帰では係数を
$$
\sum_{i=1}^n(y_i - \beta_0 - \sum_{j=1}^p\beta_j x_{ij})^2 + \lambda\sum_{j=1}^p\beta_j^2
$$
が最小となるようなものとして推定する。$\lambda$は0より大きいチューニングパラメータである。$\lambda$が０のときは最小二乗法を同じものとなり、$\infty$に近づくにつれときは全ての係数が0に近づく。$\lambda\sum_{j=1}^p\beta_j^2$は罰則項であり、これにより係数が小さいものが採用される。この手法では計算量による悩みが少なく、幾つかの$\lambda$について解いても最小二乗法の計算時間と遜色がない。
$n$と$p$がおおよそ等しいときは一つのデータが変更されたことによる影響が大きくなる(バイアスは小さいが、分散は大きい)。リッジ回帰では$\lambda$が大きいほどバイアスは大きくなり、分散が減少するので、このようなケースでも最小二乗法よりもリッジ回帰が適している。

ラッソ回帰

リッジ回帰では予測変数全てを使ったモデルを生成した。係数は0に近い値を取るが0になるわけではない。よってリッジ回帰では予測変数の種類が多いとき解くのが困難である。ラッソ回帰はこれを克服したものである。数式は以下のようになる。
$$
\sum_{i=1}^n(y_i - \beta_0 - \sum_{j=1}^p\beta_j x_{ij})^2 + \lambda\sum_{j=1}^p|\beta_j|
$$
リッジ回帰と異なる点は第二項のノルムがl2からl1になっているところだ。こうすることによって係数を0に近づけるのではなく正確に0にすることができる($\lambda$が大きければ)。これによって自動的に適切な予測変数を抽出することができる。ラッソ回帰はリッジ回帰の上位互換のような紹介をしたが、ラッソ回帰では係数が0ではないものを0としてしまうことがあるので、このような場合リッジ回帰の方が性能が良い(他の係数について考えない場合)。

チューニングパラメータの選定

チューニングパラメータ$\lambda$は交差検証などを用いて決められる。

次元削減

$m < p$となるようなm次元部分空間に予測変数を射影する。このm個の射影を予測変数としてモデルを当てはめる。つまり予測変数$X_1, ...X_p$から新たな予測変数
$$
Z_k = \sum_{j=1}^p\phi_{jk}X_j
$$
を生成する方法である(紹介したものは線型結合)。これから応答変数$Y$は
$$
Y = \theta_0 + \theta_1Z_1 + \theta_2Z_2 + \cdots + \theta_mZ_m + \epsilon
$$
のように求められる。$\theta$は新しい変数である。このように変換してあげることによってp次元の問題からm次元の問題へと単純化することができる。係数は以下のようにして求められる。

\sum_{k=1}^m\theta_k Z_k = \sum_{k=1}^m\theta_k\sum_{j=1}^p\phi_{jk}X_j\\
=\sum_{j=1}^p\sum_{k=1}^m\theta_k\phi_{jk}X_j=\sum_{j=1}^p\beta_j X_j\\
\beta_j=\sum_{k=1}^m\theta_k\phi_{jk}

主成分回帰

主成分回帰(PCV)は高次元の説明変数から低次元の特徴を抽出するために使われる手法である。この手法によって抽出した説明変数によるモデルは、元の説明変数によるモデルと比較して良いモデルを得られると考えられる。なぜなら元の説明変数の情報をほとんど全て抽出した説明変数には含まれており、かつ次元が低いので過学習を回避できるからである。新たな説明変数の数$m$は交差検証によって求められる。この手法によって求められたモデルはmが小さければラッソ回避やリッジ回帰に比べ性能が良くなるが、大きければ性能は悪くなる傾向がある。

説明変数について

説明変数を用いて特徴量の結合を行う前にそれぞれの特徴量のスケールを同じスケールにすると良い。スケール化はこのようにして行う。
$$
\tilde{x_{ij}} = \frac{x_{ij}}{\sqrt{\frac{1}{n}\sum(x_{ij}-\bar x_j)^2}}
$$
リッジ回帰なのでもこれを行うとより良い結果を得ることができる。これを行わなければ分散の大きな変数が大きな影響を持つようになってしまう。