多重共線性の仕組みとシミュレーション
多重共線性の確認
マーケティングリサーチで重回帰分析を用いる際、各説明変数の係数に関心が行くことは多いと思う。「これこれの説明変数は他に比べて影響度が強いですね、コイツを上げる(下げる)施策を考えましょう*1」などと言って施策の方向性を定めていくという具合である。
そこで必ず問題になるのがこの多重共線性で、一般に説明変数間に強い相関がある際に偏回帰係数の推定が不安定になるという状況を指す。この数理については『統計的方法のしくみ』の解説が抜群にわかりやすい。そこで多重共線性の仕組みについて、シミュレーションで確かめつつ記事に残しておくことにする。
- 作者: 永田靖
- 出版社/メーカー: 日科技連出版社
- 発売日: 1996/10/01
- メディア: 単行本
- 購入: 27人 クリック: 226回
- この商品を含むブログ (7件) を見る
推定の流れ
この記事では2つの説明変数で重回帰分析を行うことを考える。回帰モデルは
で表現される。最小二乗法では残差平方和を最小にするようパラメータを求めるので、残差を
としたとき、残差平方和
を最小にするを求める。各パラメータで偏微分して=0と置いて求めれば良いのでそのまま云々とやっていくと、行列を用いて以下のように表現できる。
ここではそれぞれの偏差平方和、はとの偏差積和、はそれぞれと、との偏差積和。あとは逆行列を左から掛けてやれば良いので、
ふぅ*2。これで偏回帰係数が求められた。
多重共線性のメカニズムと相関との関係
先の式を見ると行列にが掛かっており、これが0のとき(すなわち偏差平方和・積和行列に逆行列が存在しないとき)は解が無数に存在する(不定)か解が存在しない(不能)ということになる。この状況を多重共線性が存在するという。
では、このという状況と説明変数間の相関の関係がどうなっているかというと、
これにより、説明変数間の相関が±1という状況で多重共線性が発生するということがわかった*3。またが±1に極めて近ければが0に極めて近いということになるので、解の変化幅が大きくなり安定しないということになる*4。
Rでシミュレーション
ここまでで多重共線性の仕組みが分かったので、実際回帰式の推定にどのように影響するのかをシミュレーションで確かめてみる。シミュレーション内容はこんな感じ。
- N=1,000で目的変数と1つ目の説明変数の正規乱数を作る(今回はで固定しておく*5)
- 2つ目の説明変数についてはを0, 0.05, 0.1, ..., 1と0.05ずつ刻み、各々の値に対し100個のセットを作る
- 毎にパラメータがどのようにぶれるのか確認
#シミュレーション用データの作成 set.seed(221) #目的変数と説明変数1の作成 y <- rnorm(1000, mean = 0, sd = 1) #目的変数 e <- rnorm(1000, mean = 0, sd = 1) #説明変数作成のために必要なデータ p=0.4 #yとx1の相関係数を固定する x1 <- p * y + sqrt(1 - p ^ 2) * e #yとの相関が0.4の説明変数x1を作成 #説明変数2を作成し、シミュレーションを実施 r = seq(0, 1, by = 0.05) #r_12 = 0, 0.05, ..., 1で設定(21段階) b0 <- matrix(nrow = 100, ncol = 21) #β0を格納する箱。各相関係数に対し100個ずつ b1 <- matrix(nrow = 100, ncol = 21) #β1を格納する箱。各相関係数に対し100個ずつ b2 <- matrix(nrow = 100, ncol = 21) #β2を格納する箱。各相関係数に対し100個ずつ for(i in 1:21){ stock <- matrix(nrow = 1000, ncol = 100) for(j in 1:100){ e2 = rnorm(1000, mean = 0, sd = 1) x2 <- r[i] * x1 + sqrt(1 - r[i] ^ 2) * e2 stock[, j] <- x2 dat <- as.data.frame(cbind(y, x1, x2 = stock[, j])) res <- lm(y ~ ., data = dat) b0[j, i] <- res$coefficients[1] b1[j, i] <- res$coefficients[2] b2[j, i] <- res$coefficients[3] } }
これでについてそれぞれ、100(ある下でのデータセット数)×21(の段階数)の値が得られた。例えばの箱にはこんな感じで値が入っている。
> head(b1) r=0 r=0.05 r=0.1 r=0.15 r=0.2 r=0.25 r=0.3 r=0.35 r=0.4 r=0.45 r=0.5 1 0.4306679 0.4325078 0.4283519 0.4317270 0.4361706 0.4260468 0.4184288 0.4508224 0.4204306 0.4261973 0.4178097 2 0.4322390 0.4346633 0.4311487 0.4308968 0.4269444 0.4343129 0.4436880 0.4380596 0.4259308 0.4228499 0.4535279 3 0.4313652 0.4297881 0.4267221 0.4289261 0.4415001 0.4260835 0.4101387 0.4224972 0.4297717 0.4291471 0.4383937 4 0.4304963 0.4330273 0.4283483 0.4204201 0.4253410 0.4492216 0.4335296 0.4235167 0.4709619 0.4294702 0.4160259 5 0.4322628 0.4335394 0.4328325 0.4331517 0.4354849 0.4317349 0.4314967 0.4338696 0.4191118 0.4356418 0.3976430 6 0.4334116 0.4313089 0.4324587 0.4292935 0.4345609 0.4290567 0.4176773 0.4438420 0.4313572 0.4149576 0.4371784 r=0.55 r=0.6 r=0.65 r=0.7 r=0.75 r=0.8 r=0.85 r=0.9 r=0.95 r=1 1 0.4111970 0.4093910 0.4339163 0.4184000 0.3734337 0.4756565 0.3922978 0.4544657 0.4121455 0.4322047 2 0.4391379 0.4516133 0.3873439 0.4027791 0.4100673 0.4379951 0.3401908 0.3585950 0.3741827 0.4322047 3 0.4076027 0.4311014 0.4450212 0.4278456 0.4481529 0.4148706 0.4109423 0.5119509 0.4214167 0.4322047 4 0.4324386 0.4308051 0.4209028 0.4485870 0.4614509 0.3932755 0.4713389 0.5392622 0.4512087 0.4322047 5 0.4360566 0.4087040 0.4727114 0.4316713 0.4444006 0.4488128 0.4205936 0.4277047 0.4467345 0.4322047 6 0.4354719 0.4178111 0.4086015 0.4483683 0.4363689 0.3716729 0.3623556 0.4071001 0.6429823 0.4322047 > nrow(b1) [1] 100
どのような結果になっているか箱ひげ図で確認してみると、が1に近づくほどブレ幅が大きくなることがわかる。
#b1の偏回帰係数の図示 b1 <- as.data.frame(b1) colnames(b1) <- colnames_r b1_long <- b1 %>% tidyr::gather(key = r, value = b1) ggplot(b1_long, aes(x = r, y = b1)) + geom_boxplot() + ggtitle("説明変数間の相関の強弱に伴うb1のバラツキ") #b2の偏回帰係数の図示 b2 <- as.data.frame(b2) colnames(b2) <- colnames_r b2_long <- b2 %>% tidyr::gather(key = r, value = b2) ggplot(b2_long, aes(x = r, y = b2)) + geom_boxplot() + ggtitle("説明変数間の相関の強弱に伴うb2のバラツキ")
●の結果図示
●の結果図示
ちなみにについて同じように図示するとこうなる。
●の結果図示
これは、については
を解くことで求まるが、についてはどのデータセットでも同じ値であること、そのうえでの組み合わせ(足し合わせ)が
> head(b1 + b2) r=0 r=0.05 r=0.1 r=0.15 r=0.2 r=0.25 r=0.3 1 0.4062166 0.4154052 0.3839790 0.4027693 0.4471804 0.4732325 0.4681899 2 0.4664639 0.4493094 0.4712963 0.4394332 0.4249584 0.4497007 0.4273785 3 0.4267857 0.4824163 0.4557214 0.4181341 0.4777092 0.4289483 0.4333954 4 0.4379343 0.4109607 0.4335442 0.4378974 0.4731251 0.4555382 0.4145050 5 0.4744823 0.4237386 0.4455729 0.4360844 0.4437070 0.4832751 0.4168798 6 0.4561966 0.4234361 0.4466039 0.3999932 0.4673920 0.4392135 0.4219600 r=0.35 r=0.4 r=0.45 r=0.5 r=0.55 r=0.6 r=0.65 1 0.4214690 0.3901941 0.4546184 0.4290851 0.4504682 0.4186020 0.4276473 2 0.4738807 0.4074830 0.4196819 0.4178540 0.4350354 0.4318700 0.4474259 3 0.4628155 0.4458085 0.4428629 0.4416364 0.4580332 0.4516505 0.4254354 4 0.4202259 0.4383358 0.4085167 0.4249698 0.4341524 0.4211597 0.4240621 5 0.4645901 0.4214547 0.4467117 0.4362892 0.4431544 0.4008554 0.4301498 6 0.4132844 0.4117565 0.4340837 0.4367685 0.4224682 0.4045233 0.4415290 r=0.7 r=0.75 r=0.8 r=0.85 r=0.9 r=0.95 r=1 1 0.4308617 0.4302392 0.4172367 0.4228712 0.4474329 0.4418713 NA 2 0.4246184 0.4306662 0.4290404 0.4310337 0.4285018 0.4379655 NA 3 0.4430312 0.4253368 0.4116594 0.4210978 0.4446660 0.4302211 NA 4 0.4359602 0.4359433 0.4402329 0.4316962 0.4229854 0.4388992 NA 5 0.4174786 0.4173395 0.4187569 0.4235048 0.4437448 0.4364285 NA 6 0.4327213 0.4251063 0.4441639 0.4309610 0.4294942 0.4313048 NA
とrによらずほぼ一定であることから理解できる。
終わりに
ここまで多重共線性の仕組みとその影響を確認してきた。今回は値自体に意味はないため結果の実際的な解釈は行っていないが、多重共線性は係数の大きさや符号を不安定なものにする危険があることがよくわかった。オーソドックスな回帰でも多重共線性をはじめいくつか落とし穴や確認すべき点があるので、結果に基づいた意思決定を行う際は十分注意して分析していこうと思う次第です。
中心極限定理に頼ってイノベーター理論をそれっぽく使う
イノベーター理論はなぜ正規分布なのか?
マーケターなら誰でも知っているだろう古典的な理論に、イノベーター理論がある。セグメントの設定とその切り方(分布)について実務的、概念的にはわかりやすく納得感があるのでよく出てくる考え方だが、統計をかじる人間からすると「なんで正規分布なん?」という疑問が出てくるだろう。実際この理論を根拠にターゲティングの提案をしたとき「なんで正規分布なん?」と思いいろいろググったことがあり、その時見つけた記事が以下。
よく参考にさせてもらっているブログで、特に仮説検定のシリーズ
A/Bテストのガイドライン:仮説検定はいらない(Request for Comments|ご意見求む) - 廿TT
はとても面白く、また参考になった。
で今回も参考にさせていただいたのだが、やはりとても面白く、続きの記事に書かれている実務マーケターのフィードバックまで鮮明にイメージできた(笑)
そういうわけで基本的には記事に沿ってマーケティングを考えていけば良いと思うのだけれど、普及した理論や概念を修正したりそれらの範疇外のアプローチをとろうとすると予想以上に混乱を招き、現場の停滞を生むことが往々にしてある*1。
ということで本記事では観点を変えて、イノベーター理論に「当てはまるような」アプローチをとってみる*2。ビジネスが停滞するくらいなら、「馴染みのある理論からも実際の観点からもそんなに間違ってない」提案をしてコトを進めようという試みである*3。
着想
ということで思いついたのが、「中心極限定理」を利用するというアプローチ。統計学の中核をなす定理でその内容は、「x1, x2, …, xnが独立に同一の分布にしたがうときnが十分大きいとき、その平均は近似的に正規分布に従う」というものである。これを基礎に「経済数学の直感的方法 確率・統計編」では、
この世の中には多種多様な確率分布があって、その分布パターンは当然どれも違う形をしている。(中略)ところが実はそれらを全部集めて混ぜると、それらの様々な確率分布を合成したものは、驚いたことに結局はきれいな正規分布になってしまうというのである。
と説明しており*4、ピンときた。以下はたぶんマーケターの方にはウケが良いだろうアプローチ。
Rでシミュレーション
まず、イノベーター理論の5セグメントの順序に比例して強く(or弱く)なる変数を定義する。今回は以下の6つの変数を定義。
①先進性:新しい商品は進んで採用する
②流行感度:流行に敏感である
③情報収集:情報収集を自ら積極的に行う
④波及性:他者への影響力を持っている
⑤自分本位:他人の評価よりも自分の判断基準に従う
⑥情動性:費用対効果を重視せずに商品を購入する傾向がある
これらについて順序尺度でアンケートとってスコアとみなして足せばイノベーター理論になるよ、という話。やってみる。
まずは6カテゴリそれぞれについて「あなた自身に当てはまる度合いはどの程度ですか?」と5件法でデータを得るケースを想定。各6カテゴリが5値のいずれかをとる確率分布に従うよう、適当に確率分布を作成する。
#イノベーター理論のためのダミーデータ作成 set.seed(123) #N=1,000で4桁のユニークIDを作成 sampleID <- sample(1000:9999, 1000, replace = F) #5値の確率分布を6つ作成 tekito <- sample(20:70, 30) tekito <- matrix(tekito, c(6,5)) p <- prop.table(tekito, margin = 1) #pの各行を6カテゴリそれぞれの5段階値の出る確率として与える ----------------------------------------------- > p [,1] [,2] [,3] [,4] [,5] [1,] 0.1358025 0.2222222 0.2798354 0.1440329 0.2181070 [2,] 0.2462312 0.1608040 0.1959799 0.1407035 0.2562814 [3,] 0.1764706 0.1437908 0.1307190 0.2222222 0.3267974 [4,] 0.2752294 0.1743119 0.2155963 0.1972477 0.1376147 [5,] 0.2468619 0.2552301 0.2384937 0.1506276 0.1087866 [6,] 0.1518519 0.2592593 0.2407407 0.1148148 0.2333333 -----------------------------------------------
次に上記で作った確率分布をもとにサイズ1,000の回答データを作成し、アンケートローデータを得る。
#5段階尺度のアンケート回答データを適当に作成 for (i in 1:6){ x <- sample(1:5, 1000, replace = T, prob = p[i, ] ) assign(paste("category", i, sep = ""), x) } score <- as.data.frame(cbind(sampleID, category1, category2, category3, category4, category5, category6)) ----------------------------------------------- > head(score) sampleID category1 category2 category3 category4 category5 category6 1 3588 3 3 1 5 1 4 2 8093 5 3 4 3 3 3 3 4679 2 3 4 1 3 5 4 8944 3 5 5 4 2 5 5 9460 5 5 2 3 2 4 6 1409 5 2 4 5 4 5 -----------------------------------------------
ここで各カテゴリについてヒストグラムを書いて確認すると、こんな感じ。
#6カテゴリーのヒストグラムは par(mfrow = c(2, 3)) hist(score$category1, breaks = seq(0.5, 5.5, 1), ylim = c(0, 350), xlab = "自分には当てはまらない ⇔ 自分には当てはまる", col = "slateblue2", labels = T, main = "先進性") hist(score$category2, breaks = seq(0.5, 5.5, 1), ylim = c(0, 350), xlab = "自分には当てはまらない ⇔ 自分には当てはまる", col = "slateblue2", labels = T, main = "流行感度") hist(score$category3, breaks = seq(0.5, 5.5, 1), ylim = c(0, 350), xlab = "自分には当てはまらない ⇔ 自分には当てはまる", col = "slateblue2", labels = T, main = "情報収集") hist(score$category4, breaks = seq(0.5, 5.5, 1), ylim = c(0, 350), xlab = "自分には当てはまらない ⇔ 自分には当てはまる", col = "slateblue2", labels = T, main = "波及性") hist(score$category5, breaks = seq(0.5, 5.5, 1), ylim = c(0, 350), xlab = "自分には当てはまらない ⇔ 自分には当てはまる", col = "slateblue2", labels = T, main = "自分本位") hist(score$category6, breaks = seq(0.5, 5.5, 1), ylim = c(0, 350), xlab = "自分には当てはまらない ⇔ 自分には当てはまる", col = "slateblue2", labels = T, main = "情動性")
これらの値は順序尺度なので扱い的に問題ありそうだが、やむを得ずイノベーターとしての強さを表現するスコアの要素(=間隔尺度とみなす)とし、足し合わせた値をイノベータースコアとして定める。
#6カテゴリの回答を足し合わせた「イノベータースコア」を作成 score$innovator.score <- apply(score[2:7], 1, sum) ----------------------------------------------- > head(score) sampleID category1 category2 category3 category4 category5 category6 innovator.score 1 3588 3 3 1 5 1 4 17 2 8093 5 3 4 3 3 3 21 3 4679 2 3 4 1 3 5 18 4 8944 3 5 5 4 2 5 24 5 9460 5 5 2 3 2 4 21 6 1409 5 2 4 5 4 5 25 ----------------------------------------------- #イノベータスコアのヒストグラム par(mfrow = c(1, 1)) hist(score$innovator.score, breaks = seq(6.5, 30.5, 1), ylim = c(0, 150), xlab = "ラガード ⇔ イノベーター", col = "#ff00ff40", labels = T, main = "イノベータースコア") #Q-Qプロット qqnorm(score$innovator.score)
分布を見てみるとこんな感じで、正規分布っぽい。正規性を確認するためのQ-Qプロットを描いてみると、割と綺麗な直線になっていて正規分布とみなしてよさそう。
せっかくなのでもう一つ、正規性を検定するためのシャピロ・ウィルク検定をやってみると、
shapiro.test(score$innovator.score) ----------------------------------------------- Shapiro-Wilk normality test data: score$innovator.score W = 0.99149, p-value = 1.59e-05 -----------------------------------------------
え…?高度に有意な結果。シャピロ・ウィルク検定の帰無仮説は「対象が正規分布に従う」なので、これを棄却すると「対象が正規分布に従わない」と結論づけることになる。
これではすっきりしないので、今回のイノベータースコアの平均と標準偏差に従う正規乱数を同じ1,000サイズで出し、それの結果を見比べてみることにする。
score.mean <- mean(score$innovator.score) #イノベータースコアの平均 score.sd <- sd(score$innovator.score) #イノベータースコアの標準偏差 #イノベータースコアで算出した平均と標準偏差に従う正規分布をプロット n <- rnorm(1000, mean = score.mean, sd = score.sd) hist(n, breaks = seq(6.5, 30.5, 1), ylim = c(0, 150), xlab = "ラガード ⇔ イノベーター", col = "#0000ff40", labels = T, main = "イノベータースコア")
#ヒストグラムを重ね合わせる par(mfrow = c(1, 1)) hist(score$innovator.score, breaks = seq(6.5, 30.5, 1), ylim = c(0, 150), col = "#ff00ff40") hist(n, breaks = seq(6.5, 30.5, 1), ylim = c(0, 150), col = "#0000ff40", add = T) #正規性の検定 shapiro.test(n) ----------------------------------------------- Shapiro-Wilk normality test data: n W = 0.99928, p-value = 0.9745 -----------------------------------------------
分布に関してはよく重なっている(図中紫の部分)にも関わらず、データnの方は検定の結果p=0.97で帰無仮説を棄却できない(=正規分布に従っていないとはいえない≒正規分布に従うとみなす)。ということでもうひと手間、正規乱数の小数点以下を丸めたデータn2を作り、同じことをやってみると、
#nの小数点を丸める n2 <- round(n, digits = 0) hist(n2, breaks = seq(6.5, 30.5, 1), ylim = c(0, 150), col = "#00FF7F44", add = T) #正規性の検定 shapiro.test(n2) ----------------------------------------------- Shapiro-Wilk normality test data: n2 W = 0.99151, p-value = 1.616e-05 -----------------------------------------------
分布はnと同一なのに検定結果は有意という結果になる。これはどうも質的データを量的データとみなして進めてしまったことが原因みたい。
シャピロ・ウィルク検定
シャピロ・ウィルクの検定における検定統計量はwiki(シャピロ–ウィルク検定 - Wikipedia)をざっと見た感じ、順序統計量ってのが出てきてるからこれに25値しかとらないイノベータースコアを当てはめてしまったのが良くない気がする。
終わりに
ということで、4尺度に関して勉強するとき必ず出てくる「順序尺度において四則演算は意味がない」というのはこういうときに効いてくる話なんだと思った。シャピロ・ウィルクの検定を材料にしたこのあたりの数理の検討はまた別の記事で頑張ってみよう。
一方でもともとの主題だったイノベータースコアの妥当化?というか作り方については、発想としてはそんなに間違ってないと思うし実際正規分布に従う分布を描けることを確認した*5。
ビジネスの理論が先走って後戻りできないときに、ひとまずの妥協案として今回みたいなやり方でセグメントを定義しておいて、その後集めたデータで修正を図っていくという進め方も時にはありなんじゃないかな、と思った次第です。
このブログについて
このブログは統計学をほぼ独学で学び、日々の業務においてそれを使っている(あるいは使うチャンスを常にうかがっている)人間による、統計学の学習内容とそれに付随するよしなしごとのノートです。何分大部分において独学なので、間違いなど多々あるかもしれませんが、そこはご了承ください。また出来ましたら、そのような間違いについてはぜひご指摘・ご指導いただけると大変うれしいです。以下、今のところ思っていることのサマリ。
■ブログの内容
主には、統計学系の教科書やRに関する技術書で学んだことをまとめていく予定です。
またその知識を日常業務でどう活かすか、あるいは社会に(ほんの少しだけでも)価値提供するにはどうすればよいかなどをはじめとして、統計学(データ分析といっても良いか)と人・社会を結び付けたときに考えたこと、感じたことなども書いていければと思っています。
なのでデータを扱う人間として、というような観点で伝記とか小説とかも持ち出したりするかも。
■ブログの目的
自分にとっては学習内容をまとめることで知識の定着を図ることが目的です。
読んでいただいた方にはほんの少しでも暇つぶしになったり、また願わくはご自身の学習に少しだけ貢献することができればとっても嬉しいです。
なのでコメント、ご指摘、叱咤激励、愚痴、なんでもござれ。
■著者の関心
究極的にはデータを分析するとは、一体何なのか、を探求することです。
統計学は人の意思決定にどう貢献するのか?その貢献の仕方とは?こういったような「人が考えるときのツールとしてのデータ/分析の正体」について私は確信のある答えを持っていません。日々勉強し、使ってみて、これについて考えていきたいと思っています。
■著者について
ITベンダーでデータアナリストとして4年勤め、この春からはコンサルティングファームに環境を移しました。やることは変わらず、データを分析して、少しでも人のより良い意思決定に貢献することです。
学部ではエネルギー系のお勉強、そこからもっとソフトなところで技術を活かすことを考えたいと思うようになり、大学院ではまちづくりの分野へ。ヒトとモノ、ビジネスとアカデミア、理系と文系、そういったパッと見対立概念にあるような二つの狭間にいるのが何となく性に合う人間です。蝙蝠にならないようもっと勉強せねば。。。