リッジ回帰(ridge)で解決できる多重共線性問題と注意すべき点

リッジ回帰の記事でも書こうと思って、キーワード検索で「リッジ回帰」と打ち込んでみたら、

「リッジ回帰多重共線性」

という検索候補が出て来たのでちょっとかいてみることにしました。

多重共線性記事について書いていたところだったので、丁度いいですね。
多重共線性とは~回避の方法として相関を見るだけでは..... - バナナでもわかる話
 多重共線性を回避するためのVIF基準の解釈について - バナナでもわかる話

先に結論

先に私の結論ですが、

「このあたりの議論は、やりたいことが機械学習なのか、統計学なのかで対応が異なる」

という結論になりました。

多重共線性問題

ここでまず強い多重共線性と弱い多重共線性について定義しておこうと思います。

多重共線性には、大きく分けて2種類あって、
①説明変数間に完全な線形関係があり、従来の線形回帰分析では解が求まらないケース
②説明変数間に線形関係があるものの、確率的な誤差が混じり、一応解は求まるが、推定量の分散が大きくなってしまうケース、特に不偏推定量の場合、分散は推定量の質に直結するので、これは問題です。

①の、そもそも解が求まらない状況を強い多重共線性

②の、推定結果が不安定になる状況を弱い多重共線性

と呼んでおきます。
詳しくは過去記事にて、もう少し詳しく説明しています。
多重共線性とは~回避の方法として相関を見るだけでは..... - バナナでもわかる話

リッジ回帰と過去記事

検索してて見つけたのは次のような記事です。
強い意味での多重共線性問題とリッジ回帰の関係がわかりやすく説明されています。
リッジ回帰による多重共線性の問題回避について - 統計学と疫学と時々、助教生活

この方が一番わかりやすく短めに記事にされていたので、載せさせていただきました。

強い意味での多重共線性問題が発生している際はそもそも解が不定になるので、何かしらの意味で推定結果を求めることの出来るリッジ回帰は確かに有用そうです。

ただ、この記事の内容自体は正しいとは思うのですが、①の話しかしていません。②はどうでしょう。
そして、仮に②もうまく解決できたとして、このリッジ回帰はどういう点で有用なのでしょう。

私が感じた疑問

確かにridge回帰は解を計算する際の固有値が0より大きくなるため、解は求まるわけですが、

その解がどういう意味で適当な解なのかは理解して使わなければいけないと思います。

これは、恐らく分析の文脈が「機械学習」なのか「統計学」なのかによると思っていて、

機械学習における線形回帰の目的は非説明変数の予測であって、
パラメータの正確な推定値ではなく、MSE(平均二乗誤差)の小ささが重要になります。

リッジ回帰では、コストパラメータ $\lambda$ を使って各変数の影響を均等に潰します。

これにより、例えば外れ値の影響も潰れるため、もし外れ値があった場合に変な予測をする過学習を回避することが可能になります。
※必ずしも有効とは限らない。言うてリッジ回帰では外れ値の影響を潰しきれないことも知られている

また、機械学習の文脈では、②の意味での多重共線性は正直致命的な問題ではありません。要は良い予測が出来ればよいわけなので、 $\lambda$ でうまく潰せて、結果を出力出来て、更に良い予測(例えばMSEや予測誤差が小さい)が出来れば問題ないわけです。

統計学における線形回帰の目的の一つには非説明変数に与える説明変数の影響の推定もあるため、パラメータの推定値も重要になります。
よって、この文脈での多重共線性は①は当然のことながら②も重要になります。

更にリッジ回帰を使って値は出せても推定値にバイアスが出てしまっては、解釈が難しくなりそうです。一致性を備えて、サンプルサイズが大きい場合しかまともに解釈出来なさそうです。

つまり、ネット上の記事をあさっていて思った感想として、

確かに①の問題を解決出来るという点で、機械学習の文脈であればリッジ回帰は有用そうだが、統計学の文脈で語るなら②の意味での多重共線性も解決して、なおかつバイアスが小さい状態が望ましい。

そう思うと、統計学の文脈においては、わざわざリッジ回帰を持ち出さずとも、他の対応法を取れば良いわけで

「リッジ回帰で多重共線性を解消できるよー」

という触れ込みは、統計学と機械学習の区別のない初心者相手には少々危ないというか、誤解を招く恐れがあるのではと思いました。

リッジ回帰による①の意味での多重共線性の解決

①の状態の説明変数を用意して、実験してみましょう。

X1=as.matrix(sample(1:1000,100,replace=TRUE))
X2=as.matrix(sample(1:1000,100,replace=TRUE))
X3=0.8*X1

