同じ対象者から複数時点の有害事象の有無や受診回数を集めると、観測値同士は独立ではなくなり、通常のロジスティック回帰やポアソン回帰をそのまま当てはめれば係数は得られるものの、対象者内の似かよりを無視した標準誤差や、研究目的と異なる効果を読んでしまうおそれがあります。一般化線形混合モデルは、二値とカウントの分布を保ちながら、この依存をランダム効果で表す方法です。
そのため、この記事を通じて、二値縦断データにロジスティック混合モデル、件数の縦断データにポアソン混合モデルを使う理由と、固定効果を条件付きオッズ比または条件付き率比としての読み方の理解が進むと幸いです。さらに、ランダム切片が表せる依存、母集団平均効果との区別、過分散や観察時間への対応までを判断材料にします。
反復測定では独立性が崩れる
例えば、600人について5回ずつ治療反応の有無を記録すれば、観測数は3,000です。ただし、同じ人の5個の値は、体質や重症度など測定されていない対象者固有の特徴を共有するため、3,000個の独立な情報とは扱えません。
この依存を無視すると、係数だけでなく不確実性の評価にも影響します。そこで、一般化線形混合モデル（generalized linear mixed model: GLMM）では、応答の分布とリンク関数に加えて、対象者ごとのランダム効果をモデルに組み込みます。
GLMMは分布と個人差を結ぶ
BreslowとClayton（1993）が扱ったGLMMの枠組みでは、対象者$i$の時点$j$における応答を$Y_{ij}$、治療群の指標を$T_i$、時点を$t_{ij}$、対象者ごとのランダム切片を$b_i$とすると、基本形は次のように書けます。
$$
g\{\mathbb{E}(Y_{ij}\mid b_i)\}=\beta_0+\beta_1T_i+\beta_2t_{ij}+b_i,\qquad b_i\sim N(0,\sigma_b^2). \tag{1}
$$
式(1)は、同じ対象者の観測値が共通の$b_i$を持つことで、対象者内の似かよりを表すモデルです。$g$はリンク関数、$\beta_0$は切片、$\beta_1$は治療群の係数、$\beta_2$は時点の係数、$\sigma_b$はランダム切片の標準偏差です。
リンク関数を変えると、同じ式から二値とカウントを扱えます。二値応答にはロジットリンクと二項分布、件数には対数リンクとポアソン分布を組み合わせるため、別々の手法名を暗記するより、応答分布に合わせてGLMMの部品を選ぶ方が筋道を追いやすくなります。
二値データでは条件付きオッズを読む
ロジスティック混合モデルでは、$\exp(\beta_1)$は同じランダム切片を持つ対象者同士を比べた対象者別（subject-specific）の条件付きオッズ比です。今回のランダムに生成したデータセットでは、治療群の条件付きオッズ比は0.547、95%信頼区間は0.436〜0.684でした。
ただし、この0.547を母集団全体の平均的なオッズ比とは読めません。ロジット変換は非線形なので、ランダム効果を同じ値に固定した条件付きオッズ比と、ランダム効果を積分して得る周辺オッズ比は一般に一致しないからです。
この違いは単なる計算方法の差ではなく、推定対象の差です。
この違いについて、Zegerら（1988）は対象者別の条件付き効果と母集団平均効果を区別しており、Parzenら（2011）も通常のランダム効果ロジスティック回帰では周辺平均が同じロジット形を一般には保たないことを扱っていますが、この例の治療群は対象者内で変化しないため、0.547は同一対象者の治療時と対照時を比べた値ではありません。ランダム効果を同じ値に固定した条件付き効果が研究目的ならGLMM、集団への施策の平均的効果が目的ならGEEなどの母集団平均モデルというように、求める推定対象から枠組みを選ぶ必要があります。
カウントでは条件付き率比になる
再入院回数や有害事象件数を同じ期間で繰り返し記録する場合、ポアソン混合モデルの$\exp(\beta_1)$は、同じランダム切片を持つ対象者同士の対象者別の条件付き率比です。生成データでは治療群の条件付き率比が0.656、95%信頼区間が0.589〜0.732でしたが、ここでも時間不変の治療群については、同一対象者内の治療時と対照時の比較を意味しません。
一方、対象者によって観察日数が異なるなら、件数だけを比較してはいけません。対数リンクの線形予測子へ観察時間の対数を係数1のオフセットとして加え、単位時間当たりの発生率をモデル化します。
また、ポアソン分布は条件付き平均と条件付き分散が等しい形を置くため、過分散やゼロ過剰の確認が欠かせません。ThallとVail（1990）は縦断カウントデータで対象者内依存と過分散を扱っており、分散が大きいデータにポアソン分布を機械的に当てはめる進め方は適切ではありません。
ランダム切片が表す範囲
ランダム切片は、全時点を通じてイベントが起きやすい人、件数が多い人という対象者間の異質性を表します。今回のデータ生成では、標準偏差を二値モデル0.90、カウントモデル0.60としました。
もっとも、ランダム切片を入れただけで、近い時点ほど強く似る相関の減衰や、個人ごとに異なる時間傾向まで表せるわけではなく、ランダム傾きを加えるか、別の相関構造を選ぶかは測定計画とデータの形で変わり、この区別には共分散構造：compound symmetryからunstructuredまでの説明もつながります。
独立モデルとの差はバイアス量ではない
同じデータへ通常のロジスティック回帰とGLMMを当てはめたとき、得られる治療係数の差を、ランダム切片が除いたバイアスの量と呼ぶのは誤りです。
その理由は、独立モデルとGLMMで依存の扱いが異なるうえ、二値モデルでは条件付き効果と周辺効果という推定対象も異なるためです。係数が動いた事実だけで片方を正解とせず、研究目的に対応する効果と標準誤差を分けて比べます。
カウントデータでも、独立モデルに残る過分散がランダム切片による対象者間異質性だけに由来するとは限りません。GLMMを使った後も、分散統計量だけで分布仮定が成立したと判断せず、残差やモデル構造を検討します。
研究報告では推定対象まで書く
GLMMを使ったという記載だけでは、結果を再現できません。応答分布、リンク関数、固定効果、ランダム効果の構造、推定方法を示し、オッズ比や率比がランダム効果を同じ値に固定した対象者別の条件付き効果であることまで説明する必要があります。
この位置づけについて、BreslowとClayton（1993）は二項応答とポアソン応答を含むGLMMの近似推測を整理し、Batesら（2015）はlme4による混合モデルの推定実装を示しています。収束警告や特異適合が出た場合は数値だけを報告せず、ランダム効果構造、データ量、最適化の状態を確かめます。
この報告内容を考える時、私なら係数表より先に、研究目的が条件付き効果を求めているかを確認します。母集団平均効果を知りたいのに、ソフトウェアで使いやすいという理由だけでGLMMを選ぶ順序は推奨しません。
Rで二つの応答を比べる
600人を各5時点で観測し、治療群と対照群を各300人とした3,000観測のデータセットをランダムに生成しました。二値応答には治療オッズ比0.60、カウント応答には治療率比0.72を置き、対象者ごとのランダム切片を加えています。
set.seed(42)
n_subjects <- 600L; subject <- factor(rep(seq_len(n_subjects), each = 5))
visit <- rep(0:4, times = n_subjects); treatment_subject <- sample(rep(c(0L, 1L), each = n_subjects / 2L))
treatment <- rep(treatment_subject, each = 5)
binary_random <- rnorm(n_subjects, mean = 0, sd = 0.90)
binary_eta <- -1 + log(0.60) * treatment + log(1.15) * visit +
binary_random[as.integer(subject)]
binary_outcome <- rbinom(length(subject), size = 1, prob = plogis(binary_eta))
count_random <- rnorm(n_subjects, mean = 0, sd = 0.60)
count_eta <- log(2.20) + log(0.72) * treatment + log(1.10) * visit +
count_random[as.integer(subject)]
count_outcome <- rpois(length(subject), lambda = exp(count_eta))解析にはR 4.4.2とlme4 1.1-38を使いました。治療群の係数と95%信頼区間は次のコードと一致します。
analysis_data <- data.frame(subject, visit, treatment,
binary_outcome, count_outcome)
binary_fit <- lme4::glmer(binary_outcome ~ treatment + visit + (1 | subject),
data = analysis_data, family = binomial(link = "logit"), nAGQ = 9)
count_fit <- lme4::glmer(count_outcome ~ treatment + visit + (1 | subject),
data = analysis_data, family = poisson(link = "log"), nAGQ = 9)
binary_ci <- exp(confint(binary_fit, parm = "treatment", method = "Wald"))
count_ci <- exp(confint(count_fit, parm = "treatment", method = "Wald"))
binary_row <- setNames(c(unname(exp(lme4::fixef(binary_fit)["treatment"])), binary_ci), c("estimate", "lower", "upper"))
count_row <- setNames(c(unname(exp(lme4::fixef(count_fit)["treatment"])), count_ci), c("estimate", "lower", "upper"))
round(rbind(binary_or = binary_row, count_irr = count_row), 3)#> estimate lower upper
#> binary_or 0.547 0.436 0.684
#> count_irr 0.656 0.589 0.732出力の読み方は、同じランダム切片と同じ時点を持つ対象者同士で、治療群のオッズまたは平均件数が低いというものです。読み違えやすい点は、これらを母集団平均の確率差や件数差へ置き換えられず、p値だけから医学的な重要性も決められないことです。
さらに、同じモデルから時点ごとの条件付き予測値を描くと、二値確率と平均件数のどちらでも群間差が保たれながら、時点とともに値が上がる設定を確認できます。
prediction <- expand.grid(visit = 0:4, treatment = 0:1)
prediction$binary <- predict(binary_fit, newdata = prediction,
type = "response", re.form = NA)
prediction$count <- predict(count_fit, newdata = prediction,
type = "response", re.form = NA)
png("glmm-conditional-trajectories.png", width = 1100, height = 688,
res = 200, type = "cairo")
par(mfrow = c(1, 2), las = 1, mar = c(4.2, 4.2, 1.5, 2.2), cex.lab = 0.8, cex.axis = 0.8)
colors <- c("#3568A8", "#D05A4E")
for (outcome in c("binary", "count")) {
matplot(0:4, matrix(prediction[[outcome]], nrow = 5), type = "o",
col = colors, lty = 1, pch = c(16, 17), xlab = "Visit",
ylab = if (outcome == "binary") "Conditional probability" else "Conditional mean count")
legend("topleft", legend = c("Control", "Treatment"), col = colors,
pch = c(16, 17), lty = 1, bty = "n", cex = 0.7)
}
dev.off()
選択は応答より推定対象から始まる
二値ならロジスティック、件数ならポアソンという対応だけでは、縦断データの解析方針は決まりません。条件付き効果と母集団平均効果のどちらが研究目的に合うか、ランダム切片で依存を十分に表せるか、カウントでは観察時間と過分散をどう扱うかが選択を左右します。
実際のデータでは、ゼロ過剰、欠測、ランダム傾きの要否も結果へ影響します。こうした条件を含めて解析方針を決める必要がある場合は、Dr.データサイエンスで、研究目的とデータ構造に沿った統計解析の相談を受け付けています。
参考文献
- Breslow NE, Clayton DG. Approximate Inference in Generalized Linear Mixed Models. Journal of the American Statistical Association. 1993;88(421):9-25.
- Zeger SL, Liang KY, Albert PS. Models for Longitudinal Data: A Generalized Estimating Equation Approach. Biometrics. 1988;44(4):1049-1060.
- Parzen M, Ghosh S, Lipsitz S, Sinha D, Fitzmaurice GM, Mallick BK, Ibrahim JG. A Generalized Linear Mixed Model for Longitudinal Binary Data with a Marginal Logit Link Function. The Annals of Applied Statistics. 2011;5(1):449-467.
- Thall PF, Vail SC. Some Covariance Models for Longitudinal Count Data with Overdispersion. Biometrics. 1990;46(3):657-671.
- Bates D, Mächler M, Bolker B, Walker S. Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software. 2015;67(1):1-48.
