一般線形回帰（OLS）でコスト分析をしたら係数がマイナスの項目が出てきた、工数をかけるほど費用が減っていく･･･と言う嘆きを見かけたのですが、説明変数の相関が高いと推定される係数の標準誤差が大きくなって推定量がおかしくなる、多重共線性による問題だと考えられます*1。
こういう問題を解決する素朴な方法として、主成分回帰（Principal Component Regression）で多重共線性を制御する技があるので紹介します。試してみると望んだ結果が得られるかも知れないので、学部生で卒論に使う回帰分析の係数が気に入らない人にp-hacking的にお勧めです。なお、べたべたとコードを書いていきますが、plsパッケージで手軽に実行でき、さらに同様の目的に使えるPLS回帰も試せるので、基本的にはplsパッケージを使う方が賢いです。
さて、まず解決すべき多重共線性をつくってみましょう。最後に比較するときに、真のモデルの切片項と係数はそれぞれ1であることに注意してください。
set.seed(20201115) n <- 150 # 十分なサンプルサイズ（e.g. 1500）にすると多重共線性が問題なくなる x1 <- runif(n) x2 <- x1 + rnorm(n, sd=0.1) y <- 1 + x1 + x2 + rnorm(n)
VIFは8.02675です。
計算結果を格納する行列をつくります。
beta <- matrix(0, 3, 2) se <- matrix(0, 3, 2) colnames(beta) <- colnames(se) <- c("OLS", "PCR") rownames(beta) <- rownames(se) <- c("Const.", "x1", "x2")
比較のためにOLSで回帰します。
lm_r <- lm(y ~ x1 + x2) beta[, "OLS"] <- coef(summary(lm_r))[, 1] se[, "OLS"] <- coef(summary(lm_r))[, 2]
主成分分析を行います。scale=TRUEを忘れないようにしましょう。
r_prcomp <- prcomp(~ x1 + x2, scale=TRUE)
寄与率を低い主成分を省略し、回帰を行います。今回は、summary(r_prcomp)してCumulative Proportionが0.9を超えているので、第1主成分だけで行います*2。回帰に用いる主成分の数と説明変数の数が一致すると、OLSと同じ予測力になります。
npc <- 1 # 回帰に用いる主成分の数 lm_prc_r <- lm(y ~ predict(r_prcomp)[,1:npc])
第1主成分を、例えば平均全体工数のような抽象的な概念に対応するとして分析を進めてもよいのですが、そのような概念をあまり説明したくない場合は、主成分の係数を、元の説明変数の係数に分解しましょう。第1主成分の回帰係数と第1主成分の固有ベクトル（説明変数を標準化した値の係数）から、元の説明変数の回帰係数を求めます*3。
beta["x1", "PCR"] <- sum(r_prcomp$rotation[1,][1:npc]*coef(lm_prc_r)[2:(npc + 1)]/sd(x1)) beta["x2", "PCR"] <- sum(r_prcomp$rotation[2,][1:npc]*coef(lm_prc_r)[2:(npc + 1)]/sd(x2)) beta["Const.", "PCR"] <- coef(lm_prc_r)[1] - beta["x1", "PCR"]*mean(x1) - beta["x2", "PCR"]*mean(x2)
続けて分散の加法公式を使って、上の計算手順をもとに、元の説明変数の擬似的な標準誤差を求めます。
se["x1", "PCR"] <- sqrt(sum((r_prcomp$rotation[1, ][1:npc]/sd(x1) * coef(summary(lm_prc_r))[, 2][2:(npc + 1)])^2)) se["x2", "PCR"] <- sqrt(sum((r_prcomp$rotation[2, ][1:npc]/sd(x2) * coef(summary(lm_prc_r))[, 2][2:(npc + 1)])^2)) se["Const.", "PCR"] <- sqrt( coef(summary(lm_prc_r))[, 2][1]^2 + sum((( r_prcomp$rotation[1,][1:npc]*mean(x1)/sd(x1) + r_prcomp$rotation[2,][1:npc]*mean(x2)/sd(x2) ) * coef(summary(lm_prc_r))[, 2][2:(npc + 1)])^2) )
t検定に使う場合の自由度は、標本サイズから、回帰に使った主成分の数と、切片項の分の1を引いたものになるのに注意してください。
OLSとPCRの推定量と残差平方和を比較しておきましょう。
# 推定量を比較する cat("\nestimated coeffcients:\n") print(beta) cat("\nstandard errors:\n") print(se) # 残差平方和を比較する cat( sprintf("\nOLSの残差平方和: %.3f\n主成分回帰の残差平方和: %.3f\n", sum((y - predict(lm_r))^2), sum((y - model.matrix(~ x1 + x2) %*% beta[, 2])^2) ) )
| 係数 | 真のモデル | OLS | PCR |
|---|---|---|---|
| Const. | 1 | 1.10562969 | 1.0806319 |
| x1 | 1 | -0.01407916 | 0.9303908 |
| x2 | 1 | 1.80575174 | 0.8960958 |
| 標準誤差 | OLS | PCR |
|---|---|---|
| Const. | 0.1775728 | 0.09022349 |
| x1 | 0.8366178 | 0.05424751 |
| x2 | 0.8057794 | 0.05224790 |
OLSの残差平方和は143.823、PCRの残差平方和は145.11と予測力は僅かに悪化しますが、係数を見るとPCRの推定量の方が真のモデルに近い事がわかります。
*1:さらにomitted variable biasで、見えない成分を代理してしまっている可能性もあります。例えば、常連客との取引では低コストになる仕組みがあり、常連客によく提供するサービス項目の工数であったりすると、常連客のコスト低減効果を代理することもありえます。
*2:省略する成分は、偶発的に生じた見かけ上の、もしくは想定外だが従属変数への影響が大きくない構造（e.g. 上述の常連客のコスト低減効果）になります。
*3:を主成分回帰の切片工、
を主成分回帰の回帰式の係数、
を主成分の説明変数
の係数、主成分回帰の誤差項を
とすると、主成分回帰は
と書けるので、式を
の項で分けて整理すると、OLSと同様の説明変数とその係数の式に直せます。
