多重共線性とサンプルサイズの関係
はじめに
多重共線性、通常 “マルチコ" について、名前は良く聞くし、なんだかよくないものだということは知っているけど、それが起きる原理や結果の解釈をどうすればよいのかよくわからない、という方は多いのではないでしょうか。私もその一人です。
Nospareさんの記事を拝見した所、多重共線性の起きる理由・問題点として、
・変数間の相関が高い事により、係数の推定の際に逆行列の計算が不安定になって推定値の分散が大きくなってしまう。
・多重共線性はOLSE（最小二乗法）の推定精度が下がることが問題（＝推定値の標準誤差が大きくなることが問題、と解釈しています）であり、推定値の標準誤差が小さければ、多重共線性は気にしなくてよい。
と整理されていました。
2つ目の点について、「標準誤差が小さく、変数が有意になっていれば、係数をそのまま解釈してしまって良いということなのか…？」などという点がわからなかったこともあり、シミュレーション的に確認してみようと思います。
分析の流れとして、まずは5つの説明変数X1~X5の線形結合によって目的変数Yが生成されるという仮定でシミュレーションデータを作成します。その後、回帰分析を行い、生成の際の係数をうまく推定できるかを確認するというフローです。
マルチコの発生していない回帰分析
まずは手始めに、通常の回帰分析を試してみましょう。今回の検証ではすべての分析設計において、説明変数の数は5個、サンプルサイズは200と2,000の2パターンを検証します。
マルチコの設定はなしで、シミュレーションデータを生成しましょう。
set.seed(0)
n <- 200
p <- 5
# シミュレーションデータ生成
X1 <- matrix(rnorm(n * p), ncol = p)
beta <- seq_len(5)
y1 <- X1 %*% beta + rnorm(n)
df1 <- data.frame(y1, X1)目的変数でYの生成のためにXを使いますが、その係数は1~5で設定しています。
Y = X1 + 2X2 + 3X3 + 4X4 + 5X5 + ε
という式ですね。εは誤差項です。
回帰分析の結果、推定されるXの係数がそれぞれ上記の値になれば成功です。確認してみましょう。
model1 <- lm(y1 ~ ., data = df1)
summary(model1)結果は下記です。
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.06937 0.07163 -0.969 0.334
X1 0.85091 0.07879 10.799 <2e-16 ***
X2 2.05427 0.07212 28.483 <2e-16 ***
X3 2.96640 0.07016 42.279 <2e-16 ***
X4 4.12497 0.07714 53.473 <2e-16 ***
X5 4.90851 0.06684 73.432 <2e-16 ***推定した係数の結果は “Estimate Std"の列ですが、おおむね1~5でうまく推定できていると言ってしまって良さそうです。p値からも、係数が有意であると判断されます。
n=2,000に変更して上から再度実行するとどうなるでしょうか。
結果は下記です。
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.01362 0.02229 0.611 0.541
X1 1.03618 0.02310 44.863 <2e-16 ***
X2 1.97913 0.02216 89.302 <2e-16 ***
X3 3.01615 0.02252 133.908 <2e-16 ***
X4 4.00351 0.02227 179.808 <2e-16 ***
X5 4.96349 0.02191 226.535 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1推定した係数の値がより真値に近くなっていること、2列目の標準誤差が小さくなって推定の信頼性が高まっていることが見て取れますね。
ここまでで、マルチコが発生していない回帰分析において、サンプルサイズを増やすことで、最小二乗法の推定精度を上げられる（係数の標準誤差を小さくできる）こと、真値との誤差を小さくできること、が確認できました。
マルチコのあるデータで回帰分析
では次に、マルチコが発生しているデータで回帰分析を行ってみましょう。分析設計は先ほどと同じですが、一部マルチコを発生させてみます。
set.seed(0)
n <- 200
p <- 5
# シミュレーションデータ生成
X2 <- matrix(rnorm(n * p), ncol = p)
X2[,5] <- X2[,1] + X2[,2] + rnorm(n, sd = 0.1)
beta <- seq_len(5)
# 相関係数の確認
cor(X2)
y2 <- X2 %*% beta + rnorm(n)
df2 <- data.frame(y2, X2)Yを生成させる式は先ほどと同様ですが、X5の作り方を先ほどと変えています。このデータでは、X5がX1とX2の和で生成されていることから、X1-X5、X2-X5、の間にマルチコが発生しています。
途中で相関係数の確認も入れました。X1-X5・X2-X5の相関係数はそれぞれ0.69と0.73と、しっかりとマルチコが発生していますね。このデータで回帰分析を実行してみましょう。
model2 <- lm(y2 ~ ., data = df2)
summary(model2)結果は下記です。
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.04524 0.06928 -0.653 0.515
X1 0.56244 0.68478 0.821 0.412
X2 1.58705 0.69770 2.275 0.024 *
X3 3.07531 0.06771 45.421 < 2e-16 ***
X4 4.09608 0.07466 54.860 < 2e-16 ***
X5 5.40666 0.69069 7.828 3.14e-13 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1X1・X2・X5の係数の推定値が、大幅にずれてしまっていることがわかります。X2は有意水準5%、X5は有意水準0.1%、でそれぞれ有意となっていますが、係数の推定はうまくいっているとはいえません。
この時点で、最初の疑問にあった、「標準誤差が小さく、変数が有意になっていれば、係数をそのまま解釈してしまって良い」と言ってしまうのは危険であるように思えます。
では次に、サンプルサイズを増やすことで、マルチコを回避して係数の推定をうまくできるのかどうかを確認してみましょう。コードは上の"n <- 200″の部分を適当に変えて上から実行し直すだけです。
下記が、n=2000にした際の推定結果です。
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.00996 0.02169 0.459 0.646
X1 1.20806 0.21239 5.688 1.48e-08 ***
X2 2.20981 0.21162 10.442 < 2e-16 ***
X3 2.99762 0.02133 140.546 < 2e-16 ***
X4 4.00987 0.02217 180.866 < 2e-16 ***
X5 4.80950 0.21085 22.810 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1先ほどと比較すると、推定した係数が真値にかなり近くなっていることがわかりますね。先程は有意になっていなかったX1の推定値は、有意になっていることが見て取れます。
それでも、係数の推定値は真の値から少しずれています。やはり、「標準誤差が小さく、変数が有意になっていれば、係数をそのまま解釈してしまって良い」と言い切ってしまうのは危険なようです。
サンプルサイズと二乗誤差の推移の可視化
ここまでで、サンプルサイズを増やすことにより、係数の推定値が真の値に近づくことを示すことができました。
最後に、多重共線性が発生しているデータに対してサンプルサイズを増やしていった際に、係数の真値と推定値の二乗誤差がどのように推移していくのか確認してみましょう。コードは下記です。
library(ggplot2)
library(tidyr)
library(dplyr)
beta_true <- seq_len(5)
sample_sizes <- c(10, 100, 200, 300, 400)
errors <- matrix(NA, nrow = length(sample_sizes), ncol = length(beta_true))
for (i in 1:length(sample_sizes)) {
set.seed(0)
n <- sample_sizes[i]
p <- length(beta_true)
X3 <- matrix(rnorm(n * p), ncol = p)
X3[,5] <- X3[,1] + X3[,2] + rnorm(n, sd = 0.1)
y3 <- X3 %*% beta_true + rnorm(n)
df3 <- data.frame(y3, X3)
model <- lm(y3 ~ ., data = df3)
# 誤差計算
beta_estimated <- coef(model)[-1] # 切片を除外
errors[i, ] <- abs(beta_estimated - beta_true)^2
}
colnames(errors) <- c("Beta1", "Beta2", "Beta3", "Beta4", "Beta5")
errors_df <- data.frame(sample_size = sample_sizes, errors) %>%
pivot_longer(cols = -sample_size, names_to = "Variable", values_to = "Value")
ggplot(errors_df, aes(x = as.factor(sample_size), y = Value)) +
geom_bar(stat = "identity", position = position_dodge(), aes(fill=as.factor(sample_size))) +
facet_wrap(. ~ Variable, nrow = 1) +
scale_fill_manual(values = c("lightblue", "blue", "darkblue", "royalblue", "steelblue")) +
theme_minimal() +
labs(x = "", y = "2乗誤差", fill="サンプルサイズ") +
theme(text = element_text(size = 10))
データの生成過程はこれまでと同様、
Y = X1 + 2X2 + 3X3 + 4X4 + 5X5 + ε
上記の設定でデータを生成させています。
このデータに対し、
Y = β1X1 + β2X2 + β3X3 + β4X4 + β5X5 + ε
という式で推定したβ1 ~ β5の推定精度を、係数の真値（1~5）との差の二乗誤差で確認しています。
実行結果はこちら。
X1-X5・X2-X5のそれぞれにマルチコを発生させているので、β1・β2・β5の推定精度が悪い（2乗誤差が大きい）傾向にあることが見て取れます。サンプルサイズが小さい間はこの誤差が顕著ですが、サンプルサイズを増やすことで、マルチコの影響を回避し、係数の推定値が真の値に近づいていることがわかりますね。
一方で、マルチコを発生させていないβ3・β4については、サンプルサイズが10と小さい場合でも、誤差がかなり小さくなっており、推定が正確であることが見て取れます。マルチコが発生していない状況では、サンプルサイズが小さめでも大きな問題なく回帰分析ができそうです。
最後に、上のグラフを1つの標本と見立て、シミュレーションデータの乱数seedを変えることで標本を増やし、真値との誤差の動きを可視化してみましょう。
beta_true <- seq_len(5)
sample_sizes <- c(10, 100, 200, 300, 400)
p <- length(beta_true)
errors = c()
for (i in 1:length(sample_sizes)) {
error_temp = c()
n <- sample_sizes[i]
for (seed in seq_len(30)) {
set.seed(seed)
X3 <- matrix(rnorm(n * p), ncol = p)
X3[,5] <- X3[,1] + X3[,2] + rnorm(n, sd = 0.1)
y3 <- X3 %*% beta_true + rnorm(n)
df3 <- data.frame(y3, X3)
model <- lm(y3 ~ ., data = df3)
# 誤差計算
beta_estimated <- coef(model)[-1] # 切片を除外
error_temp <- rbind(error_temp, abs(beta_estimated - beta_true) / beta_true)
}
error_temp = error_temp |> as.data.frame() |> mutate(n = n)
errors = rbind(errors, error_temp)
}
errors_summary = errors |> group_by(n) |> summarise(across(everything(), list(mean=mean, sd=sd)))
colnames(errors_summary) <- c("sample_size", "Beta1_mean", "Beta1_sd", "Beta2_mean", "Beta2_sd", "Beta3_mean", "Beta3_sd", "Beta4_mean", "Beta4_sd", "Beta5_mean", "Beta5_sd")
# 平均値と標準偏差の列を区別し、縦長にする
errors_summary_long <- errors_summary |>
pivot_longer(cols = -sample_size, names_to = c("variable", "stat"), names_pattern = "(.*)_(.*)")
# 平均値と標準偏差を別々のデータフレームに分割してから、標準誤差に計算して結合
errors_summary_mean <- errors_summary_long |> filter(stat == "mean")
errors_summary_sd <- errors_summary_long |> filter(stat == "sd")
errors_summary_mean$se = errors_summary_sd$value / sqrt(errors_summary_sd$sample_size)
ggplot(errors_summary_mean, aes(x=as.factor(sample_size), y=value)) +
facet_wrap(~variable, nrow=1) +
scale_fill_manual(values = c("lightblue", "blue", "darkblue", "royalblue", "steelblue")) +
geom_point() +
geom_errorbar(aes(ymin = value - 2*se, ymax = value + 2*se), position=position_dodge(.9), width=0.4, show.legend=FALSE) +
geom_bar(aes(color=NULL, fill=as.factor(sample_size)), alpha=0.5, position='dodge', stat='identity') +
theme_minimal(base_family="HiraKakuPro-W3") +
labs(x = "", y = "絶対パーセント誤差（信頼区間付き）", fill="サンプルサイズ") +
theme_minimal()
サンプルサイズが小さいときは推定値の信頼区間もブレブレですが、サンプルサイズを増やすことで、マルチコが発生させている変数についても、信頼区間を狭めて精緻な分析ができるようになっていますね。
まとめと考察
実務でよく問題になるマルチコの問題について、サンプルサイズを増やすことで、係数の推定精度を上げられることを示しました。
ちなみに、サンプルサイズを増やして係数の推定精度が高くなっている場合でも、VIFの値はサンプルサイズが小さい場合とほとんど変化していないので、VIFの解釈として、「VIFが高い → 係数推定精度が悪くなる」という直接的な解釈にはつながらないこともわかりますね。
サンプルサイズが足りないことは、マーケティングの実務でも頻繁に問題になります。不用意にサンプルサイズを増やした分析により有意性を主張することは問題ですが、回帰分析係数の際の回帰係数の推定を真値に近づける目的でいうと、サンプルサイズを増やすことは手段の1つになりそうです。
サンプルサイズはどの程度を基準とすれば良いのか、リッジ回帰でどこまでマルチコを回避できるか、という点は、今後の研究課題とします。