再入院の有無を来院ごとに記録した縦断研究では、同じ患者から得た4回分の観測が互いに似た値をとるため、観測を独立と見なした解析は標準誤差を小さく見積もります。一般化推定方程式（GEE）は、この相関の形をあらかじめ完全に言い当てなくても、集団全体での平均的な効果を推定できる枠組みで、Liang and Zeger（1986）が提案しました。
ところが、同じデータに変量切片つきの混合モデルを当てはめて得た治療効果と、GEEで得た治療効果は一致しません。乱数で生成した400例1,600観測のデータでこの差を数値にし、どちらを主解析として報告するのかを研究目的から決める順序まで扱います。
400例を4回追う再入院データの設定
ここで使うのは、治療群200例と対照群200例が4回ずつ来院し、各回で再入院の有無を記録した、ランダムに生成したデータセットです。個人差は標準偏差1.5の変量切片として置き、個人内の治療効果を対数オッズ比0.80、来院1回あたりの変化を-0.15としました。
この設定では潜在尺度で測った級内相関が0.406になり、同じ患者の4観測はかなり強く似た値をとります。記号は、$y_{ij}$を被験者$i$の$j$回目の再入院の有無（1が再入院）、$\mathbf{x}_{ij}$を説明変数のベクトル、$\boldsymbol{\beta}$を回帰係数、$\mu_{ij}=\mathbb{E}[y_{ij}]$を周辺平均とします。
なお、同じ縦断データでも、繰り返し測定した値と生存時間を同時に扱う場面では別の枠組みが要ります。その整理は縦断データと生存データの結合モデルで扱っているので、ここでは二値アウトカムの繰り返し測定に絞り、GEEが何を推定していて、混合モデルの係数とどこで食い違うのかを読み取れる状態を目指します。
一般化推定方程式が解いている式
GEEは尤度を立てません。周辺平均のモデルと作業相関だけを置いて、次の推定方程式を解きます。
$$\sum_{i=1}^{N}\mathbf{D}_i^{\top}\mathbf{V}_i^{-1}\left(\mathbf{y}_i-\boldsymbol{\mu}_i\right)=\mathbf{0}\tag{1}$$
式(1)は、各被験者の観測値と周辺平均の差を作業共分散行列$\mathbf{V}_i$で重み付けて足し合わせ、その合計が0になる$\boldsymbol{\beta}$を探すことを表しています。$\mathbf{D}_i$は$\boldsymbol{\mu}_i$を$\boldsymbol{\beta}^{\top}$で微分した行列、$\mathbf{V}_i$は作業相関$\mathbf{R}_i(\alpha)$と周辺分散から組み立てた行列です。
ここで効いてくるのは、$\mathbf{V}_i$が真の相関構造と一致している必要がない点です。一致していなくても$\mathbb{E}[\mathbf{y}_i-\boldsymbol{\mu}_i]=\mathbf{0}$が成り立つ限り推定方程式の期待値は0のままなので、$\hat{\boldsymbol{\beta}}$は周辺平均モデルのパラメータへ一致します。
もっとも、作業相関の誤りを引き受けるのは標準誤差のほうであり、残差から分散を組み直すサンドイッチ推定量が使われます。
$$\widehat{\mathrm{Var}}(\hat{\boldsymbol{\beta}})=\mathbf{B}^{-1}\left(\sum_{i=1}^{N}\mathbf{D}_i^{\top}\mathbf{V}_i^{-1}\hat{\mathbf{r}}_i\hat{\mathbf{r}}_i^{\top}\mathbf{V}_i^{-1}\mathbf{D}_i\right)\mathbf{B}^{-1}\tag{2}$$
式(2)は、モデルから計算した$\mathbf{B}=\sum_{i=1}^{N}\mathbf{D}_i^{\top}\mathbf{V}_i^{-1}\mathbf{D}_i$で挟む形にして、実際に観測された残差の広がりを分散へ反映させることを表しています。$\hat{\mathbf{r}}_i=\mathbf{y}_i-\hat{\boldsymbol{\mu}}_i$は被験者$i$の残差ベクトルです。
作業相関を外しても点推定は動かない
今回のデータで交換可能構造を指定すると、推定された共通相関は0.244、尺度パラメータは1.002で、反復5回で収束しました。治療効果の推定値は対数オッズ比0.5711、ロバスト標準誤差は0.1368でした。
実際、同じデータで作業相関を独立に置き替えても、点推定は0.5713とほとんど動かず、通常のロジスティック回帰の係数とは小数点以下6桁まで一致します。動いたのは標準誤差のほうで、相関を無視したモデル依存の標準誤差は0.1038と、ロバスト標準誤差0.1368の0.759倍まで縮みました。
0.759という比の意味は、独立性を仮定したまま反復測定を解析したときに何が起きるか、という話そのものです。信頼区間は実際より狭くなり、有意になりやすい向きへ偏ります。
一方、交換可能構造を指定した場合は、ロバスト標準誤差とモデル依存の標準誤差の比が1.002で、ほとんど差がありません。作業相関が実際の相関に近いほど2つの標準誤差は近づくので、この乖離の大きさ自体が作業相関の当てはまりを見る材料になります。
集団平均の効果と個人内の効果は別の量
GEEと混合モデルの結果が食い違ったとき、どちらかが偏っているのだろうと読まれがちです。ただ、二値アウトカムの場合、両者はそもそも別のパラメータを推定しており、値の違いは偏りではなく推定対象の違いから生じます（Neuhaus et al. 1991）。
ここで、ある患者本人が治療を受けた場合と受けなかった場合を比べた対数オッズ比、つまり個人内の効果として0.80を置きました。これに対し、治療群全体の再入院割合と対照群全体の再入院割合から計算した対数オッズ比は、変量切片について数値積分すると0.5707になります。
同じデータ生成の仕組みから、2つの違う値が出てきます。
この縮みはロジットリンクの非線形性から来ており、個人ごとのロジスティック曲線を個人差の分布で平均すると、平均された曲線は元の曲線より傾きが緩やかになります。傾きが緩むぶん、対応する対数オッズ比は0へ寄ります（図1）。
Zeger et al.（1988）は、この関係を近似式の形で与えています。
$$\beta_{\mathrm{PA}}\simeq\frac{\beta_{\mathrm{SS}}}{\sqrt{1+0.346\,\sigma^{2}}}\tag{3}$$
式(3)は、変量切片の標準偏差$\sigma$が大きいほど、集団平均の対数オッズ比$\beta_{\mathrm{PA}}$が個人内の対数オッズ比$\beta_{\mathrm{SS}}$より0へ縮むことを表しています。$\sigma=1.5$を入れると縮小係数は1.3336になり、0.80を割った値は0.5999です。
この近似を数値積分で求めた0.5707と比べると、近似式の値は0.029だけ大きくなりました。$\sigma$が1を超える設定ではこの程度の誤差が出るため、論文へ書く数値は近似式からではなく、推定値か積分の結果から取ります。
どちらを報告するかは研究目的で決まる
集団平均の効果が答えているのは、対象集団全体へ治療を広げたときに再入院割合がどれだけ動くか、という研究目的です。個人内の効果が答えているのは、同じ患者が治療を受けた場合と受けなかった場合でオッズがどれだけ変わるか、という別の研究目的になります。
そのため、医療政策の評価や集団への介入効果を扱う研究では、GEEを勧めます。Hubbard et al.（2010）は、近隣環境の要因と健康アウトカムの関連を例に、集団平均モデルと混合モデルの選択が推定対象と仮定の置き方で決まることを整理しています。
逆に、個人の経過を予測したい場面や、個人差の大きさそのものが関心の対象になる場面では、変量効果を明示的に置く混合モデルのほうが向きます。変量切片の標準偏差は混合モデルでは推定されるパラメータですが、GEEの出力に対応する量は現れません。
なお、今回のように個人差が大きい設定では、推定対象の違いだけで0.80と0.571という開きが出るので、値が近いかどうかは選択の根拠にならないため、両方を当てはめて値の近いほうを採る、という進め方は適切ではありません。
GEEを勧めない条件
クラスター数が少ない研究になると、話が変わります。サンドイッチ推定量は残差の広がりからクラスター間のばらつきを推定するため、クラスター数が小さいと分散を小さめに見積もり、検定のサイズが名目の5%より膨らみます。
この偏りに対して、Mancl and DeRouen（2001）は、残差を当てはめ値の影響で補正した共分散推定量を提案しました。施設や病棟をクラスターとする研究のように、クラスター数が数十を下回る解析では、補正なしのロバスト標準誤差をそのまま報告するのは避けたほうがよいでしょう。
もう一点、欠測の扱いにも条件が付きます。GEEの一致性は欠測が完全にランダム（MCAR）かどうかに依存するので、脱落が予後と関連する研究では、逆確率重み付けを加えるか、尤度に基づく混合モデルへ移るかの判断が要ります。
なお、今回のデータは4回分が全て揃った設計なので、この問題は現れていません。実際の縦断研究では、来院の抜けがどの程度あるかを先に数えておくほうが、作業相関の構造をどれにするかより結果に効きます。
論文でGEEの結果をどう書くか
報告で最低限そろえるのは、周辺平均モデルの形（リンク関数と説明変数）、指定した作業相関の構造、標準誤差がロバストかどうか、そしてクラスターの定義と数です。作業相関の構造だけが書かれていて標準誤差の種類が分からない報告では、読者が推定量の性質を判断できません。
また、係数の解釈も集団平均の量だと分かる書き方にします。今回の結果なら、治療群のオッズが対照群の1.77倍（95%信頼区間1.35〜2.31）という言い方は集団平均の効果の説明ですが、ある患者が治療を受けるとオッズが1.77倍になる、と書けば個人内の効果の説明へすり替わります。
私なら、混合モデルとGEEの両方の結果が手元にある段階で、まず変量切片の標準偏差の推定値を見ます。この値が小さければ2つの係数はほぼ重なるので解釈は割れませんが、今回のように1を超えていれば、どちらを主解析にするかを解析計画の段階へ戻って確認します。
Rで確かめる
ここからは、再入院の有無を4回分並べたデータを乱数で生成し、geepackパッケージのgeeglm関数で交換可能構造のGEEを当てはめます。個人差の標準偏差1.5、個人内の対数オッズ比0.80という設定は、先に書いたとおりです。
set.seed(42)
n_subject <- 400
subject_id <- rep(seq_len(n_subject), each = 4)
visit_index <- rep(0:3, times = n_subject)
treatment <- rep(rep(c(0, 1), each = n_subject / 2), each = 4)
random_intercept <- rnorm(n_subject, mean = 0, sd = 1.5)
linear_predictor <- -0.80 + 0.80 * treatment - 0.15 * visit_index +
random_intercept[subject_id]
readmission <- rbinom(length(linear_predictor), size = 1,
prob = plogis(linear_predictor))そのうえで、被験者を単位にしたクラスターを指定して周辺平均モデルを当てはめ、治療効果の推定値とロバスト標準誤差を取り出します。
analysis_data <- data.frame(readmission = readmission, treatment = treatment,
visit = visit_index, subject_id = subject_id)
gee_fit <- geepack::geeglm(readmission ~ treatment + visit, id = subject_id,
data = analysis_data, family = binomial(),
corstr = "exchangeable")
round(summary(gee_fit)$coefficients["treatment", c("Estimate", "Std.err")], 4)#> Estimate Std.err
#> 0.5711 0.1368治療効果の推定値は対数オッズ比0.5711、ロバスト標準誤差は0.1368で、95%信頼区間は0.303〜0.839でした。オッズ比に直すと1.77（1.35〜2.31）、z値は4.17です。
ここで読み違えやすいのは、0.5711がデータ生成で置いた0.80を外しているように見える点です。GEEが推定しているのは集団平均の効果なので、比べる相手は0.80ではなく、変量切片について積分した周辺の対数オッズ比のほうになります。
marginal_probability <- function(eta, sigma) {
integrate(function(b) plogis(eta + b) * dnorm(b, mean = 0, sd = sigma),
lower = -8 * sigma, upper = 8 * sigma)$value
}
marginal_logit <- function(trt) qlogis(vapply(0:3, function(visit) {
marginal_probability(-0.80 + 0.80 * trt - 0.15 * visit, 1.5)
}, numeric(1)))
round(mean(marginal_logit(1) - marginal_logit(0)), 4)#> [1] 0.5707その相手にあたる値は0.5707で、推定値0.5711との差はロバスト標準誤差の0.03倍ほどに収まりました。観測された周辺割合から直接計算した対数オッズ比0.5719も同じ水準です。
ここで曲線を描くと、ここまでの数値の背景にある形が読み取れます。個人ごとのロジスティック曲線と、それらを変量切片の分布で平均した曲線を、同じ軸へ並べます。
eta_grid <- seq(-6, 6, length.out = 401)
marginal_curve <- vapply(eta_grid, marginal_probability, numeric(1), sigma = 1.5)
plot(eta_grid, plogis(eta_grid), type = "n", ylim = c(0, 1),
xlab = "Fixed-effect linear predictor", ylab = "Probability of readmission")
for (b in seq(-3, 3, by = 0.75)) {
lines(eta_grid, plogis(eta_grid + b), col = "grey78", lwd = 1)
}
lines(eta_grid, marginal_curve, col = "#1f4e79", lwd = 3)
もう一枚は、当てはめた周辺平均モデルが来院ごとの再入院割合をどの程度再現しているかです。観測された割合の点と、あてはめ値から作った曲線を重ねます。
observed_probability <- tapply(readmission, list(treatment, visit_index), mean)
gee_fitted <- function(trt) plogis(coef(gee_fit)[1] + coef(gee_fit)[2] * trt +
coef(gee_fit)[3] * 0:3)
plot(0:3, observed_probability["1", ], ylim = c(0.2, 0.62), pch = 16,
col = "#b23a48", xlab = "Visit", ylab = "Proportion readmitted")
points(0:3, observed_probability["0", ], pch = 17, col = "#1f4e79")
lines(0:3, gee_fitted(1), col = "#b23a48", lwd = 2)
lines(0:3, gee_fitted(0), col = "#1f4e79", lwd = 2)
図2で確認できるのは、GEEが個々の患者の軌跡ではなく、各時点での群の割合を再現しているという点です。個人の変動を説明したい研究でこの当てはまりの良さを根拠にするのは、見ている対象がズレています。
選ぶ前に決めておくこと
GEEと混合モデルのどちらを使うかは、当てはまりの良さや有意性ではなく、報告したい効果が集団の割合の変化なのか、個々の患者のオッズの変化なのかで決まります。二値アウトカムで個人差が大きい研究では、この違いが今回の0.5711と0.80のような開きになって現れました。
そのうえで、GEEを選んだ場合はクラスター数と欠測の状況を確認します。クラスター数が数十を下回るならMancl and DeRouen（2001）の補正を考え、脱落が予後と関連するなら重み付けか尤度に基づく方法へ移る判断が必要になります。
ここまで読んで、自分の研究ではどちらを主解析にすべきか迷った方もいるかもしれません。実際のデータでは、変量効果の分散、クラスター数、欠測の起こり方、そして事前に定めた解析計画まで含めて判断することになります。
そうした判断の続きは、Dr.データサイエンスでご相談を承っています。
参考文献
- Liang KY, Zeger SL. Longitudinal data analysis using generalized linear models. Biometrika. 1986;73(1):13-22.
- Zeger SL, Liang KY, Albert PS. Models for longitudinal data: a generalized estimating equation approach. Biometrics. 1988;44(4):1049-1060.
- Neuhaus JM, Kalbfleisch JD, Hauck WW. A comparison of cluster-specific and population-averaged approaches for analyzing correlated binary data. International Statistical Review. 1991;59(1):25-35.
- Mancl LA, DeRouen TA. A covariance estimator for GEE with improved small-sample properties. Biometrics. 2001;57(1):126-134.
- Hubbard AE, Ahern J, Fleischer NL, Van der Laan M, Lippman SA, Jewell N, Bruckner T, Satariano WA. To GEE or not to GEE: comparing population average and mixed models for estimating the associations between neighborhood risk factors and health. Epidemiology. 2010;21(4):467-474.
