患者の収縮期血圧を解析するとき、同じ医師が診た患者は互いに似やすく、その医師も同じ施設の方針や対象集団を共有するため、患者1,800人がいても、医師単位の曝露を比較する情報が独立な1,800人分あるとは考えられません。患者が医師に、医師が施設に入る三階層を無視すると、標準誤差や施設間・医師間のばらつきの読み方が崩れるため、ここでは単純な三階層ランダム切片モデルに範囲を絞ります。
患者・医師・施設という三階層
多階層データとは、観測単位が上位の集団にまとめられ、同じ集団内の観測値に依存が生じうるデータです。Diez-Roux（2000）は、個人レベルと集団レベルの要因を同時に扱い、両階層のばらつきを分けて考える枠組みとして多階層解析を説明しています。
また、患者を第1階層、医師を第2階層、施設を第3階層とすると、同じ医師の患者は医師と施設の両方を共有し、同じ施設で別の医師に診られた患者は施設だけを共有します。この違いがあるため、全ての患者間に同じ相関を置く共分散構造だけでは、二種類の共有関係を表し切れません。
ここで最初に確かめたいのは、所属が本当に入れ子になっているかです。患者が複数の医師を受診する、医師が複数施設に勤務するといった交差分類があれば、この記事の単純な三階層モデルをそのまま当てはめる方法は適切ではありません。
ランダム切片で依存を分ける
施設$j$の医師$k$が診る患者$i$の収縮期血圧を$Y_{ijk}$とし、医師単位の曝露$E_{jk}$と患者年齢$A_{ijk}$を説明変数にします。三階層ランダム切片モデルは、次の形で書けます。
$$Y_{ijk}=\beta_0+\beta_1E_{jk}+\beta_2A_{ijk}+u_j+v_{jk}+\varepsilon_{ijk}.\tag{1}$$
式(1)は、収縮期血圧を固定効果と、施設ごとの偏差$u_j$、施設内の医師ごとの偏差$v_{jk}$、患者ごとの残差$\varepsilon_{ijk}$に分けています。$\beta_1$は年齢を一定としたときの曝露群と非曝露群の平均差であり、この観察的な関連を因果効果とは呼びません。
ランダム切片は、施設や医師によって基準となる血圧水準が動くモデルです。一方、曝露と血圧の関連そのものが施設ごとに異なる状況ならランダム傾きが候補になりますが、その詳細は今回の範囲に含めません。
ただし、ランダム切片を加えるだけで施設レベルの交絡が消えるわけではなく、施設の診療方針が医師の曝露と患者の血圧の両方に関係する場合、その方針を測定して固定効果に含めるか、施設内の比較と施設間の比較を分ける必要があり、ランダム効果は未測定交絡の代わりにはなりません。
また、患者年齢のような第1階層の変数では、施設平均との差と施設平均を分けることで、同じ施設内で年齢が異なる患者の比較と、平均年齢が異なる施設同士の比較を区別できます。中心化を行う場合は、どの階層の平均を引いたかによって係数の意味が変わるため、式と方法欄の両方に基準を記載します。
分散分割係数が表す共有範囲
正規応答のランダム切片モデルでは、施設分散を$\sigma_u^2$、医師分散を$\sigma_v^2$、患者残差分散を$\sigma_e^2$とします。Goldsteinら（2002）が扱った分散分割の考え方に沿うと、施設の分散分割係数と同じ医師の患者間相関は次のようになります。
$$\mathrm{VPC}_{\mathrm{facility}}=\frac{\sigma_u^2}{\sigma_u^2+\sigma_v^2+\sigma_e^2},\qquad
\mathrm{ICC}_{\mathrm{same\ physician}}=\frac{\sigma_u^2+\sigma_v^2}{\sigma_u^2+\sigma_v^2+\sigma_e^2}.\tag{2}$$
式(2)の左は全分散のうち施設差が占める割合で、右は同じ医師に属する患者が施設差と医師差の両方を共有するために生じる相関です。同じ施設でも医師が異なる患者は医師効果を共有しないので、その相関は施設の分散分割係数と等しくなります。
ランダムに生成したデータセットでは、施設、医師、患者残差の推定分散がそれぞれ17.14、9.26、66.03でした。したがって、施設の分散分割係数は0.185、同じ医師の患者間相関は0.286となり、後者を施設内の全患者へ一律に当てはめると共有範囲を読み違えます。
ただし、応答尺度とリンク関数が変われば残差分散の扱いも変わるため、この式を二値アウトカムや計数アウトカムへそのまま移すことはできません。あくまで式(2)は正規応答のモデルに限った説明であることにご留意ください。
医師単位の曝露は医師数にも依存する
例では30施設に各5人の医師、各医師に12人の患者を置き、合計150人の医師と1,800人の患者を生成しました。曝露は医師単位で割り当てられており、曝露医師90人に患者1,080人、非曝露医師60人に患者720人が属します。
この生成標本で三階層モデルを当てはめると、曝露群と非曝露群の推定平均差は-3.77mmHg、標準誤差は0.64mmHg、95%信頼区間は-5.04～-2.51mmHgでした。生成時には差を-4mmHgと置いており、今回の一標本では推定値が0の方向へ0.23mmHg動いています。
一方、患者1,800人を独立として扱った線形回帰では、推定平均差は-3.78mmHgと近いものの、標準誤差は0.46mmHgまで小さくなりました。標準誤差の比は0.72ですが、階層を無視すれば常にこの割合になるという結果ではなく、この生成条件と一標本での観察です。
また、推定値が近かったことを根拠に、独立な線形回帰で十分だったとは判断できません。医師単位の曝露に関する精度は患者数だけで決まらず、曝露が変化する医師数と、その医師が属する施設数にも支えられます。
標本数は各階層で数える
Moineddinら（2007）は、患者が医師に、医師が診療組織に属する三階層ロジスティック回帰のシミュレーションを行い、各階層の標本数が推定へ影響することを報告しました。ただし、その論文で得られた標本数の推奨を、連続アウトカムを扱う本例の一般的な閾値にはしません。
同じ1,800人でも、少数の医師へ患者が集中したデータと、150人の医師へ分散したデータでは、医師レベル曝露を比較する情報量が異なります。そのため、患者総数だけで解析可能性を判断せず、施設数、施設当たり医師数、医師当たり患者数、曝露の階層を並べて確認します。
さらに、施設が30ある本例でも、施設レベルの説明変数を多数入れれば上位階層の情報は急に薄くなります。患者数を増やすだけでは上位階層の変数の変動は増えないので、係数の数やランダム効果の複雑さは上位単位数に合わせる必要があります。
そして、施設や医師ごとの患者数が均等でない場合は、大規模な施設が全体の推定へ強く寄与しますが、施設レベルの効果を支える単位数そのものは増えません。特に、患者が数人しかいない医師が多いと医師分散と患者残差を分けにくくなるため、総患者数に加えて各階層の単位数と人数分布を示します。
固定効果と分散を一緒に報告する
曝露係数だけを報告すると、施設間・医師間の異質性がどの程度残ったのかが見えません。Merloら（2005）は集団特性の固定効果と集団間分散を併せて解釈する重要性を述べており、多階層モデルを通常の回帰係数だけに縮めると、この利点が失われます。
その報告内容について、Monsalvesら（2020）のLEVELは、階層構造の図、各階層の標本数、固定効果、ランダム効果、分散分割係数を明示するよう推奨しています。本例なら30施設、150医師、1,800患者という構造と、曝露が医師レベルであることを方法に記載し、平均差と95%信頼区間に加えて各分散成分を結果へ載せます。
さらに、推定法が最尤法か制限付き最尤法か、ランダム効果にどの分布を置いたか、信頼区間をどの方法で求めたかも結果の再現に必要です。収束警告や特異な当てはまりが生じた場合は、複雑なモデルの結果だけを残さず、どの構造を簡略化したかと、その変更で固定効果がどう動いたかを記録します。
なお、p値が小さくても-3.77mmHgに臨床的な意味があるかは、この解析だけでは決まりません。臨床的に意味のある差は研究領域の根拠と研究目的に照らして別に定め、統計的な不確実性は95%信頼区間で表します。
Rで三階層モデルを確かめる
Rのnlmeパッケージを使い、同じデータを再現できる乱数生成条件の下で、30施設、150医師、1,800患者のデータセットを生成しました。施設、医師、患者残差の標準偏差は順に4、3、8mmHgと置き、年齢を調整した三階層ランダム切片モデルを制限付き最尤法で当てはめています。
set.seed(42)
n_facility <- 30; n_physician <- 150; n_patient <- 1800
facility <- factor(rep(1:n_facility, each = 60))
physician <- factor(rep(1:n_physician, each = 12))
physician_exposure <- unlist(lapply(1:n_facility, function(index) sample(c(0, 0, 1, 1, 1))))
exposure <- rep(physician_exposure, each = 12)
age <- stats::rnorm(n_patient, mean = 65, sd = 8)
facility_effect <- stats::rnorm(n_facility, mean = 0, sd = 4)
physician_effect <- stats::rnorm(n_physician, mean = 0, sd = 3)
sbp <- 135 - 4 * exposure + 0.25 * (age - 65) +
facility_effect[facility] + physician_effect[physician] + stats::rnorm(n_patient, sd = 8)
analysis_data <- data.frame(sbp, exposure, age_centered = age - 65, facility, physician)
fit <- nlme::lme(sbp ~ exposure + age_centered,
random = ~ 1 | facility / physician, data = analysis_data, method = "REML")variance_table <- nlme::VarCorr(fit)
variance <- as.numeric(variance_table[c(2, 4, 5), "Variance"])
c(exposure = nlme::fixef(fit)["exposure"],
facility_vpc = variance[1] / sum(variance),
same_physician_icc = sum(variance[1:2]) / sum(variance))#> exposure facility_vpc same_physician_icc
#> -3.7735922 0.1854307 0.2856027出力のexposureは年齢を調整した曝露群と非曝露群の平均差で、負の値は生成時の設定どおり曝露群の収縮期血圧が低い方向を表します。ここで読み違えやすいのは、facility_vpcを同じ医師の患者間相関とみなすことであり、同じ医師の患者は施設分散と医師分散の両方を共有します。
さらに、推定されたランダム切片を階層別に並べる図は、同じモデルから次のコードで再現できます。点の横方向の広がりは各単位の推定された偏差であり、分散成分そのものを点の範囲だけから読み取る図ではありません。
facility_random <- sort(nlme::ranef(fit, level = 1)[, 1])
physician_random <- sort(nlme::ranef(fit, level = 2)[, 1])
png("draft/figures/random-effects.png", width = 1100, height = 688, res = 200, type = "cairo")
par(mfrow = c(1, 2), mar = c(4, 4, 2, 1), cex = 0.72)
plot(facility_random, seq_along(facility_random), pch = 16,
xlab = "Random intercept (mmHg)", ylab = "Ordered facility", main = "Facility level")
abline(v = 0, lty = 2, col = "grey50")
plot(physician_random, seq_along(physician_random), pch = 16,
xlab = "Random intercept (mmHg)", ylab = "Ordered physician", main = "Physician within facility")
abline(v = 0, lty = 2, col = "grey50"); dev.off()
階層と曝露の単位を揃える
三階層データの解析では、患者数の多さより先に、誰が誰に属し、曝露がどの階層で変化するかを確定する必要があります。所属が入れ子なら施設と医師のランダム切片で共有関係を分け、固定効果の95%信頼区間と各階層の分散を一緒に読む方法が研究目的に合います。
逆に、交差分類を単純な入れ子へ押し込むことや、医師レベル曝露を患者レベルの独立な観測として扱うことは避けた方がよいでしょう。モデルを複雑にすること自体ではなく、データが生まれた階層と、知りたい関連の単位を一致させることが判断の中心です。
参考文献
- Diez-Roux AV. Multilevel analysis in public health research. Annual Review of Public Health. 2000;21:171–192.
- Goldstein H, Browne WJ, Rasbash J. Partitioning variation in multilevel models. Understanding Statistics. 2002;1(4):223–231.
- Merlo J, Chaix B, Yang M, Lynch J, Råstam L. A brief conceptual tutorial on multilevel analysis in social epidemiology: interpreting neighbourhood differences and the effect of neighbourhood characteristics on individual health. Journal of Epidemiology and Community Health. 2005;59(12):1022–1028.
- Moineddin R, Matheson FI, Glazier RH. A simulation study of sample size for multilevel logistic regression models. BMC Medical Research Methodology. 2007;7:34.
- Monsalves MJ, Bangdiwala AS, Thabane A, Bangdiwala SI. LEVEL (Logical Explanations & Visualizations of Estimates in Linear mixed models): recommendations for reporting multilevel data and analyses. BMC Medical Research Methodology. 2020;20(1):3.
