多重共線性
多重共線性とは、ある回帰モデルにおいて、独立変数間の相関が強い時、その係数の推定値の分散が大きくなってしまい、推定値が安定しない問題である。
最小2乗法による回帰分析で実際に試してみた。
手順は以下の通り。
①従属変数yを正規分布から発生
②yと相関係数0.5の関係にある独立変数x1を作成
＜以下1000回繰り返し＞
③独立変数x1と相関係数rの関係にある独立変数x2を作成
④y~x1+x2の形で回帰分析
⑤x1とx2の偏回帰係数の推定値を取り出す
＜以上1000回繰り返し＞
以上の工程を、rを-1から0.1ずつ+1まで（合計21個）用意して行う。
よって、21×1000回分の偏回帰係数の推定値を作成した。
この時の各rの時の1000個分の偏回帰係数がどれくらい分散しているのかをヴァイオリンプロットで図示した。
なお、n=1000である。
結果は以下の通り。
まず、x1の偏回帰係数から。
これをみると相関係数の絶対値が大きくなるにつれて、その値が変動していることがわかる。
なお、 r=-1、r=1の時はx2の偏回帰係数はNAとなり推定値が出なかったため、x1の偏回帰係数はどのパターンでも同じになる。
ただし、全てプラスの値を取っているので、その強弱の解釈はできないものの、今回の結果ではプラスかマイナスかの判断は過ちは起きないことになる。
x2の偏回帰係数の結果は以下の通り。
こちらは、特に相関係数の絶対値が0.6~0.7になると急に分散が大きくなっているようにみえる。
そもそもx2はyに対してランダムに数値を発生させているので、相関が低くとも初めから分散しており、値も0を中心に正規分布に近いかたちをとっている。
なお、そのp値を取り出してみて同じようにプロットしてみると、以下のようになった。
まず、x1のp値。
なぜか、相関係数がの絶対値が0.9の時だけおかしなことになっているが、縦軸をみると非常に小さい値であり、通常p=0.05で帰無仮説を棄却することを考えれば、p値には影響はないようである。
次に、x2のp値。
これはいずれの相関であっても0から1までほとんど同じように分散している。
以上から、今回の実験では、偏回帰係数にはほとんど相関が強くなればなるほどばらつく現象がみられたものの、p値にはほとんど影響を及ぼさないことがわかった。
なお、以上の実験のコードは以下の通り。
まずは偏回帰係数の方。
n<-1000 #サンプル数
N<-1000 #試行回数
r<-seq(-1,1,0.1) #-1~+1まで0.1ずつ増加させた相関係数
beta1<-matrix(nrow=N,ncol=length(r)) #x1の偏回帰係数の収納
beta2<-matrix(nrow=N,ncol=length(r)) #x2の偏回帰係数の収納
a<-0.5 #従属変数とx1間の相関係数
y<-rnorm(n) #従属変数
u<-rnorm(n)
x1<-a*y+sqrt(1-a^2)*u #独立変数x1
k<-1 #収納する列
for(j in 1:length(r)){
for(i in 1:N){
z<-rnorm(n)
x2<-r[j]*x1+sqrt(1-r[j]^2)*z #x1と相関係数rの関係にある独立変数x2
lm<-lm(y~x1+x2) #最小二乗法による回帰分析
result<-lm$coefficient #偏回帰係数の取り出し
beta1[i,k]<-result[2] #独立変数x1の偏回帰係数
beta2[i,k]<-result[3]#独立変数x2の偏回帰係数
}
k<-k+1
}
#図示
library(ggplot2)
data1<-data.frame(beta1)
data2<-data.frame(beta2)
g<-ggplot()
g<-g+geom_violin(mapping=aes(x="r=-1",y=X1),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.9",y=X2),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.8",y=X3),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.7",y=X4),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.6",y=X5),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.5",y=X6),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.4",y=X7),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.3",y=X8),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.2",y=X9),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.1",y=X10),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0",y=X11),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.1",y=X12),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.2",y=X13),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.3",y=X14),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.4",y=X15),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.5",y=X16),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.6",y=X17),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.7",y=X18),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.8",y=X19),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.9",y=X20),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=1",y=X21),data=data1,position="identity")
print(g)
g<-ggplot()
g<-g+geom_violin(mapping=aes(x="r=-1",y=X1),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.9",y=X2),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.8",y=X3),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.7",y=X4),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.6",y=X5),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.5",y=X6),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.4",y=X7),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.3",y=X8),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.2",y=X9),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.1",y=X10),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0",y=X11),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.1",y=X12),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.2",y=X13),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.3",y=X14),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.4",y=X15),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.5",y=X16),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.6",y=X17),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.7",y=X18),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.8",y=X19),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.9",y=X20),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=1",y=X21),data=data2,position="identity")
print(g)
以下はp値
n<-1000 #サンプル数
N<-1000 #試行回数
r<-seq(-1,1,0.1) #-1~+1まで0.1ずつ増加させた相関係数
p.value1<-matrix(nrow=N,ncol=length(r)) #x1の偏回帰係数の収納
p.value2<-matrix(nrow=N,ncol=length(r)) #x2の偏回帰係数の収納
a<-0.5 #従属変数とx1間の相関係数
y<-rnorm(n) #従属変数
u<-rnorm(n)
x1<-a*y+sqrt(1-a^2)*u #独立変数x1
k<-1 #収納する列
for(j in 1:length(r)){
for(i in 1:N){
z<-rnorm(n)
x2<-r[j]*x1+sqrt(1-r[j]^2)*z #x1と相関係数rの関係にある独立変数x2
lm<-summary(lm(y~x1+x2)) #最小二乗法による回帰分析のsummary
result<-lm$coefficient[,"Pr(>|t|)"] #p値の取り出し
p.value1[i,k]<-result[2] #独立変数x1のp値
p.value2[i,k]<-result[3]#独立変数x2のp値
}
k<-k+1
}
#図示
library(ggplot2)
data1<-data.frame(p.value1)
data2<-data.frame(p.value2)
g<-ggplot()
g<-g+geom_violin(mapping=aes(x="r=-1",y=X1),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.9",y=X2),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.8",y=X3),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.7",y=X4),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.6",y=X5),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.5",y=X6),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.4",y=X7),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.3",y=X8),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.2",y=X9),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.1",y=X10),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0",y=X11),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.1",y=X12),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.2",y=X13),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.3",y=X14),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.4",y=X15),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.5",y=X16),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.6",y=X17),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.7",y=X18),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.8",y=X19),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.9",y=X20),data=data1,position="identity")
g<-g+geom_violin(mapping=aes(x="r=1",y=X21),data=data1,position="identity")
print(g)
g<-ggplot()
g<-g+geom_violin(mapping=aes(x="r=-1",y=X1),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.9",y=X2),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.8",y=X3),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.7",y=X4),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.6",y=X5),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.5",y=X6),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.4",y=X7),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.3",y=X8),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.2",y=X9),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=-0.1",y=X10),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0",y=X11),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.1",y=X12),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.2",y=X13),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.3",y=X14),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.4",y=X15),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.5",y=X16),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.6",y=X17),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.7",y=X18),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.8",y=X19),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=0.9",y=X20),data=data2,position="identity")
g<-g+geom_violin(mapping=aes(x="r=1",y=X21),data=data2,position="identity")
print(g)
それでは、サンプル数を変えてみたらどうなるのだろうか。
冗長になってしまうので結果は省略するが、nを減らすごとに偏回帰係数の分散は大きくなる傾向にあった。
よって、多重共線性を回避するためにはnをできるだけ大きくするということが重要であると言えるかもしれない。
