生活習慣病のコホートでHbA1cを説明するとき、BMIと腹囲を同じ回帰モデルへ入れることがあります。どちらも体格を反映するため、同じ集団では強く相関することがありますが、この状態で個々の係数を見ると、モデル全体としてはHbA1cをある程度説明できているのに、BMIと腹囲の95%信頼区間は広く、推定値の符号まで標本ごとに変わることがあります。
ここでVIFが10を超えたという理由で片方を削除すると、標準誤差は小さくなるかもしれませんが、同時に、係数が表している比較も変わってしまいます。多重共線性で分けて考えたいのは、係数の偏り、係数の精度、モデルの当てはめ、そして研究で知りたい量の4つです。乱数で生成した300例と2,000回の反復生成を使い、どこに影響が出るのかを確認します。
BMIと腹囲を同時に入れると何が起きるか
重回帰のBMI係数は、腹囲と年齢を同じ値にしたうえでBMIが1kg/m^2高いことに対応するHbA1cの平均値の差です。この比較条件については単回帰分析と重回帰分析の違いで扱いました。回帰係数そのものの位置づけは回帰分析の目的につながります。
BMIと腹囲が強く相関していると、腹囲と年齢を調整したあとに残るBMIのばらつきが小さくなります。BMIそのものに十分な個人差があっても、その大部分を腹囲から予測できるなら、「同じ腹囲なのにBMIだけが違う」という情報はわずかしか残りません。多重共線性で係数の精度が落ちるのは、この残った情報だけを使ってBMI固有の係数を推定するためです。
なお、完全に線形従属している場合と、強く相関している場合は分けて考えます。ある説明変数が他の説明変数から完全に再現できれば、通常の最小二乗法では個々の係数を一意に決められませんが、一方で、今回のように相関が0.953でも完全一致ではなければ係数は計算でき、その代わり標準誤差が大きくなります。
300例のデータでは、腹囲を平均88cm、標準偏差10cmで生成し、BMIは腹囲と強く関連するように作りました。HbA1cの生成式にはBMI1kg/m^2あたり0.040%、腹囲1cmあたり0.012%、年齢1歳あたり0.012%の係数を置いています。
set.seed(42)
library(car)
generate_cohort <- function(n, collinear = TRUE) {
age <- rnorm(n, mean = 60, sd = 10)
waist <- rnorm(n, mean = 88, sd = 10)
if (collinear) bmi <- 3.0 + 0.22 * waist + rnorm(n, sd = 0.70)
if (!collinear) bmi <- 25.0 + rnorm(n, sd = sqrt(0.22^2 * 100 + 0.70^2))
hba1c <- 4.2 + 0.040 * bmi + 0.012 * waist + 0.012 * age + rnorm(n, sd = 0.60)
data.frame(hba1c, bmi, waist, age)
}
cohort <- generate_cohort(300)
fit_full <- lm(hba1c ~ bmi + waist + age, data = cohort)
vif(fit_full)#> bmi waist age
#> 10.946948 10.953857 1.005291BMIの係数は0.0808、標準誤差は0.0547、95%信頼区間は-0.0269から0.1885でした。腹囲の係数は0.0029で、95%信頼区間は-0.0214から0.0273です。一方、BMIや腹囲とはほぼ独立に生成した年齢のVIFは1.005で、係数0.0129、標準誤差0.0038となりました。
この1回のデータでは、強く重なったBMIと腹囲で区間が広くなり、年齢では同じ現象がほとんど見られません。ここから言えるのは、多重共線性がモデル内の全係数を一様に不安定にするのではなく、他の説明変数からどれだけ予測されるかによって係数ごとに影響が異なるということです。
VIFは標準誤差の拡大を表している
分散拡大係数（VIF）は、ある説明変数を残りの説明変数でどこまで説明できるかから計算します。$j$番目の説明変数を他の説明変数へ回帰したときの決定係数を$R_j^2$とすると、VIFは次の式です。
$$\mathrm{VIF}_j=\frac{1}{1-R_j^2}\tag{1}$$
式(1)では、他の説明変数で説明されずに残った割合$1-R_j^2$が小さいほどVIFが大きくなります。VIFはアウトカムを使わず、説明変数同士の関係だけから決まります。
今回、BMIを腹囲と年齢で説明したときの$R_j^2$は約0.909でした。したがってVIFは約10.95になります。これは「BMIの90.9%が腹囲と年齢によって決定される」という意味ではなく、通常の線形回帰で得られた決定係数が0.909だったという記述にとどめます。
通常の最小二乗法で切片を含み、標準的な分散推定を使う場合、係数の標準誤差は次の形に分けられます。$\hat{\sigma}$を残差標準偏差、$s_j$を説明変数$j$の標本標準偏差、$n$を標本サイズとします。
$$\mathrm{SE}(\hat{\beta}_j)=\frac{\hat{\sigma}}{s_j\sqrt{n-1}}\sqrt{\mathrm{VIF}_j}\tag{2}$$
式(2)では、残差の大きさ、説明変数自身のばらつき、標本サイズが同じなら、共線性による標準誤差の倍率は$\sqrt{\mathrm{VIF}_j}$です。今回のBMIでは$\sqrt{10.95}=3.31$なので、共線性がない同じ条件と比べた標準誤差0.0165が0.0547まで拡大します。
vif_bmi <- unname(vif(fit_full)["bmi"])
sigma_hat <- summary(fit_full)[["sigma"]]
se_ideal <- sigma_hat / (sd(cohort[["bmi"]]) * sqrt(nrow(cohort) - 1))
c(VIF = vif_bmi, se_ideal = se_ideal, se_actual = se_ideal * sqrt(vif_bmi))
r_grid <- seq(0, 0.99, by = 0.01)
plot(r_grid, 1 / (1 - r_grid^2), type = "l", log = "y",
xlab = "Correlation", ylab = "Variance inflation factor")
abline(h = 10, lty = 2)#> VIF se_ideal se_actual
#> 10.946948 0.01654249 0.05473278
ここで、VIFが大きいことと係数に偏りが生じることは分けます。通常の線形モデルの仮定が成り立っているなら、説明変数同士が相関していることだけを理由に最小二乗推定量の期待値がずれることはありません。今回問題になっているのは主に分散です。測定誤差、モデルの誤指定、交絡などがあれば別の問題が加わりますが、それらをVIFの影響としてまとめると原因を取り違えます。
VIF=10も境界線ではありません。Marquardt（1970）は10を大きなVIFの目安として扱いましたが、O’Brien（2007）は、係数の分散がVIFだけでは決まらないことを理由に、固定的なしきい値から変数削除へ進む方法に注意を促しています。VIFが10を超えたかどうかより、その研究で必要な係数の信頼区間がどれだけ広がっているかを見る方が適切です。
VIFだけで変数を落とすと比較条件が変わる
腹囲をモデルから外すと、BMIの標準誤差は小さくなります。ただし、そこで推定しているBMI係数は、腹囲を入れたモデルのBMI係数と同じ量ではありません。
今回のデータで、腹囲をBMIと年齢へ回帰したときのBMI係数を$\hat{\delta}$とします。腹囲を含む完全モデルの係数を$\hat{\beta}_{\mathrm{bmi}}$、$\hat{\beta}_{\mathrm{waist}}$とすると、腹囲を外したモデルのBMI係数$\tilde{\beta}_{\mathrm{bmi}}$には次の関係があります。
$$\tilde{\beta}_{\mathrm{bmi}}=\hat{\beta}_{\mathrm{bmi}}+\hat{\delta}\hat{\beta}_{\mathrm{waist}}\tag{3}$$
式(3)では、腹囲を外したBMI係数に、腹囲との関連と腹囲の係数が混ざります。今回の$\hat{\delta}$は4.22なので、BMIが1kg/m^2高い人では、年齢を同じ値にしたとき腹囲も平均4.22cm高いというデータ上の関係が反映されます。
delta <- coef(lm(waist ~ bmi + age, data = cohort))["bmi"]
contrast_vec <- c(0, 1, delta, 0)
combo_est <- as.numeric(contrast_vec %*% coef(fit_full))
combo_se <- sqrt(as.numeric(t(contrast_vec) %*% vcov(fit_full) %*% contrast_vec))
c(delta = unname(delta), estimate = combo_est, se = combo_se,
reduced = unname(coef(lm(hba1c ~ bmi + age, data = cohort))["bmi"]))#> delta estimate se reduced
#> 4.21672846 0.09310672 0.01654260 0.09310672腹囲を外したBMI係数0.0931は、腹囲を固定してBMIだけを1kg/m^2変えた比較ではありません。BMIが1kg/m^2異なるときに腹囲もデータ上の関係に沿って異なる比較が混ざった値です。標準誤差が小さくなったから、BMI固有の係数を精密に推定できるようになった、と読むと対象が入れ替わります。
この点は、交絡調整を目的とするモデルでは特に重要です。ある変数を交絡調整に必要だと事前に定めているなら、VIFが大きいことだけを理由に削除すると、目的としていた調整そのものを変えてしまいます。Schistermanら（2009）が示しているように、変数を調整するかどうかは、その変数が因果構造のどこにあるかによって、偏りを減らす場合も増やす場合もあります。VIFにはその因果的な役割を判定する機能はありません。
したがって、「VIFが10を超えたら片方を削除する」という規則は適切とは言えず、BMIと腹囲のどちらが研究上必要なのか、両方を別々の係数として解釈したいのか、それとも体格を表す一つの指標として扱いたいのかを先に決めるべきです。
2,000回繰り返すと係数だけが大きく散らばる
1回のデータだけを見ると、BMIの係数0.0808が生成時に置いた0.040より大きかったことに目が向きます。そこで同じ条件から300例を2,000回生成し直し、係数がどこを中心に散らばるかを確認します。
sim <- matrix(NA_real_, nrow = 2000, ncol = 4)
for (i in seq_len(2000)) {
d_high <- generate_cohort(300)
b <- coef(lm(hba1c ~ bmi + waist + age, data = d_high))
delta_i <- coef(lm(waist ~ bmi + age, data = d_high))["bmi"]
d_low <- generate_cohort(300, collinear = FALSE)
b_low <- coef(lm(hba1c ~ bmi + waist + age, data = d_low))
sim[i, ] <- c(b["bmi"], b["waist"], b_low["bmi"], b["bmi"] + delta_i * b["waist"])
}
rbind(mean = apply(sim[, c(1, 3, 4)], 2, mean),
sd = apply(sim[, c(1, 3, 4)], 2, sd))#> [,1] [,2] [,3]
#> mean 0.03974788 0.04000130 0.08975577
#> sd 0.05010088 0.01505049 0.01504796相関の強いデータでのBMI係数は平均0.0397で、生成式に置いた0.040の近くにあります。一方、反復ごとの標準偏差は0.0501でした。BMIと腹囲の相関をなくし、BMI自身の分散はほぼ同じにしたデータでは、BMI係数の平均は0.0400、標準偏差は0.0151です。
つまり今回の生成条件では、共線性によってBMI係数の中心が一定方向へ移ったのではなく、標本ごとの散らばりが約3.3倍になっています。1回目の0.0808も、その広い標本分布から出てきた一つの値として読む方が状況に合っています。
plot(density(sim[, 1]), lwd = 2, xlab = "Estimated BMI coefficient",
ylab = "Density", main = "Sampling distributions")
lines(density(sim[, 3]), lwd = 2, lty = 2)
plot(sim[, 1], sim[, 2], pch = 16, cex = 0.4,
xlab = "BMI coefficient", ylab = "Waist coefficient")
abline(v = 0.040, h = 0.012, lty = 2)
BMIと腹囲の係数は、2,000回の反復で強い負の相関を示しました。BMIの係数が上へ動いた標本では腹囲の係数が下へ動き、その逆も起きます。BMI係数は21.2%の反復で0を下回りました。生成式ではBMI係数を0.040と置いているので、符号の反転は「BMIがHbA1cを下げる」という別の関係を発見したことを意味しません。2つの似た説明変数へ効果を分配する位置が、標本ごとに大きく動いています。
式(3)で作った組み合わせの標準偏差は0.0150でした。ただし、この結果から「共線性があれば、この組み合わせを報告すればよい」と一般化はできません。ここでは$\delta$を同じデータから求めているので、何を比較する量なのかがデータの相関構造に依存しています。研究上意味のある複合指標を使うのであれば、BMIと腹囲をどうまとめるかは解析結果を見る前に定める方が解釈しやすくなります。
予測を目的とする場合には、Harrell（2015）が扱う変数クラスタリングや主成分など、アウトカムを使わずに関連の強い変数群を縮約する方法も候補になります。
予測では係数とは別の不安定さを見る
多重共線性について、「個々の係数は不安定でも予測には影響しない」と説明されることがありますが、ここは条件を付ける必要があります。
BMIと腹囲を、同じ情報を保つ別の線形結合へ置き換えただけなら、モデルが張る空間は変わりません。そのため、同じデータに対する当てはめ値も変わりません。これは共線性に強いという性質ではなく、同じ線形モデルを別の座標で表現しているためです。
一方、強く相関したBMIと腹囲のデータでは、観測値は説明変数空間の細い帯の中に集まります。その帯に沿った組み合わせなら予測が比較的安定していても、帯から離れた組み合わせでは事情が変わります。たとえばBMIが高いのに腹囲だけが極端に小さい人を予測すると、学習データではほとんど観測されていない方向へモデルを延長することになります。
したがって、予測モデルでVIFが高い場合には、個々の係数の標準誤差だけを見て終わりにはできません。新しい対象者の説明変数の組み合わせが学習データの範囲に入っているか、予測性能が内部検証や外部検証で保たれているかを別に確認します。係数推定の問題と予測性能の問題は、同じ「共線性」という言葉で一つにまとめない方が判断しやすくなります。
研究ではVIFより変数の役割を先に決める
Vatchevaら（2016）は疫学研究における多重共線性を扱い、BMIや腹囲を含むデータも使って、相関やVIFなどの診断が結果の解釈に関わることを示しています。大切なのは、VIFを計算したあと自動的に変数を削ることではなく、どの係数がどの程度不安定になっているのかを把握することです。
私なら、VIFが大きかったときには、その変数をモデルへ入れた目的から確認します。研究上の主要な曝露なら、信頼区間が広いという結果そのものを受け入れる場合がありますが、交絡調整に必要な変数なら、VIFだけを理由に削除しません。同じ構成概念を測った代替指標が複数あるだけなら、どの指標を使うかを研究目的から決めるか、事前に定めた縮約方法を使います。予測が目的なら、個別係数の解釈より検証データでの予測性能を優先します。
また、カテゴリー変数やスプラインのように1つの項が複数の回帰係数を持つ場合、通常の1自由度のVIFをそのまま当てはめることはできません。FoxとMonette（1992）は、複数の係数からなる項へ分散拡大の考え方を拡張した一般化VIF（GVIF）を提案しています。carパッケージでは複数自由度の項にGVIFが表示され、項の自由度を考えた$\mathrm{GVIF}^{1/(2\mathrm{df})}$も確認できます。
区間が広いからといって、ロバスト標準誤差へ切り替えれば共線性が解消することもありません。ロバスト標準誤差が扱うのは分散不均一などによる分散推定の問題であり、説明変数同士に独立した情報が少ないことはそのまま残ります。施設内相関に対するクラスター頑健な標準誤差も同様です。
多重共線性が見つかったときに必要なのは、VIFを小さくする操作ではなく、知りたい係数を現在のデータがどの程度の精度で切り分けられるのかを確認し、その係数自体が研究目的に必要なのかを判断します。VIFは、その判断を始めるための診断値として使う方が役に立ちます。
実際の研究では、説明変数の役割、因果構造、測定指標の重複、必要な推定精度、予測か効果推定かという目的まで合わせてモデルを決めることになります。Dr.データサイエンスでは、こうした変数選択と回帰モデルの組み立てを含む解析方針のご相談を承っています。
参考文献
- Marquardt DW. Generalized Inverses, Ridge Regression, Biased Linear Estimation, and Nonlinear Estimation. Technometrics. 1970;12(3):591-612.
- Fox J, Monette G. Generalized Collinearity Diagnostics. Journal of the American Statistical Association. 1992;87(417):178-183.
- O’Brien RM. A Caution Regarding Rules of Thumb for Variance Inflation Factors. Quality & Quantity. 2007;41(5):673-690.
- Schisterman EF, Cole SR, Platt RW. Overadjustment bias and unnecessary adjustment in epidemiologic studies. Epidemiology. 2009;20(4):488-495.
- Vatcheva KP, Lee M, McCormick JB, Rahbar MH. Multicollinearity in Regression Analyses Conducted in Epidemiologic Studies. Epidemiology (Sunnyvale). 2016;6(2):227.
- Harrell FE Jr. Regression Modeling Strategies: With Applications to Linear Models, Logistic and Ordinal Regression, and Survival Analysis. 2nd ed. Springer International Publishing; 2015.
