同じ人を繰り返し測ると測定値の間に相関が生じるため、一般化推定方程式（GEE）は集団平均の関係を表すモデルと、同じ人の測定値がどの程度似るかを仮に表す「作業相関構造」を組み合わせます。相関構造を実際のデータと違う形に指定したとき、回帰係数と標準誤差のどちらに何が起こるかは分けて考える必要があります。
ここでは240人を5時点で測ったランダムに生成したデータセットを使い、独立、交換可能、一次自己回帰という三つの相関指定を比較します。通常のモデルに基づく標準誤差とロバスト標準誤差を並べ、後者が補う誤指定の範囲と、平均モデル、欠測、独立した人数に関する条件を確認します。
作業相関は何のために置くか
GEEの基本を扱った記事では、個人ごとの効果ではなく、集団全体の平均に関する係数を推定しました。今回も対象者$i$の5時点の応答を一組として扱い、その平均$\mu_i$に対する回帰係数$\beta$を推定しますが、同一人の測定値が似る度合いは平均とは別に指定します。
$$V_i=A_i^{1/2}R_i(\alpha)A_i^{1/2},\qquad A_i=\operatorname{diag}\{\operatorname{Var}(Y_{it})\}\tag{1}$$
式(1)の$V_i$は対象者$i$の測定値に仮に与える共分散行列で、$A_i$は各時点の分散、$R_i(\alpha)$は作業相関行列です。相関パラメータ$\alpha$を使って同一人内の関連を表し、この行列を各人の測定値にどの重みを置くかという推定計算に用います。
「作業」と呼ぶのは、$R_i(\alpha)$を測定値の真の相関行列だと確定するためではなく、平均を効率よく推定するための候補として使うからです。ただし、候補を変えると重みが変わり、有限の標本では回帰係数自体も動きうるため、「相関指定は標準誤差だけに影響する」ということではない点にご留意ください。
そして、この違いは共分散構造の記事で扱う、測定時点間の関連の置き方ともつながります。GEEでは応答の同時分布を一つに決めずに集団平均の係数を求めるため、相関の指定と、その指定が不完全だった場合の分散推定を別々に解釈することが重要です。
独立・交換可能・AR(1)の違い
独立という指定では、同じ人の異なる時点間の作業相関を0に置きます。実際に測定値が独立だと主張する設定ではなく、同一人の値を一組としてロバスト分散を求める場合にも、係数を推定するための重みを独立型に置けます。
一方、交換可能な構造では異なる2時点なら間隔にかかわらず同じ相関$\alpha$を置きます。5時点なら時点0と1も0と4も同じ値になるため、時間が離れた測定ほど似なくなるデータでは、計算上の相関構造と観察される関連に差が残ります。
$$R_{st}^{\mathrm{ind}}=\mathbb{1}(s=t),\quad R_{st}^{\mathrm{ex}}=\begin{cases}1&s=t\\\alpha&s\ne t\end{cases},\quad R_{st}^{\mathrm{AR(1)}}=\alpha^{|s-t|}\tag{2}$$
式(2)は、独立型、交換可能型、一次自己回帰型の相関を時点$s,t$の組合せごとに示します。AR(1)型では等間隔の測定を前提に、時点差が1なら$\alpha$、4なら$\alpha^4$と置き、時点の実際の間隔が大きく異なる研究ではそのまま当てはめると適切な分析が行えなくなってしまいます。
今回のデータ生成では、標準偏差6の個人共通の差と、隣接時点で相関0.55の測定誤差を足しました。したがって、残差の0.55を測定値全体の時点差1の相関と読むのは誤りで、個人共通の差も含めて計算すると時点差1の相関は0.862、時点差4では0.720です。
この生成構造は交換可能型にも純粋なAR(1)型にも一致せず、作業相関は両方とも近似です。実際の推定で交換可能型の$\alpha$は0.779、AR(1)型の$\alpha$は0.885になりましたが、後者は個人共通の変動も含む作業相関であり、生成時の残差相関0.55とは推定対象が異なります。
二種類の標準誤差が表すもの
作業相関を含む$V_i$を正しい共分散とみなして求めるのが、モデルに基づく標準誤差です。その前提とデータの関連が合わない場合には、係数の標本間の揺れを小さく、または大きく見積もることがあり、点推定値の比較だけではこの違いを見落とします。
そこで、ロバスト分散では、各人について平均の予測からの差を求め、そのまとまりが人と人の間でどれほど変動するかを使います。対象者を独立な単位とし、一人の中の測定値が相関していても、その人の寄与をまとめて扱うため、作業相関が完全に合わない場合の分散推定に利用できます。
$$\widehat{\operatorname{Var}}_{\mathrm{rob}}(\widehat\beta)=B^{-1}\left\{\sum_{i=1}^{n}D_i^{\mathsf T}V_i^{-1}(Y_i-\widehat\mu_i)(Y_i-\widehat\mu_i)^{\mathsf T}V_i^{-1}D_i\right\}B^{-1},\quad B=\sum_{i=1}^{n}D_i^{\mathsf T}V_i^{-1}D_i\tag{3}$$
式(3)は、各人の残差の積から作る中央の部分を、平均と係数の関係を表す$B$で左右から挟むロバスト分散です。$D_i$は平均$\mu_i$の係数$\beta$による変化を表し、この形から、分散の計算単位が測定値1200個ではなく独立した人240人であることが分かります。
ただし、ロバスト分散が作業相関の誤指定に耐えるという性質は、平均モデルの設定、対象者間の独立性、独立した対象者数などの条件の下での漸近的な結果です。相関構造を変えると有限の標本では係数も標準誤差も動きうるため、今回のように両方を比較して報告します。
240人の再現例で見える差
再現例では対照112人と比較群128人をそれぞれ5時点測定し、群と時間の交互作用に真の係数3を置きました。GEEでは各群の平均的な時間変化の差をこの係数で表すため、以下の数値は特定の一人の変化速度ではなく、集団平均の差として読みます。
このデータでは、独立型と交換可能型の交互作用係数がともに2.959、AR(1)型では2.893となりました。前二者が一致したのは欠測のない均衡した時点配置と共変量の形による結果であり、別の測定配置では係数が異なり得ます。
一方、通常のモデルに基づく標準誤差は独立型0.282、交換可能型0.133、AR(1)型0.196と大きく異なります。同じ順序のロバスト標準誤差は0.180、0.180、0.171で、作業相関が異なるときの分散推定の変化を、実際の標本で確認できます。
この図の交換可能型では黒い区間が青い帯より狭く、独立型では逆に黒い区間が広がっています。この方向の違いは、ロバスト標準誤差を「保守的に大きくする補正」とだけ説明すると見落とす部分で、どの共分散をモデルに仮定したかによって比較結果が変わります。
また、青い区間は三つとも生成時の係数3を含みました。ただし、1回の標本で真の値を含んだ結果と被覆率は異なるため、後者を調べるには同じ条件のデータ生成を多数回繰り返し、区間が3を含む割合を求める必要があります。
ロバスト分散で直せない誤指定
平均が時間とともに曲がるデータに直線の平均モデルを当てはめた場合、ロバスト標準誤差だけを選び直しても、曲がり方を落とした平均の推定対象は変わりません。したがって、群と時間の交互作用を読む前に、時間の関数形、群ごとの平均の推移、残差の時点別の傾向を確かめます。
時間とともに変わる説明変数が過去の応答にも影響される場合には、さらに別の問題があります。PepeとAnderson（1994）は、そのような変数を含む周辺平均モデルで非対角の作業相関を使うと、通常の推定方程式に必要な条件が崩れうると論じており、相関だけを変えた分析で係数が動いたときは共変量の時間的な関係も確認します。
この再現例の説明変数は開始時点に決まる群と予定時点だけなので、時間とともに変わる説明変数の問題は入れていません。一般の研究で、以前の状態に応じて治療や測定値が変わる場合には、作業相関の選び方と、何を集団平均の効果として推定するかを切り離せません。
この欠測の問題も同様で、今回は全員から5回の値が得られるようにしました。過去の測定値によって次の測定の有無が変わると、観測された値だけを用いる通常の未加重GEEは偏りうるため、Robinsら（1995）が扱った逆確率重み付けなど、欠測の仕組みに対応した解析が必要になります。
そのため、ロバスト標準誤差で相関の誤指定を補えることと、欠測による選択や平均モデルの誤指定を補えることは区別します。途中の脱落がある研究では、時点別の解析人数と欠測の理由を示し、推定したい対象に対してどの仮定を置いたかを記載します。
独立した対象者が少ない場合
ロバスト分散の中央の項は、独立した人の寄与を足して求めます。240人を5回測った今回の例では独立な単位は240人ですが、同じ1200個の測定値でも24人を50回ずつ測った研究なら、独立な単位は24人に減ります。
この違いについてManclとDeRouen（2001）は、独立した集まりの数が少ないと通常のロバスト分散が小さくなり、検定の実際の偽陽性率が名目水準より高くなりうることを示し、少標本向けの補正を提案しました。頑健という名称だけで標本数の条件を省略せず、人や施設など実際に独立する単位の数を示す必要があります。
また、施設内の人が共通の影響を受ける研究では、人を独立単位とする式(3)の前提も変わります。個人内の反復だけをまとめるのか、施設まで同じ集まりとして扱うのかはデータの階層に依存するため、多階層データの記事で扱う単位を確認します。
公表研究と報告事項
LiangとZeger（1986）が作業相関とロバスト分散を含むGEEの枠組みを示した後、PanとConnett（2002）はLung Health Studyの肺機能と慢性的な咳嗽の解析で、作業相関の選択が推定の効率と時間依存共変量の扱いに関わることを検討しています。
このような観察研究を報告する際には、von Elmら（2007）のSTROBEが統計手法、欠測の扱い、解析人数の明示を求めます。GEEの結果を再現可能に示すには、それらに加えて、平均モデルの説明変数と時間の形、独立単位の定義、作業相関構造、標準誤差をモデルに基づく方法と頑健な方法のどちらで求めたかを記載します。
この比較では、同じ「95%信頼区間」でも分散の求め方が異なると幅が変わります。今回の交換可能型で0.133と0.180という標準誤差の差が生じたように、係数だけでなく区間の計算法まで示さなければ、研究間の結果を比べる条件がそろいません。
Rで係数と区間を再現する
今回のデータでは、対象者ごとの共通の変動と、近い時点ほど似る測定誤差を重ねています。以下のコードは同じ1200個の測定値を使い、作業相関だけを独立、交換可能、AR(1)に切り替えて係数と二種類の標準誤差を比べるものです。
set.seed(42)
library(geepack)
n <- 240; m <- 5; id <- rep(seq_len(n), each = m); time <- rep(0:4, n)
arm_i <- rbinom(n, 1, 0.5); arm <- arm_i[id]
b <- rnorm(n, sd = 6); phi <- 0.55; sigma <- 4
z <- matrix(rnorm(n * m), nrow = n); e <- matrix(NA_real_, n, m); e[, 1] <- sigma * z[, 1]
for (j in 2:m) e[, j] <- phi * e[, j - 1] +
sigma * sqrt(1 - phi^2) * z[, j]
y <- 100 + 2 * time + 3 * arm * time + b[id] +
as.vector(t(e))
d <- data.frame(id, time, arm, y)
fit <- function(cs) geeglm(y ~ arm * time, id = id, data = d, family = gaussian, corstr = cs, std.err = "san.se")
fits <- lapply(c("independence", "exchangeable", "ar1"), fit); names(fits) <- c("independence", "exchangeable", "ar1")このデータでは対象者を識別する番号ごとに5回の測定値を連続して並べ、各モデルへ同じ平均の式を当てはめました。群と時間の交互作用の係数を取り出し、モデルに基づく分散とロバスト分散の対角成分の平方根を比べます。
extract <- function(f) {
j <- match("arm:time", names(coef(f)))
c(beta = unname(coef(f)[j]),
model_se = sqrt(diag(f$geese$vbeta.naiv))[j],
robust_se = sqrt(diag(f$geese$vbeta))[j])
}
mat <- t(vapply(fits, extract, numeric(3))); round(mat, 3)#> beta model_se robust_se
#> independence 2.959 0.282 0.180
#> exchangeable 2.959 0.133 0.180
#> ar1 2.893 0.196 0.171この出力では独立型と交換可能型の係数・ロバスト標準誤差が一致しますが、モデルに基づく標準誤差は大きく異なります。また、AR(1)型では係数自体が2.893へ動いており、相関構造の指定と分散推定法の両方を報告する理由が数値から分かります。
png("gee-se-comparison.png", width = 1400, height = 950, res = 160)
par(mar = c(5, 10, 3, 2), mgp = c(2.8, .8, 0))
plot(NA, xlim = c(2, 4), ylim = c(.4, 3.8), yaxt = "n",
xlab = "Treatment-by-time coefficient (95% Wald interval)", ylab = "")
axis(2, at = 3:1, labels = names(fits), las = 1)
abline(v = 3, col = "gray70")
for (i in 1:3) {
yy <- 4 - i
rect(mat[i, 1] - 1.96 * mat[i, 3], yy - .10,
mat[i, 1] + 1.96 * mat[i, 3], yy + .10,
col = adjustcolor("#2B6E91", alpha.f = .25), border = NA)
points(mat[i, 1], yy, pch = 16, col = "#2B6E91", cex = 1.5)
segments(mat[i, 1] - 1.96 * mat[i, 2], yy - .13,
mat[i, 1] + 1.96 * mat[i, 2], yy - .13, col = "gray30", lwd = 1.5)
}
legend("topright", legend = c("Robust 95% CI", "Model-based 95% CI",
"Generating value"), col = c("#2B6E91", "gray30", "gray70"),
lty = c(NA, 1, 1), lwd = c(NA, 1.5, 1),
pch = c(15, NA, NA), pt.cex = c(1.8, NA, NA), bty = "n", cex = .8)
dev.off()青い半透明の帯はロバスト分散、黒い線は作業相関の仮定に基づく分散から計算しています。表示する区間の幅だけを比べるのではなく、どの仮定で得た区間かを区別して読むことが、GEEの結果報告には必要です。
作業相関の選択を結果へ結び付ける
今回の5時点では、AR(1)型が時間差に応じて相関を弱める一方、個人共通の変動のために遠い時点間にも相関が残ります。どの候補も生成構造に完全には一致しませんが、平均モデルが正しく、対象者が独立し、十分な人数がある状況では、ロバスト分散を併記して相関の誤指定による不確かさを扱えます。
したがって、結果には作業相関構造と分散推定法、独立単位の数をそろえて示します。群と時間の係数が相関指定を変えて大きく動く場合には、効率の差だけで片付けず、時間依存共変量、欠測、平均の形を確かめることが、数値の解釈を支えます。
参考文献
- Liang KY, Zeger SL. Longitudinal data analysis using generalized linear models. Biometrika. 1986;73(1):13–22.
- Mancl LA, DeRouen TA. A covariance estimator for GEE with improved small-sample properties. Biometrics. 2001;57(1):126–134.
- Pan W, Connett JE. Selecting the working correlation structure in generalized estimating equations with application to the Lung Health Study. Statistica Sinica. 2002;12(2):475–490.
- Pepe MS, Anderson GL. A cautionary note on inference for marginal regression models with longitudinal data and general correlated response data. Communications in Statistics – Simulation and Computation. 1994;23(4):939–951.
- von Elm E, Altman DG, Egger M, et al. The STROBE statement: guidelines for reporting observational studies. PLoS Medicine. 2007;4(10):e296.
- Robins JM, Rotnitzky A, Zhao LP. Analysis of semiparametric regression models for repeated outcomes in the presence of missing data. Journal of the American Statistical Association. 1995;90(429):106–121.
