HC3を指定しても直らないものがあります。平均モデルの形が間違っていれば、分散をどう計算しても推定している量がズレたままですし、少数の外れ値が係数を引っ張っている状態も変わりません。分散の補正が扱うのは、係数の周りの不確実性だけです。
そして、10施設から40人ずつ集めたデータでは、不等分散とは別の依存が入ります。同じ施設の患者は診療手順や患者背景を共有していて、400人分の観測が400個の独立した情報にはなりません。HC3は対象者ごとに残差の大きさを使いますが、対象者同士が関連していることは扱いません。
この記事では、分散補正が直さない問題を先に確認してから、施設内相関を持つ400人のデータで、クラスター頑健な標準誤差がどこまで区間を広げるか、自由度をどう数えるか、施設数が少ない場合のCR2までを扱います。不等分散の補正は、HC0からHC3までで扱いました。
平均モデルのズレと推定対象
severityと測定値の平均の関係が強く曲がっているのに直線だけを当てはめた場合、標準誤差をHC3へ変えても、当てはめた直線は曲がりません。また、必要な交絡因子をモデルから外していたり、説明変数に測定誤差があったりすれば、分散の補正とは別の問題が残ります。
ここで、平均モデルが完全に正しくなければロバスト標準誤差に意味がない、とまで言い切るのも正確ではありません。対象者を母集団からの独立なランダム標本と考える設定では、OLSが推定する母集団の最良線形近似を、推定対象として定義できます。
$$\boldsymbol{\beta}^{*} =\arg\min_{\mathbf{b}} \mathbb{E}\!\left[(Y-\mathbf{x}^{\top}\mathbf{b})^{2}\right] \tag{1}$$
式(1)は、真の条件付き平均が直線でなくても、母集団で予測の二乗誤差を最も小さくする直線の係数を推定対象として定めています。
このランダムな説明変数を含む標本抽出の設定で、モーメントや行列の非退化性などの条件が満たされれば、サンドイッチ分散は$\boldsymbol{\beta}^{*}$を推定する不確実性を扱えます。ただし、その係数が、研究者の知りたい一定の治療効果や、全ての患者に共通する傾きになるとは限りません。説明変数を固定した条件付き推論と、ランダム標本についての推論を混ぜないことも必要です。
推定対象を変えずに分散だけ直したつもりが、もともとの係数の意味を曖昧なままにしている場合があります。KingとRobertsは、通常の標準誤差とロバスト標準誤差の差を、モデルを見直す手掛かりとして読むことを論じています。ただし、差が小さいことだけで平均モデルの正しさが証明されるわけではありません。
正規性と外れ値への対応
不等分散の例では、説明変数で条件付けた誤差は正規分布から生成しました。それでも通常の標準誤差とHC3は異なっています。正規性と等分散性は別の仮定で、正規性の検定を通ったことだけでは、通常の回帰出力を選ぶ根拠になりません。
逆に、アウトカムの分布が正規分布でなくても、平均を推定する目的に線形回帰が合い、独立性やモーメントなどの条件が満たされるなら、大標本近似による推論が候補になります。Lumleyらは、公衆衛生データで正規性を機械的な手法選択基準にすることの問題を論じています。
ただし、極端に重い裾を持つ分布や、少数の観測値が推定量の変動をほぼ決めるデータでは、通常の中心極限定理に基づく近似が不安定になります。ロバストという名前から、極端な観測値への耐性まで備わると考えるのは注意が必要です。
外れ値の影響を減らすロバスト回帰は、係数を求める目的関数や重みを変更する方法です。これに対して、OLSにHC3を組み合わせるだけなら、外れ値に引かれた係数はそのまま残ります。データ入力の誤り、対象集団の違い、平均の非線形性を確認する作業を、標準誤差の変更で置き換えることは推奨しません。
施設内で相関する400人のデータ
同じ施設の患者が共通した測定環境や診療体制の影響を受ける場合、患者同士の誤差が独立とは限りません。施設差をモデルの側で表す方法は、frailtyモデルと多施設データで扱っています。これは残差の分散が患者ごとに異なる問題とは別で、個人単位のHC3を指定しただけでは扱えません。
次の例では40施設から各10人、合計400人のデータを生成します。施設単位で治療群と対照群へ20施設ずつ割り付け、同じ施設の患者には共通する成分を加えます。施設内相関の意味を数式で確認するための生成モデルで、治療方針に関する臨床的な提案ではありません。
$$Y_{gi}=100+2A_g+b_g+e_{gi}, \qquad b_g\sim N(0,9), \qquad e_{gi}\sim N(0,4) \tag{2}$$
式(2)は、治療群の平均を対照群より2大きくし、施設$g$に共通する変動$b_g$と、患者$i$に固有の変動$e_{gi}$を加えています。
$A_g$は施設の割付で、対照群が0、治療群が1です。異なる施設の$b_g$は独立とし、患者固有の誤差も互いに独立で、割付と各誤差成分も独立に生成します。同じ施設では$b_g$を共有するため、患者の測定値が一緒に高くなったり低くなったりします。
cluster_id <- rep(seq_len(40), each = 10)
treatment <- rep(sample(rep(c(0, 1), each = 20)), each = 10)
cluster_measurement <- 100 + 2 * treatment + rep(rnorm(40, sd = 3), each = 10) + rnorm(400, sd = 2)
cluster_fit <- lm(cluster_measurement ~ treatment)
cluster_se <- sqrt(c(OLS = vcov(cluster_fit)[2, 2], HC3 = vcovHC(cluster_fit, type = "HC3")[2, 2], Cluster_HC1 = vcovCL(cluster_fit, cluster = cluster_id, type = "HC1", cadjust = TRUE)[2, 2]))
cluster_df <- c(398, 398, 39)
cluster_results <- cbind(Estimate = coef(cluster_fit)[2], SE = cluster_se, df = cluster_df, Lower = coef(cluster_fit)[2] - qt(0.975, df = cluster_df) * cluster_se, Upper = coef(cluster_fit)[2] + qt(0.975, df = cluster_df) * cluster_se)
print(round(cluster_results, digits = 3))#> Estimate SE df Lower Upper
#> OLS 2.792 0.380 398 2.044 3.540
#> HC3 2.792 0.381 398 2.043 3.542
#> Cluster_HC1 2.792 1.042 39 0.685 4.899治療群の係数は2.792で、通常の標準誤差は0.380、個人単位のHC3は0.381、施設単位のvcovCL()では1.042となりました。この結果では、HC3へ変更しても、施設内相関を無視したことによる差がほとんど残っています。
出力のCluster_HC1は、施設ごとにまとめた寄与に、後述する有限標本補正を加えたものです。個人単位のHC1と同じ呼び方の部分がありますが、中央の分散計算でまとめる単位が異なります。また、このコードの信頼区間は、通常のOLSと個人単位HC3には自由度398、施設単位の補正には自由度39のt分布を使っており、分散の補正だけでなく参照分布も変更した計算です。
png("robust-cluster-residuals.png", width = 1100, height = 800, res = 200, type = "cairo")
plot(cluster_id, y = residuals(cluster_fit), pch = 16, col = "#0072B280", xlab = "Cluster ID", ylab = "OLS residual", main = "Shared variation within each cluster")
abline(h = 0, lty = 2)
dev.off()
図1は、残差の広がりだけでなく、どの観測値同士が一緒に動くのかを見る図です。ただし、実際の研究でこの形が目立たない場合も、施設単位の割付や標本抽出があるなら、その構造を図の見た目だけで無視することはできません。frailtyモデルと多施設データでも扱うように、施設内の依存をどう扱うかは、解析前に確認したい事項です。
施設ごとに寄与を足す意味
独立した患者のHC0では、各患者の$\mathbf{x}_i\widehat{\varepsilon}_i$を使いました。クラスター・ロバスト分散では、この寄与を施設内で先に足し、その施設全体の寄与がどれくらい変動するかを使います。
施設数を$G$、施設$g$の説明変数行列と残差ベクトルを$\mathbf{X}_g$、$\widehat{\boldsymbol{\varepsilon}}_g$とすると、補正前の分散は次の形です。
$$\widehat{\mathbf{V}}_{\mathrm{CL0}} =(\mathbf{X}^{\top}\mathbf{X})^{-1} \left\{\sum_{g=1}^{G}\widehat{\mathbf{s}}_g \widehat{\mathbf{s}}_g^{\top}\right\} (\mathbf{X}^{\top}\mathbf{X})^{-1}, \qquad \widehat{\mathbf{s}}_g =\mathbf{X}_g^{\top}\widehat{\boldsymbol{\varepsilon}}_g \tag{3}$$
式(3)は、患者ごとに独立した寄与を足す代わりに、施設ごとにまとめた寄与の変動から、係数の分散を推定しています。
施設内で足してから二乗するので、その施設の異なる患者の寄与同士を掛けた項も残ります。この項が、個人単位のHC0には含まれていなかった部分です。施設内の相関構造を特定の形へ固定せずに扱えることが利点ですが、通常の大標本理論は、独立した施設が十分に増えることや、一部の施設が情報の大半を占めないことなどに依存します。
「施設内の相関は自由に扱える」と「施設間の相関も無視してよい」は同じ意味ではありません。複数施設が同じ地域の影響を強く共有する場合には、施設単位でまとめても、施設同士の独立性が残らない可能性があります。CameronとMillerは、何をクラスターとするかを、研究の構造から判断する必要性を詳しく論じています。
分散補正と信頼区間の自由度
今回使ったvcovCL()のtype="HC1"、cadjust=TRUEでは、式(3)に次の係数を掛けます。ここで$p$は切片を含む係数の数です。
$$\widehat{\mathbf{V}}_{\mathrm{CL1}} =\frac{G}{G-1}\frac{n-1}{n-p} \widehat{\mathbf{V}}_{\mathrm{CL0}} \tag{4}$$
式(4)は、施設数と観測数に応じた有限標本補正を、補正前のクラスター分散へ加えています。
この定義は今回使用したsandwichパッケージの仕様です。CR1という名称の補正係数はソフトウェアによって異なる場合があるため、CR1とだけ記載するより、関数、型、補正設定を残した方が再現しやすくなります。
係数$\widehat{\beta}_j$について、自由度$\nu$のt分布を使う両側95%信頼区間は、次の形です。
$$\widehat{\beta}_j \pm t_{0.975,\nu}\, \sqrt{\widehat{\mathbf{V}}_{jj}} \tag{5}$$
式(5)は、係数の推定値を中心として、標準誤差にt分布の97.5%点を掛けた幅を左右へ取っています。
分散推定量を変更する処理と、式(5)の$\nu$を決める処理は別です。ロバスト分散行列を計算した後に、通常の回帰出力のp値や信頼区間をそのまま使えば、補正前の推論が残ることがあります。また、vcovCL()は分散共分散行列を返す関数で、自由度や信頼区間を自動的に決める関数ではありません。
今回の施設単位の区間には、よく使われる近似として$\nu=G-1=39$を指定しました。ただし、これによって厳密なt分布が成立するわけではなく、施設数が少ない場合や施設間の情報量が偏っている場合には、別の小標本対応が必要になります。
png("robust-cluster-intervals.png", width = 1100, height = 800, res = 200, type = "cairo")
plot(cluster_results[, "Estimate"], y = 3:1, xlim = range(cluster_results[, c("Lower", "Upper")], 2), ylim = c(0.5, 3.5), yaxt = "n", pch = 16, xlab = "Treatment coefficient and 95% interval", ylab = "", main = "Same estimate, different uncertainty")
segments(cluster_results[, "Lower"], y0 = 3:1, x1 = cluster_results[, "Upper"], y1 = 3:1, lwd = 3, col = c("#777777", "#009E73", "#0072B2"))
axis(2, at = 3:1, labels = c("OLS", "HC3", "Cluster"), las = 1, cex.axis = 0.8)
abline(v = 2, lty = 2, col = "#D55E00")
dev.off()
この標本では、通常のOLSと個人単位HC3の区間は真の効果2を含まず、施設単位の区間は含んでいます。ただし、1回の標本で真の値を含んだことを、被覆率が95%に達した証拠として使うことはできません。全ての区間で下限は0より大きいため、この例の主題は有意・非有意の切り替えではなく、独立した情報をどの単位で数えるかによる不確実性の差です。
施設内相関と真の標準誤差
式(2)では、施設共通成分の分散を9、患者固有成分の分散を4と置きました。この生成モデルなら、同じ施設の患者同士の級内相関係数、ICCを計算できます。また、全ての施設の人数を$m$とした均等な設計では、施設内相関による平均の分散の増加を表すデザイン効果も求められます。
$$\rho=\frac{\tau^2}{\tau^2+\sigma^2}, \qquad \mathrm{DE}=1+(m-1)\rho \tag{6}$$
式(6)は、施設共通成分の分散$\tau^2$が全体の分散に占める割合をICCとし、1施設に$m$人いる場合に、独立な観測と比べて平均の分散が何倍になるかを表しています。
ここでの$\sigma^2$は患者固有成分の分散です。今回の均等な施設単位割付では、同じデザイン効果を、群間平均差の分散の比較にも使えます。ただし、一般の調整回帰、施設人数が大きく異なるデータ、施設内で治療が変わるデータへ、この単純な倍率をそのまま移すことはできません。
治療群と対照群の施設数をそれぞれ$G_1$、$G_0$とし、施設ごとの人数が全て$m$なら、今回のモデルでの群間平均差の分散は次の形になります。
$$\operatorname{Var}(\widehat{\beta}_{A}) =\left(\tau^2+\frac{\sigma^2}{m}\right) \left(\frac{1}{G_1}+\frac{1}{G_0}\right) \tag{7}$$
式(7)は、患者を施設内で平均した際に残る分散を、治療群と対照群の施設数に応じて足し合わせています。
患者固有成分の分散は施設内の人数を増やすと小さくなりますが、施設共通成分はそのまま残ります。患者数が増えても、施設数を増やした場合と同じ精度が得られるとは限らない理由が、この式から分かります。
print(round(c(ICC = 9 / 13, Design_effect = 1 + 9 * 9 / 13, Cluster_oracle_SE = sqrt((9 + 4 / 10) * (1 / 20 + 1 / 20))), digits = 3))#> ICC Design_effect Cluster_oracle_SE
#> 0.692 7.231 0.970ICCは0.692、デザイン効果は7.231、生成モデルから求めた係数の真の標準誤差は0.970です。施設単位の補正で推定した1.042は、この真の値そのものではなく、40施設の標本から推定した値なので、標本によって変動します。
この数値は、実際の研究でICCが0.692になることを想定すべきだという主張ではありません。相関の影響を見やすくするために置いた生成条件です。また、患者数をデザイン効果で割った値を求めても、それを全ての回帰係数に共通する自由度として使えるわけではありません。
少数クラスターとCR2
40施設の例で計算が動いたことと、40施設あればどの研究でも十分ということは別です。同じ施設数でも、人数が偏っていたり、治療群に属する施設が少なかったり、特定の係数がほぼ1施設の情報で決まったりすれば、推論は不安定になります。
患者数が多くても、治療群の割付が施設単位なら、その治療効果について独立した割付情報の数は施設数のままです。したがって、30施設や50施設を一律の合格基準にするより、群別の施設数、施設人数の分布、施設ごとのレバレッジ、推定したい係数への寄与を確認します。MacKinnonらは、クラスター数だけでは捉えきれない情報の偏りを含めて、実際の推論方法を論じています。
小標本への対応として、CR2とSatterthwaite近似を組み合わせる方法があり、CR2は、作業用の分散モデルの下で分散推定の偏りを減らすように、残差の寄与を補正します。通常のOLSで独立等分散を作業モデルとし、必要な逆行列が存在する場合、その形は次のように書けます。
$$\begin{aligned} \widehat{\mathbf{V}}_{\mathrm{CR2}} &=(\mathbf{X}^{\top}\mathbf{X})^{-1} \left\{\sum_{g=1}^{G} \mathbf{X}_g^{\top}\mathbf{A}_g \widehat{\boldsymbol{\varepsilon}}_g \widehat{\boldsymbol{\varepsilon}}_g^{\top} \mathbf{A}_g^{\top}\mathbf{X}_g \right\} (\mathbf{X}^{\top}\mathbf{X})^{-1},\\ \mathbf{A}_g&=(\mathbf{I}_{m_g}-\mathbf{H}_{gg})^{-1/2}. \end{aligned} \tag{8}$$
式(8)は、施設内の残差ベクトルに行列による補正を加え、その施設が回帰の当てはめへ与えた影響を反映して、分散推定の偏りを減らしています。
$m_g$は施設$g$の人数、$\mathbf{I}_{m_g}$はその大きさの単位行列、$\mathbf{H}_{gg}$はハット行列の施設内部分です。施設固定効果を含む場合など、単純な逆行列が存在しない構造では、この式をそのまま計算するのではなく、対応する一般化が必要です。
また、CR2の標準誤差を計算しただけで小標本推論が完成するわけではありません。Satterthwaite近似では、係数ごとの情報の分布から自由度を近似するため、同じ回帰モデルでも係数によって自由度が異なり、$G-1$よりかなり小さくなる場合があります。PustejovskyとTiptonの方法と、clubSandwichパッケージのvcovCR()、coef_test()は、この部分を区別して扱います。
今回掲載した計算は、sandwichによる施設単位HC1と自由度39の近似で、CR2を計算した結果ではありません。少数施設の研究へ移す際には、CR2と小標本自由度、研究デザインに合うワイルド・クラスター・ブートストラップなどを比較候補として解析計画に組み込みます。ただし、治療群がごく少数の施設にしかない場合など、独立した情報が乏しい問題を、補正方法だけで解消できるとは限りません。
広がった区間が表しているもの
治療群の係数2.792は、OLSでもHC3でもクラスター頑健でも同じでした。動いたのは区間の幅で、0.380から1.042へ、およそ2.7倍になっています。
この差は、400人を400個の独立した情報として数えていたことから来ています。施設内相関0.692を考えると、実質的に使えている情報は400人分よりはるかに少なく、その事実を区間の幅が受け取っただけです。
次の記事では、クラスターの単位をどう決めるか、固定効果・混合モデル・GEEとの関係、そしてサンドイッチ分散の一般形から二値アウトカムの修正Poisson回帰までを扱います。
参考文献
Cameron AC, Miller DL. A practitioner’s guide to cluster-robust inference. Journal of Human Resources. 2015;50(2):317-372.
Pustejovsky JE, Tipton E. Small-sample methods for cluster-robust variance estimation and hypothesis testing in fixed effects models. Journal of Business & Economic Statistics. 2018;36(4):672-683.
MacKinnon JG, Nielsen MØ, Webb MD. Cluster-robust inference: a guide to empirical practice. Journal of Econometrics. 2023;232(2):272-299.
Zeileis A, Köll S, Graham N. Various versatile variances: an object-oriented implementation of clustered covariances in R. Journal of Statistical Software. 2020;95(1):1-36.
Lumley T, Diehr P, Emerson S, Chen L. The importance of the normality assumption in large public health data sets. Annual Review of Public Health. 2002;23(1):151-169.
