Rで線形回帰において総当たり法を行う関数を作成する

今回は、前回書いた0と1の組み合わせを全て出力する関数を応用して、説明変数選択の自動化(総当たり法)を行ってみたいと思います。
ちなみに前回書いた記事とはこちら
Rで任意の個数の0と1の組み合わせパターンを出力する関数を作る - バナナでもわかる話

何言ってるか分からないと思うので初めに少し説明から入ります。

説明変数の選択
総当たり法のアイデア
関数の記述

説明変数の選択

被説明変数ベクトル $y$ に対して、説明変数行列 $X$ を使って線形回帰をするというようなことを考えます。
この時、実際の分析を行うにあたって真の構造はわからないので何を説明変数とするかという説明変数選択問題が発生します。

当然、予測をするのか説明をするのかでどんな基準で変数選択をするかは変わってきて、例えばAIC基準を利用したり、lassoを使ってみたり、決定係数を見てみたり色々方法はあるわけですが(lassoを除けば)結局次のような問題が生じます。

「基準は良いがそのアルゴリズムはどう組むの？？全パターン試すの意外としんどくない？」

そこで、普通は逐次選択法といって一定のアルゴリズムに従って逐次的に確認する方法を採用します。この逐次選択法はRではパッケージ化されていて簡単に試すことが出来ます。

しかし、この方法だと全てのパターンを見ることは出来ていません。あくまで時間的なコストと精度との兼ね合いをもって、逐次的に処理をしているにすぎません。そのため、可能なら全てのパターンを見た方が良いわけです。全てのパターンを見ることを総当たり法と呼びます。

総当たり法のアイデア

まず、想定されうる全ての説明変数を含んだ行列 $X$ を用意します。
そして前回の0と1の全パターンを出力する関数を使って、個数が説明変数の種類数と等しくなるような0と1の組み合わせを出力します。
1の部分を用いる説明変数だと考えて、対応する説明変数の組み合わせで回帰を繰り返し行います。

そして特定の基準(今回は修正済み決定係数で比較する)に従って順位付けしたものを出力します。

これで総当たり法が柔軟に行えるはずです。

関数の記述

sentaku=function(X,y){
###前回のsubsetfun関数
	subsetfun=function(kosuu){
		XX=matrix(,2^kosuu,kosuu)
		for(i in 1:kosuu){
			CCC=t(rbind(rep(0,2^(kosuu-i)),rep(1,2^(kosuu-i))))
			XX[,i]=rep(CCC,2^(i-1))
		}
		return(XX)
	}
###データ
	X=as.matrix(X)
	y=as.vector(y)
###説明変数の組み合わせパターンの出力
	subset=subsetfun(ncol(X))
###箱の作成
	kaiki_list=list()
	adj_vec=c()
	X_trans=matrix(0,nrow(X),ncol(X))
###総当たり回帰の実行
	for(i in 1:nrow(subset)){
		subset_row=subset[i,]
		for(ii in 1:nrow(X)){
			X_trans[ii,]=X[ii,]*subset_row  ###掛け算をすることで、0の部分は説明変数が全て0になる
		}
		X_trans=as.data.frame(X_trans)
		adj_vec[i]=summary(lm(y~.,data=X_trans))$adj.r.squared  ###ここを変えれば基準を変えられる
		kaiki_list[[i]]=lm(y~.,data=X_trans)
	}
return(list(kaiki_list=kaiki_list,adj.r.squared=adj_vec))
}

出てきたものを決定係数が小さい物から順に並べるには次のようにすれば良い。

s=sentaku(X,y)
#これで[[2]]に修正済み決定係数の列が出来た。
#決定係数の小さい順に並べてみる
for(i in order(s[[2]])){
	print(s[[2]][i])
}

これで、例えばこの説明変数を除くと急激に決定係数が減少するだとか、逆にこの説明変数はどちらでも良いなんてものが比較できたりします。