線形回帰モデルに従うデータセットがあるものの、分散が均一とは見なし難いとします。
誤差項の大きさが何かと相関していると見做せれば、一般化最小二乗法（GLS）で有効推定量を得るのがよい方法になりますが、誤差項の構造は明確ではないことが多いです。
この場合は、係数の推定量はOLSのものを用います。不均一分散があっても一致性と不偏性があるからです。しかし、\(E[\epsilon {{}^t\!\epsilon} \mid X]\)を\(\sigma^2 I\)に置き換えるところができません。これは均一分散を意味するからです。
そこで、
\[ \text{VAR}[ b \mid X ] = ({}^t\!XX)^{-1}{{}^t\!X} E[\epsilon {{}^t\!\epsilon} \mid X] X ({}^t\!XX)^{-1} \]
の\(E[\epsilon {{}^t\!\epsilon} \mid X]\)を他で置き換えることを考えます。もっとも単純なのが残差\(e\)の二乗を対角成分にとる正方行列で置き換えることです。
\[ \approx ({}^t\!XX)^{-1}{{}^t\!X} \begin{pmatrix} e_1 e_1 & 0 & \cdots & 0 \\ 0 & e_2 e_2 & \cdots & 0 \\ 0 & 0 & \ddots & 0 \\ 0 & 0 & \cdots & e_n e_n \end{pmatrix} X ({}^t\!XX)^{-1} \]
ホワイトの方法と呼ばれており、これに沿った標準誤差はホワイト頑強標準誤差もしくはHC₀と呼ばれます。出発点となる不均一分散に頑強な標準誤差（Heteroskedasticity-Consistent Standard Errors）です。
HC₀は漸近的に（真の値に）一致する分散共分散行列になるのですが、下方推定バイアスが入る事が知られています。均一分散のときにホワイトの方法を用いるとします。OLSの分散の推定量の導出を思い出すと、
\[ E[e{}^t\!e] = E[(M \epsilon){}^t\!(M \epsilon)] = \sigma^2 (I-X({}^t\!XX)^{-1}({}^t\!X)) \]
射影行列（projection matrix）と呼ばれる\(X({}^t\!XX)^{-1}{}^t\!X\)の対角成分は梃子比（leverage）\(h_{ii}\)ですから、\(h_{ii}\)だけバイアスが入る事がわかります。漸近的には一致します。OLSの不偏推定量の導出の議論から、\(h_{ii}<1\)で\(\sum_{i=1}^{n} h_{ii} = k\)となり、係数の数\(k\)は一定ですから、サンプルサイズ\(n\)が増えると\(h_{ii}\)は小さくなっていくからです。しかし、とくに小標本ではバイアスが出ます。また、ダミー変数などにも弱いようです。
HC₀の改良として様々な方法が提案されています。HC₁は標本サイズに応じて補正をかけますが、統計学的に十分な正当化はありません。HC₂は\(1-h_{ii}\)で割ることでバイアスを補正しようとしますが、不均一分散のときの性質はあまり良くないそうです。HC₃は実用上の性質はよいようですが、過剰推定バイアスが入ることが知られています。HC₁もHC₂もHC₃も自由度\(n-k\)のt分布に従わないので、p値や信頼区間の計算には調整した自由度が必要になります。他にHC₄やブートストラッピング手法などありますが、決定版と言える手法はないそうです1。
しかし、完璧ではないとは言え、近年は不均一分散に頑強な標準誤差を求められることも多くなりました。計算方法は抑えておきましょう。
1 データセット
分散（標準偏差）が不均一なデータセットをつくります。
2 不均一分散に頑強な標準誤差の推定
2.2 HC₂の推定（自由度調整なし）と係数の検定
sandwichパッケージでHC₂標準誤差は出ます。また、HC₂だけではなく、HC₀からHC₅までサポートされています。
library(sandwich)
vcov_hc2 <- vcovHC(r_lm, type = "HC2")
matrix(c(sqrt(diag(vcov(r_lm))), sqrt(diag(vcov_hc2))), 3, 2,
dimnames = list(rownames(vcov_hc2), c("OLS", "HC2"))) OLS HC2
(Intercept) 1.7077845 1.4845115
x 0.5591708 0.5106786
z 0.4989273 0.4212183lmtestパッケージを使えば、sandwichで出した分散共分散行列を用いて係数の検定と信頼区間の計算ができます。
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.26962 1.48451 -0.1816 0.8572374
x 2.15981 0.51068 4.2293 0.0002407 ***
z -3.26220 0.42122 -7.7447 2.499e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1実践ではこの方法が使われる事が多いです。
しかし、t分布にあわせるのに必要な調整した自由度が分からないので、厳密にはp値と信頼区間が不適切になります2。サンプルサイズを\(n\)、説明変数の数を\(k\)と置いたときのOLSの自由度\(n-k\)でt検定をかけていますが、この自由度は大きすぎます。
2.3 HC₂の推定（自由度調整あり）と係数の検定
HC₂であればlmで推定した後にdfadjustパッケージを用いる方法があります。調整された自由度の計算もしてくれますし、堅牢な手法の数値計算で処理されています。
Coefficients:
Estimate HC1 se HC2 se Adj. se df p-value
(Intercept) -0.2696173 1.4301396 1.4845115 1.6784187 10.42182 8.593647e-01
x 2.1598056 0.4997741 0.5106786 0.5550900 15.08424 7.201512e-04
z -3.2621983 0.4106650 0.4212183 0.4659995 12.54645 3.973517e-06どの標準誤差を用いてp値の計算をしているのか不安になる表が出で、さらにヘルプでは明記されていないのですが、ソースコードを確認するとHC₂を用いています。実際、
(Intercept) x z
8.593647e-01 7.201512e-04 3.973517e-06 で、print(r_dfa)で表示されるp値を再現できます。
2.3.1 信頼区間
dfadjustパッケージにはconfintがないので、つくっておきます。
confint.dfadjustSE <- function(object, parm, level = 0.95, ...){
lwr_upr <- c((1 - level)/2, 1 - (1 - level)/2)
beta <- object$coefficients[, "Estimate"]
se <- object$coefficients[, "HC2 se"]
df <- object$coefficients[, "df"]
index <- 1:length(beta)
if(!missing(parm)) index <- index[names(beta) %in% parm]
confint_da <- t(sapply(index, \(i) beta[i] + se[i] * qt(lwr_upr, df[i])))
rownames(confint_da) <- names(beta)[index]
colnames(confint_da) <- sprintf("%.1f %%", 100*lwr_upr)
confint_da
}2.4 ブートストラッピングによる信頼区間の計算
分散共分散行列を計算するのではなく、ブートストラッピングで係数のサンプルを作成し、その信頼区間を求める方法もあります。バイアスが少ないWild Bootstrappingを用いて推定してみましょう。
library(lmboot)
r_wb <- wild.boot(y ~ x + z, B = 30000, data = df01, seed = NULL, bootDistn = "normal")
class(r_wb) <- c(class(r_wb), "wild.boot") # lmbootパッケージが戻り値にクラス名を指定しないので、後で総称関数confintを使うためにクラス名をつける
t(apply(r_wb$bootEstParam, 2, quantile, c(0.025, 0.975))) # 95%信頼区間 2.5% 97.5%
(Intercept) -2.938482 2.399354
x 1.237316 3.083495
z -4.021931 -2.494218fix(wild.boot)でソースコードを読めますが、残差に乱数で変化をつけてベクトル\(e_{bootstrapping}\)をつくり、説明変数の行列を\(X\)として、\(\beta_{OLS} +
({}^t\!XX)^{-1}{}^t\!Xe_{bootstrapping}\)のサンプルをつくる作業をB回繰り返しています。
サンプルから信頼区間を計算するのには、上のようにパーセンタイル法を使うのが簡便ですが、Bootstrappingに入るバイアスを減らしたBCa法をつかいたい場合はbayestestRパッケージの助けを借りると楽です。
confint.wild.boot <- function(object, parm, level = 0.95, ...){
index <- 1:ncol(object$bootEstParam)
if(!missing(parm)) index <- index[names(beta) %in% parm]
M <- t(sapply(index, \(i) bayestestR::bci(object$bootEstParam[, i], ci = level)[-1]))
lwr_upr <- c((1 - level)/2, 1 - (1 - level)/2)
rownames(M) <- colnames(object$bootEstParam)[index]
colnames(M) <- sprintf("%.1f %%", 100*lwr_upr)
M
}3 推定方法ごとの信頼区間の比較
ブートストラッピング手法があるので、標準誤差ではなく信頼区間で比較します。
2.5 % 97.5 %
(Intercept) -3.773702 3.234467
x 1.012482 3.307129
z -4.285913 -2.238484 2.5 % 97.5 %
(Intercept) -3.315583 2.776349
x 1.111980 3.207631
z -4.126467 -2.397930 2.5 % 97.5 %
(Intercept) -3.559257 3.020023
x 1.071849 3.247762
z -4.175541 -2.348856 2.5 % 97.5 %
(Intercept) -2.941929 2.398598
x 1.240602 3.087872
z -4.016813 -2.490674ばらばらですね。
詳しくは末石 (2025)「ミクロ計量分析講義資料 スライド2: 不均一分散とクラスターに頑健な統計的推測」やImbens and Kolesár (2016)を参照してください。↩︎
dfオプションがあるので自由度は渡せます。↩︎