ステップワイズ変数選択：増減法と減増法を区別せよ
井口豊（生物科学研究所，長野県岡谷市）
最終更新：2024 年 1 月 5 日
1. はじめに
回帰分析の変数選択でしばしば使われるのが，ステップワイズ法である。その中で，増減法（forward-backward stepwise）と減増法（backward-forward stepwise）を取り上げ， R を利用して，結果の違いを解説した。
増減法は，最も単純なモデル（null model，定数項のみ）に対して，説明変数を 1 個加えることから始めて，その後，何らかの基準（例えば， AIC）によって，説明変数の出し入れを行い，最適な説明変数の組み合わせを探す方法である。
一方，減増法は，最も複雑なモデル（full model，全変数投入）に対して，説明変数を 1 個減じることから始めて，その後は，増減法と同様にして，最適な説明変数の組み合わせを探す方法である。
英語名をそのまま訳すと，それぞれ，前進後退法（forward-backward）と後退前進法（backward-forward）と呼ぶべきかもしれない。
解説者によって，これらを混同している，あるいは，一方しか解説せず，単に，ステップワイズ法としか述べていない場合も多々あるので注意が必要である。例えば，吉川隆範（2018）「多変量解析と統計的因果推論」その p. 8 を見ると，変数増減法 に対して，「変数減少法からスタートする」と述べているが，これは，変数減増法である。
2. R による計算例
二種類のステップワイズ法を R で実行する場合， step 関数を使うが，その direction オプションには，どちらの方法でも， "both" （前進と後退の両方向）を指定する。前述のとおり，初期モデルの設定を変えてやれば，自動的に，増減法と減増法が実行される。 R を使ったネット上や参考書の変数選択の解説では，この点にあまり触れられていない。
ここでは，以下のような仮想的な 5 個の説明変数で， 10 組のデータ（n = 10）を考えて，増減法と減増法による結果の違いを見た。合わせて，パッケージ bestglm の bestglm 関数による総当り法の結果も示した。
#### スクリプト開始
set.seed(1)
n<- 10
y<- rnorm(n)
x1<- rnorm(n)
x2<- rnorm(n)
x3<- y + rnorm(n)
x4<- y + rnorm(n)
x5<- x1 + x2 + rnorm(n)
x6<- x3 + x4 + rnorm(n)
# 定数項モデル
mod.1 <- lm(y ~ 1)
# 全変数モデル
mod.2 <- lm(y ~ x1 + x2 + x3 + x4 + x5 + x6)
# 変数増減法
step(mod.1, direction="both",
scope=list(upper=mod.2, lower=mod.1))
# 変数減増法
step(mod.2, direction="both",
scope=list(upper=mod.2, lower=mod.1))
# 総当り法
library(bestglm)
Xy<- data.frame(x1, x2, x3, x4, x5, x6, y)
bestglm(Xy, IC="AIC")
### スクリプト終了
結果として，選択された説明変数は，次の表 1 のようになった。
| 変数選択法 | 選択変数 |
|---|---|
| 増減法 | x4 |
| 減増法 | x3, x4, x6 |
| 総当り法 | x3, x4, x6 |
標本が小さい（n = 10）という理由もあるが，この例では増減法と減増法の結果が異なることがわかる。 論文では，ステップワイズ（stepwise）としか書かないものも多い。しかし，読者の観点からすれば，どちらのステップワイズなのか，それを明記したほうが良い。
特に， R を利用した論文では，全変数からの減増法を単にステップワイズと称しているような場合がある。この場合，今回の結果のように，多めの説明変数が選択されているかもしれないので，読者も注意が必要である。
3. 特定の変数を残したい場合
上記の例では，最も単純なモデルを定数項のみ，としたが，その応用として，特定の変数を残したモデルも可能である。逆に言えば，やみくもに，全変数を対象として変数選択してはダメだ，ということである。
例えば，上記のデータで， x1 は残したい，という場合のスクリプトは以下のようになる。総当たり bestglm 関数を利用する場合は，引数 TopModels で上位モデル数を指摘して， x1 が選択されるモデルを上位から探っていく。ここでは，デフォルトにもなっているトップ 5 モデルを算出した。
#### スクリプト開始
set.seed(1)
n<- 10
y<- rnorm(n)
x1<- rnorm(n)
x2<- rnorm(n)
x3<- y + rnorm(n)
x4<- y + rnorm(n)
x5<- x1 + x2 + rnorm(n)
x6<- x3 + x4 + rnorm(n)
# 1 変数 x1 モデル
mod.1 <- lm(y ~ x1)
# 全変数モデル
mod.2 <- lm(y ~ x1 + x2 + x3 + x4 + x5 + x6)
# 変数増減法
step(mod.1, direction="both",
scope=list(upper=mod.2, lower=mod.1))
# 変数減増法
step(mod.2, direction="both",
scope=list(upper=mod.2, lower=mod.1))
# 総当り法
library(bestglm)
Xy<- data.frame(x1, x2, x3, x4, x5, x6, y)
# トップ 5 モデル
bestglm(Xy, TopModels = 5, IC="AIC")$BestModels
### スクリプト終了
まず，ステップワイズの結果は，次の表 2 のようになった。
| 変数選択法 | 選択変数 |
|---|---|
| 増減法 | x1, x3, x4, x6 |
| 減増法 | x1, x3, x4, x6 |
増減法でも，減増法でも，x1 に加えて， x3, x4, x6 が選択された。
次に，総当たりで，適合が良い順にトップ 5 のモデルを示す。
表 3. 適合が良い順にトップ 5 モデル
x1 x2 x3 x4 x5 x6 Criterion
1 FALSE FALSE TRUE TRUE FALSE TRUE -18.05885
2 TRUE FALSE TRUE TRUE FALSE TRUE -16.56425
3 FALSE TRUE TRUE TRUE FALSE TRUE -16.52335
4 FALSE FALSE FALSE TRUE FALSE FALSE -16.50985
5 FALSE FALSE TRUE TRUE TRUE TRUE -16.28970
x1 が取り込まれるモデルで最上位は， 2 番目に良いモデルであることが分かる。そして，それに加えて選択されたのは， x3, x4, x6 であった。