X=cbind(X1,X2)
X=cbind(X,X3)
Y=1.4+0.9*X1+0.01*X2+0.5*X3+rnorm(100,0,1)

$X_1$ と $X_3$ が従属関係にあるため、従来の回帰分析では当然解は求まりません。試しに、線形回帰分析を行う関数lmを使って、回帰分析してみます。

> lm(Y~X)

Call:
lm(formula = Y ~ X)

Coefficients:
(Intercept)           X1           X2           X3  
    1.49489      1.29965      0.01015           NA

内部でどういう処理をしているのかわかりませんが、まあとにかく $X_3$ の部分はNAになりました。

恐らく不定になり、計算が出来なくなったため、 $X_3$ を除いて回帰してくれたんだと思います。
実際、 $X_1$ の係数は $0.9+0.5*0.8=1.3$ に近い値をとっています。

同じデータにリッジ回帰をやってみましょう。
コストパラメータ $\lambda$ を選択するのが面倒なので、今回は0から1000まで0.1刻みで全部出力してプロットしてみます。

とりあえずリッジ回帰の各推定値を出力する関数を作成します。

ridge=function(X,y,lambda){
	n=nrow(X);
	p=ncol(X);
	X=as.matrix(X);
	x.bar=array(dim=p);
	for(j in 1:p){
		x.bar[j]=mean(X[,j]);
	}
	for(j in 1:p){
		X[,j]=X[,j]-x.bar[j];
	}
	y=as.vector(y);
	y.bar=mean(y);
	y=y-y.bar;
	beta=as.vector(solve(t(X)%*%X+n*lambda*diag(p))%*%t(X)%*%y)
	beta.0=y.bar-x.bar%*%beta;
	return(list(beta=beta,beta.0=beta.0))
}

パラメータを変化させて、値を出してみます。

pl1=vector()
pl2=vector()
pl3=vector()
for(lambda in seq(1,1000,0.1)){
	res=ridge(X,Y,lambda)
	pl1[which(seq(1,1000,0.1)==lambda)]=res$beta[1]
	pl2[which(seq(1,1000,0.1)==lambda)]=res$beta[2]
	pl3[which(seq(1,1000,0.1)==lambda)]=res$beta[3]
}

par(mfcol=c(2,2))
plot(seq(1,1000,0.1),type="l",pl1,xlab="lambda",ylab="推定値",main="beta1")
plot(seq(1,1000,0.1),type="l",pl2,xlab="lambda",ylab="推定値",main="beta2")
plot(seq(1,1000,0.1),type="l",pl3,xlab="lambda",ylab="推定値",main="beta3")

f:id:bananarian:20181027200627p:plain

$\lambda$ が大きくなると値がどんどん小さくなっています。

確かに不定にはなっていませんが....値はどうでしょう。元の値と比べて随分離れた値を取っていませんか。
これは、コストパラメータによって値が潰れている証拠です。

以上のことから、確かに値が不定になる問題は解決できたけども、統計学の意味での推定に、積極的に用いるべきかは疑問です。
確かにチューニング次第でMSEの点で最小二乗法より優れた解を導くことは出来るとは思いますが、説明変数を除いて実行したlmの方がバイアスが無く、実際近い値が出て来ています。

②の意味での多重共線性とリッジ回帰

改めてデータを作ってみます。

X1=as.matrix(sample(1:1000,100,replace=TRUE))
X2=as.matrix(sample(1:1000,100,replace=TRUE))
X3=0.8*X1+rnorm(100,0,1)

X=cbind(X1,X2)
X=cbind(X,X3)

Y=1.4+0.9*X1+0.01*X2+0.5*X3+rnorm(100,0,1)

線形回帰をやってみます。

> lm(Y~X)

Call:
lm(formula = Y ~ X)

Coefficients:
(Intercept)           X1           X2           X3  
   1.827952     0.995931     0.009426     0.379831

不定にはなっていませんね。ただ、推定値の値が大きくはずれていますね(偶然の可能性もありますが)。
弱い多重共線性によって、分散が大きくなった結果と考えるのが自然でしょう。

ではリッジをやってみます。

pl1=vector()
pl2=vector()
pl3=vector()
for(lambda in seq(1,1000,0.1)){
	res=ridge(X,Y,lambda)
	pl1[which(seq(1,1000,0.1)==lambda)]=res$beta[1]
	pl2[which(seq(1,1000,0.1)==lambda)]=res$beta[2]
	pl3[which(seq(1,1000,0.1)==lambda)]=res$beta[3]
}