多重共線性とサンプルサイズの関係

2023年11月26日2023年11月27日

はじめに

多重共線性、通常 “マルチコ" について、名前は良く聞くし、なんだかよくないものだということは知っているけど、それが起きる原理や結果の解釈をどうすればよいのかよくわからない、という方は多いのではないでしょうか。私もその一人です。

Nospareさんの記事を拝見した所、多重共線性の起きる理由・問題点として、

・変数間の相関が高い事により、係数の推定の際に逆行列の計算が不安定になって推定値の分散が大きくなってしまう。
・多重共線性はOLSE（最小二乗法）の推定精度が下がることが問題（＝推定値の標準誤差が大きくなることが問題、と解釈しています）であり、推定値の標準誤差が小さければ、多重共線性は気にしなくてよい。

と整理されていました。

2つ目の点について、「標準誤差が小さく、変数が有意になっていれば、係数をそのまま解釈してしまって良いということなのか…？」などという点がわからなかったこともあり、シミュレーション的に確認してみようと思います。

分析の流れとして、まずは5つの説明変数X₁~X₅の線形結合によって目的変数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 = X₁ + 2X₂ + 3X₃ + 4X₄ + 5X₅ + ε
という式ですね。εは誤差項です。

回帰分析の結果、推定される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を生成させる式は先ほどと同様ですが、X₅の作り方を先ほどと変えています。このデータでは、X₅がX₁とX₂の和で生成されていることから、X₁-X₅、X₂-X₅、の間にマルチコが発生しています。
途中で相関係数の確認も入れました。X₁-X₅・X₂-X₅の相関係数はそれぞれ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 ‘ ’ 1

X₁・X₂・X₅の係数の推定値が、大幅にずれてしまっていることがわかります。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

先ほどと比較すると、推定した係数が真値にかなり近くなっていることがわかりますね。先程は有意になっていなかったX₁の推定値は、有意になっていることが見て取れます。

それでも、係数の推定値は真の値から少しずれています。やはり、「標準誤差が小さく、変数が有意になっていれば、係数をそのまま解釈してしまって良い」と言い切ってしまうのは危険なようです。

サンプルサイズと二乗誤差の推移の可視化

ここまでで、サンプルサイズを増やすことにより、係数の推定値が真の値に近づくことを示すことができました。

最後に、多重共線性が発生しているデータに対してサンプルサイズを増やしていった際に、係数の真値と推定値の二乗誤差がどのように推移していくのか確認してみましょう。コードは下記です。

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 = X₁ + 2X₂ + 3X₃ + 4X₄ + 5X₅ + ε
上記の設定でデータを生成させています。

このデータに対し、
Y = β₁X₁ + β₂X₂ + β₃X₃ + β₄X₄ + β₅X₅ + ε
という式で推定したβ₁ ~ β₅の推定精度を、係数の真値（1~5）との差の二乗誤差で確認しています。

実行結果はこちら。

X₁-X₅・X₂-X₅のそれぞれにマルチコを発生させているので、β₁・β₂・β₅の推定精度が悪い（2乗誤差が大きい）傾向にあることが見て取れます。サンプルサイズが小さい間はこの誤差が顕著ですが、サンプルサイズを増やすことで、マルチコの影響を回避し、係数の推定値が真の値に近づいていることがわかりますね。

一方で、マルチコを発生させていないβ₃・β₄については、サンプルサイズが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()