多施設で集めた測定値に回帰モデルを当てはめると、治療群の係数は同じままなのに、標準誤差の計算方法を変えるだけで信頼区間が大きく動く場合があります。残差のばらつきが対象者によって異なるデータでも、通常の回帰出力だけでは推定の不確実性を表しきれません。
こうした場面で使うロバスト標準誤差は、回帰係数の分散を計算し直す方法です。ただし、不等分散に対する補正と、同じ施設や同じ患者から得た観測値の相関に対する補正は別物で、時系列データではさらに時間方向の依存が加わります。何に対してロバストなのかを決めると、使う関数だけでなく、解析単位や信頼区間の作り方まで決まります。
この記事では、不等分散だけを扱います。標準偏差と標準誤差の違いから始めて、通常の回帰が置いている仮定、HC0からHC3までの補正の違い、そして400人のランダムに生成したデータで、係数を変えずに標準誤差だけが動く様子を確かめます。
標準偏差と標準誤差の違い
患者ごとの測定値がどれくらい散らばっているかを表す標準偏差と、回帰係数をどれくらい精密に推定できているかを表す標準誤差は、対象が異なります。同じ母集団から研究を繰り返した際に、推定値が研究ごとに動く程度を表すのが標準誤差で、実際の解析では、その大きさを手元の標本から推定しています。
たとえば、施設内の患者数を増やせば測定値は増えますが、同じ施設の患者が共通した要因の影響を受けている場合、その全てを互いに独立した情報として扱うことはできません。患者間のばらつきが大きいかどうかと、治療群の平均差をどれくらい精密に推定できるかは、患者数だけでは結び付かないためです。
また、95%信頼区間は、その計算手順を同じ条件で繰り返した際に、真の値を含む割合が95%となることを目指した区間です。得られた区間に真の値が入っている確率を、そのまま95%と読むものではありません。そして、ロバスト標準誤差を使った区間にも、独立した情報の数や大標本近似などの条件が残ります。Mansourniaらは、医療・疫学研究で混同されやすい分散と相関の問題を、具体例から説明しています。
回帰係数と分散の役割
線形回帰では、回帰係数を求める計算と、その係数の分散を推定する計算を分けて考えられます。ここでは、対象者数を$n$、切片を含む推定係数の数を$p$とし、説明変数を並べた$n$行$p$列の行列を$\mathbf{X}$、アウトカムのベクトルを$\mathbf{Y}$、そして行列の各行に対応する説明変数ベクトルを、切片の1を含めて$\mathbf{x}_i$と表します。
通常最小二乗法、OLSによる係数は、$\mathbf{X}$が列フルランクなら次の形になります。
$$\widehat{\boldsymbol{\beta}} =(\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{Y} \tag{1}$$
式(1)は、観測値とモデルの予測値の差の二乗和が最も小さくなる回帰係数を求めています。
同じ対象者、同じ説明変数、同じ重みでOLSを当てはめ、その後に標準誤差の計算方法だけを変更しても、式(1)の係数は変わりません。したがって、ロバスト標準誤差を使った後で係数まで動いていたら、解析対象者の除外、重み、モデル式、あるいは推定方法も変わっていないかを確認します。最小二乗法の推定の仕組みと、係数を推定した後の分散計算は別の処理です。
条件付き平均が$\mathbf{X}\boldsymbol{\beta}$に等しい線形モデルで、誤差ベクトルを$\boldsymbol{\varepsilon}$、その条件付き分散共分散行列を$\boldsymbol{\Omega}$とすると、係数の分散は次の式で表されます。
$$\operatorname{Var}(\widehat{\boldsymbol{\beta}}\mid\mathbf{X}) =(\mathbf{X}^{\top}\mathbf{X})^{-1} \mathbf{X}^{\top}\boldsymbol{\Omega}\mathbf{X} (\mathbf{X}^{\top}\mathbf{X})^{-1}, \qquad \boldsymbol{\Omega} =\operatorname{Var}(\boldsymbol{\varepsilon}\mid\mathbf{X}) \tag{2}$$
式(2)は、対象者ごとの誤差の分散と、対象者同士の誤差の共分散が、回帰係数の不確実性へ伝わることを表しています。
独立で等分散なら、$\boldsymbol{\Omega}$は対角成分が全て$\sigma^2$、それ以外が0の行列になります。この条件を使う通常の分散推定量は、残差$\widehat{\varepsilon}_i$から求めた共通の分散を、次のように組み込みます。
$$\widehat{\mathbf{V}}_{\mathrm{OLS}} =\widehat{\sigma}^{2}(\mathbf{X}^{\top}\mathbf{X})^{-1}, \qquad \widehat{\sigma}^{2} =\frac{\sum_{i=1}^{n}\widehat{\varepsilon}_i^{\,2}}{n-p} \tag{3}$$
式(3)は、残差から求めた1つの分散を全ての対象者に共通して使い、係数の分散を計算しています。この形が成り立つ前提は、回帰分析の仮定にまとめてあります。
係数ごとの標準誤差は、この行列の対角成分の平方根です。患者によって誤差の分散が異なったり、同じ施設の患者同士で誤差が相関したりすると、式(3)で省略した部分を、そのまま無視できなくなります。
不等分散に対するHC0
説明変数の値が大きくなるにつれて残差が広がるデータでは、全員に共通の誤差分散を置くことが難しくなります。この状態の見つけ方と対処は、不均一分散性で扱っています。ただし、条件付き平均のモデルが正しく、説明変数と誤差の関係について必要な条件が満たされていれば、不等分散があるだけでOLS係数が不偏でなくなるわけではありません。まず影響を受けるのは、通常の標準誤差と、それに基づく推論です。
Whiteの分散推定量として知られるHC0では、各対象者の残差の二乗を使います。HCはheteroskedasticity-consistentの略で、不等分散の下でも、大標本で分散を整合的に推定するという意味です。
$$\widehat{\mathbf{V}}_{\mathrm{HC0}} =(\mathbf{X}^{\top}\mathbf{X})^{-1} \left\{ \sum_{i=1}^{n} \mathbf{x}_i\mathbf{x}_i^{\top}\widehat{\varepsilon}_i^{\,2} \right\} (\mathbf{X}^{\top}\mathbf{X})^{-1} \tag{4}$$
式(4)は、共通の誤差分散を使う代わりに、対象者ごとの残差の大きさを説明変数と組み合わせて、回帰係数の分散を推定しています。
残差の二乗が、その対象者の真の誤差分散を正確に復元する、という説明ではありません。1人から得た残差は1つしかなく、ばらつきます。それでも、独立した対象者について寄与を足し合わせ、推定係数の数や極端な観測値などに関する条件が満たされると、係数の分散全体を推定できます。
また、式(4)には、異なる対象者の残差を掛け合わせた項がありません。そのため、患者内や施設内の相関まで扱う式として、そのまま使うことはできません。不等分散と相関は、ここで分けておきます。
HC0からHC3までの補正
OLS残差は、真の誤差をそのまま観測したものではありません。同じデータで回帰係数を推定しているため、残差には当てはめによる縮小が入り、その影響は説明変数の配置によって異なります。
この配置を表す量がレバレッジです。ハット行列を$\mathbf{H}$、その$i$番目の対角成分を$h_{ii}$とすると、次の形になります。
$$\mathbf{H} =\mathbf{X}(\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}, \qquad h_{ii} =\mathbf{x}_i^{\top}(\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{x}_i \tag{5}$$
式(5)の$h_{ii}$は、その対象者の説明変数の組み合わせが、回帰の当てはめにどれくらい強く関わるかを表しています。
アウトカムが大きい対象者を指す量ではなく、説明変数側の位置から決まる点に注意が必要です。極端な年齢や、ほかの対象者ではほとんど見られない説明変数の組み合わせがあると、高いレバレッジになる場合があります。
HC1、HC2、HC3では、式(4)の中央で使う残差二乗を、次の$\omega_i$へ置き換えます。
$$\omega_i= \begin{cases} \widehat{\varepsilon}_i^{\,2}, & \mathrm{HC0},\\[2pt] \dfrac{n}{n-p}\widehat{\varepsilon}_i^{\,2}, & \mathrm{HC1},\\[6pt] \dfrac{\widehat{\varepsilon}_i^{\,2}}{1-h_{ii}}, & \mathrm{HC2},\\[6pt] \dfrac{\widehat{\varepsilon}_i^{\,2}}{(1-h_{ii})^{2}}, & \mathrm{HC3}. \end{cases} \tag{6}$$
式(6)は、HC1では全員に共通の自由度補正を加え、HC2とHC3では対象者ごとのレバレッジに応じて残差二乗を補正することを表しています。
HC3は、高いレバレッジの観測値に対してHC2より強く補正します。ただし、この説明はHC0との比較についてのもので、通常のOLS標準誤差よりHC3が必ず大きくなるという意味ではありません。通常の分散推定量とHC型では、各観測値のばらつきを係数へ反映する仕方そのものが異なるためです。
MacKinnonとWhiteは小標本での補正を比較し、LongとErvinもHC型の有限標本での挙動を調べています。sandwichパッケージのvcovHC()はHC3を既定値にしていますが、解析コードでは補正の種類を明記した方が、ほかのソフトウェアとの照合が容易です。
不等分散を持つ400人のデータ
ここからは、ランダムに生成したデータセットを使います。臨床研究の測定値を想定した数値ですが、実在する患者の記録でも、特定の治療に関するエビデンスでもありません。severityを連続的な基準時点の指標とし、その値が大きくなるほど測定値のばらつきが増えるように作ります。
$$Y_i=100+3S_i+\exp(0.8S_i)Z_i, \qquad S_i\sim\operatorname{Uniform}(-2,2), \qquad Z_i\sim N(0,1) \tag{7}$$
式(7)は、severityに相当する$S_i$が1増えると平均測定値が3増え、誤差の標準偏差は$\exp(0.8S_i)$に応じて変わるデータを生成しています。
対象者同士は独立で、$S_i$と$Z_i$も独立とし、平均の関係は直線のままにして、不等分散と平均モデルの誤指定を混ぜずに確認します。使用するのはRの標準機能とsandwichパッケージで、以下のコードは上から順に同じセッションで実行します。確認環境はR4.4.2、sandwich3.1-1です。
set.seed(42)
library(sandwich)
severity <- runif(400, min = -2, max = 2)
measurement <- 100 + 3 * severity + rnorm(400, sd = exp(0.8 * severity))
independent_fit <- lm(measurement ~ severity)
independent_se <- sqrt(c(OLS = vcov(independent_fit)[2, 2], HC0 = vcovHC(independent_fit, type = "HC0")[2, 2], HC1 = vcovHC(independent_fit, type = "HC1")[2, 2], HC2 = vcovHC(independent_fit, type = "HC2")[2, 2], HC3 = vcovHC(independent_fit, type = "HC3")[2, 2]))
print(round(cbind(Estimate = coef(independent_fit)[2], SE = independent_se), digits = 3))#> Estimate SE
#> OLS 3.093 0.086
#> HC0 3.093 0.112
#> HC1 3.093 0.112
#> HC2 3.093 0.112
#> HC3 3.093 0.113傾きの推定値は全て3.093で、通常の標準誤差は0.086、HC3では0.113となりました。変更したのは係数の計算ではなく分散の計算なので、推定値は揃ったままです。また、HC0、HC1、HC2は小数第3位で0.112に丸められていますが、丸める前の値は異なります。
この例で標準誤差が広がったことを、HC3が常に保守的になるという一般則へ広げることはできません。そして、真の傾き3に近い係数が得られたことと、その係数の信頼区間が繰り返し標本で目標どおりに振る舞うことも別の確認です。
png("robust-heteroskedasticity.png", width = 1100, height = 800, res = 200, type = "cairo")
plot(severity, y = residuals(independent_fit), pch = 16, col = "#0072B280", xlab = "Baseline severity", ylab = "OLS residual", main = "Independent observations, unequal variances")
abline(h = 0, lty = 2)
dev.off()
図1で右側の残差が広がったままなのは、補正が失敗したためではありません。ロバスト標準誤差は測定値も残差も書き換えず、そのばらつきを使って係数の分散を計算し直しています。分散を安定させるための変換や、観測値に重みを付けて係数を推定し直す方法とは、役割が異なります。
HC3を使う時の判断
残差プロットに扇形が見えた際に、直ちに対数変換へ進むか、それとも元の単位で平均差を推定するかは、研究の問いによって変わります。測定値の元の単位で平均の関係を説明したく、平均モデルもその目的に合っているなら、OLS係数とHC3を組み合わせる方針が候補になります。
私なら、独立した対象者の通常の線形回帰では、等分散性の検定のp値だけで標準誤差の種類を切り替えるより、平均モデル、レバレッジ、解析対象者の構造を先に確認し、分散計算を解析計画に書きます。等分散性の検定が有意でなかった場合も、等分散が成立したと確認できるわけではなく、標本が小さい時には検出できない逸脱が残るためです。
また、HC3を指定していても、係数の数が多く、少数の対象者が推定を強く左右する場合には、大標本近似が十分に働かないことがあります。HC4やHC5もレバレッジを意識した補正ですが、数字が大きい型ほど常に優れているという順序ではありません。どの補正を使ったかだけでなく、何を推定し、どれくらい独立した情報が残っているかを確認します。
補正の型より先に決めるもの
ここまでで動かしたのは、分散の計算方法だけです。係数3.093は最後まで変わらず、標準誤差が0.086から0.113へ動きました。
この差は、対象者ごとに残差のばらつきが違うという一点から出ています。HC0からHC3の選択は、その事実を認めた後の細部で、先に決まっているのは何を推定したいかと、どの対象者が独立だと考えるかです。
次の記事では、ロバスト標準誤差が直さない問題を確認したうえで、施設内で相関するデータへ進みます。同じ400人でも、依存の種類が変われば必要な補正も変わります。
参考文献
White H. A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity. Econometrica. 1980;48(4):817-838.
MacKinnon JG, White H. Some heteroskedasticity-consistent covariance matrix estimators with improved finite sample properties. Journal of Econometrics. 1985;29(3):305-325.
Long JS, Ervin LH. Using heteroscedasticity consistent standard errors in the linear regression model. The American Statistician. 2000;54(3):217-224.
Zeileis A. Object-oriented computation of sandwich estimators. Journal of Statistical Software. 2006;16(9):1-16.
Mansournia MA, Nazemipour M, Naimi AI, Collins GS, Campbell MJ. Reflection on modern methods: demystifying robust standard errors for epidemiologists. International Journal of Epidemiology. 2021;50(1):346-351.
