研究者: 反復測定の2値データでGLMM（
glmer）を当てはめました。査読対策に「ロバスト分散」も出しておこうと思ってsandwich::vcovCL()を使ったら、エラーで止まってしまいました……。そもそも、ロバスト分散を使えばGLMMの仮定が崩れていても大丈夫になるんでしょうか？
統計ER: 実は
vcovCL()はglmerの結果に対応していないんです。merDerivパッケージを足せば計算できますよ。それと、ロバスト分散は「何でも直してくれる万能薬」ではありません。守ってくれるもの・守ってくれないものを整理してから使いましょう！
個体差や反復測定などでデータ間に相関がある場合、一般化線形混合モデル（GLMM）は強力な道具である。そして、モデルの仮定が完全には正しくない場合に備えて、標準誤差を「データ自身のばらつき」から計算し直すのがロバスト分散（サンドイッチ分散）である。
本記事では、GLMMにロバスト分散を組み合わせる意味と限界を整理したうえで、Rで実際に計算する手順（lme4 + merDeriv + sandwich）と、結果の読み方、論文での書き方までを解説する。
GLMMの標準誤差が「自信過剰」になるとき
GLMMは、一般化線形モデル（ロジスティック回帰など）に変量効果を加え、同じ患者から得られたデータ同士の相関を表現するモデルである。例えば、患者ごとの「改善しやすさ」の違いを変量切片として入れると、同じ患者の測定値同士が似通っていることをモデルが理解できる。
ただし、GLMMの標準誤差（モデルベースの標準誤差）は「モデルが正しい」という前提で計算されている。特に次の部分がずれていると、標準誤差は実際のばらつきより小さく（あるいは大きく）見積もられうる。
- 相関構造（変量効果の構造）：例えば、患者ごとに改善の「スピード」も違う（変量傾きが必要）のに、変量切片だけのモデルを当てはめている
- 変量効果の分布：変量効果は正規分布に従うと仮定しているが、実際には歪んでいたり、2つの集団が混ざっていたりする
- 応答の分布：二項分布・ポアソン分布の仮定よりも、実際のばらつきが大きい（過分散）
モデルが自己申告する標準誤差は、いわば天気予報の「降水確率の自己評価」である。予報モデルが正しければ自己評価は当てになるが、モデルが現実とずれていれば、実績（実際に雨が降ったかどうか）で答え合わせをしたほうが確かだ。ロバスト分散は、この「実績による答え合わせ」を標準誤差に対して行う方法である。
ロバスト分散（サンドイッチ推定）の考え方
ロバスト分散は、Huber–White のサンドイッチ推定量とも呼ばれる（White, 1982）。パラメータ推定値の分散を次の形で計算する。
$$\widehat{\mathrm{Var}}(\hat{\theta}) = A^{-1}\, B\, A^{-1}, \qquad B = \sum_{i=1}^{N} s_i(\hat{\theta})\, s_i(\hat{\theta})^{\top}$$
- $A$（パン）：対数尤度の曲がり具合（負のヘッセ行列）。モデルが「このくらい正確に推定できるはず」と考えている情報量
- $s_i$：患者 $i$ のスコア（対数尤度の傾き）を、その患者の全測定についてまとめたもの
- $B$（具）：患者ごとのスコアが実際にどれだけばらついたか
モデルが正しければ $A$ と $B$ はほぼ等しくなり、サンドイッチ分散は通常の分散（$A^{-1}$）とほぼ一致する。モデルがずれていると両者が食い違い、その食い違いを $B$ が実データから補正する。パンでモデルの想定を挟み、具に現実を入れる、というイメージである。
※コラム：「対数尤度の曲がり具合（負のヘッセ行列）」とは
最尤法は、対数尤度という「山」の頂上（最も当てはまりのよいパラメータの値）を探す作業である。頂上が鋭くとがった山なら、少し横にずれただけで高さが大きく下がるので、頂上の位置ははっきり決まる。逆に、なだらかな丘のような山だと、頂上付近はどこも同じような高さで、頂上の位置はあいまいになる。霧の中で山頂の旗を探すとき、とがった山なら迷わないが、平らな丘では見当がつきにくいのと同じである。
この「とがり具合」を数値にしたものが、対数尤度の2階微分（曲率）である。パラメータが複数あるので、2階微分を並べた行列になり、これをヘッセ行列と呼ぶ。山の頂上では上に凸なので2階微分は負になる。そこで符号を反転させた負のヘッセ行列（観測情報行列とも呼ぶ）を「情報量」として使う。情報量が大きい（山がとがっている）ほど推定は正確で、その逆行列 $A^{-1}$ が、
summary()が通常表示する標準誤差のもとになっている。ただし、この計算は「山の形をモデルが正しく描けている」ことが前提である。サンドイッチ推定では、この $A^{-1}$ をパンとして使いつつ、具 $B$ で実際のばらつきを入れて補正する。
※コラム：なぜ「患者単位」でまとめるのか
ロバスト分散では、スコアを測定1回ごとではなく患者（クラスター）ごとに合計してからばらつきを見る。こうすると、同じ患者の中の相関がどんな形であっても、その影響が $B$ にまるごと入る。だからこそ「相関構造の指定ミス」に強い。その代わり、ばらつきを見積もる単位は「患者の数」になるため、患者数（クラスター数）が少ないと $B$ 自体が不安定になる。
ロバスト分散が守ってくれるもの・守ってくれないもの
ここが最も誤解されやすい点である。ロバスト分散が直すのは標準誤差（＝信頼区間の幅とP値）だけで、推定値そのものは1ミリも変わらない。
| 状況 | ロバスト分散で救えるか | 補足 |
|---|---|---|
| 変量効果の構造の指定ミス（変量傾きの入れ忘れなど） | 標準誤差は救える | 推定値はそのまま。本来は構造を見直すのが本筋 |
| 患者内の過分散 | 標準誤差は救える | — |
| 変量効果の分布の形のずれ（正規でない） | 標準誤差は救える。推定値の偏りは救えない | ロジスティックGLMMの回帰係数は、分布の形のずれに比較的強いことが知られる（McCulloch & Neuhaus, 2011）。切片や変量効果の分散は影響を受けやすい |
| 平均構造の誤り（交絡因子の入れ忘れ、リンク関数の誤りなど） | 救えない | 推定値そのものが偏っている |
| 患者数（クラスター数）が少ない | むしろ悪化しうる | サンドイッチ分散は過小になりやすい（Mancl & DeRouen, 2001） |
つまり、ロバスト分散は「推定値が狙っている量はおおむね正しいが、その不確かさの見積もりに自信がない」ときの保険である。偏った推定値を正しくする道具ではない。
※コラム：GLMMの「条件付き」効果とGEEの「周辺」効果
ロバスト分散の本家は、一般化推定方程式（GEE; Liang & Zeger, 1986）である。GEEは「集団全体で平均したオッズ比」（周辺効果、population-averaged）を推定し、標準誤差は最初からサンドイッチで計算する。一方、GLMMの回帰係数は「同じ変量効果をもつ患者の中でのオッズ比」（条件付き効果、subject-specific）である。ロジスティックモデルでは、周辺オッズ比は条件付きオッズ比より1に近くなる（減衰する）ため、GLMM＋ロバスト分散とGEEは、似ているが別の量を推定している。結果を比べるときは、この違いを念頭に置きたい。
通常の標準誤差とロバスト分散の使い分け
- 主解析：事前に計画したGLMM（モデルベースの標準誤差）を基本とし、変量効果の構造は、研究デザインと臨床的な知識から妥当なものを選ぶ
- 感度分析：ロバスト分散でも信頼区間を計算し、主解析と大きく変わらないかを確認する。大きく変わる場合は、まず変量効果の構造（変量傾きの必要性など）を見直す
- クラスター数が少ないとき（目安として数十未満）：素朴なサンドイッチ分散は小さく出やすい。小標本補正（Mancl–DeRouen補正など）を検討するか、ロバスト分散に頼りすぎない
- 集団平均の効果を知りたいとき：はじめからGEE（または周辺効果を出せるモデル）を選ぶ
「とにかくロバスト分散を使えば安心」ではなく、「主解析＋ロバスト分散で答え合わせ」という位置づけで使うのが実践的である。
具体例：治療法・性別と経時的な病状改善
架空の臨床試験を考える。ある疾患の患者100人を治療法A・Bにランダムに割り付け、治療開始時から1週間ごとに計4回（visit = 0, 1, 2, 3）、病状が「改善」したかどうか（1 = 改善、0 = 不変・悪化）を評価した。性別（Female / Male）も記録されている。
ここで、真のデータ生成の仕組みでは、患者ごとに「改善しやすさ」（変量切片）だけでなく「改善のスピード」（変量傾き）も異なるとする。しかし解析では、よく使われる変量切片だけのGLMMを当てはめてみる。つまり、相関構造をあえて誤指定した状況である。
患者 $i$ の $j$ 回目の評価で改善する確率を $p_{ij}$ とすると、当てはめるモデルは次のとおりである。
$$\mathrm{logit}(p_{ij}) = \beta_0 + \beta_1\,\mathrm{TreatB}_i + \beta_2\,\mathrm{Male}_i + \beta_3\,(\mathrm{TreatB}_i \times \mathrm{Male}_i) + \beta_4\,\mathrm{visit}_{ij} + u_i$$
- $\beta_0 \sim \beta_4$：固定効果の係数（治療法・性別は患者単位、visitは測定ごとに変わる共変量）
- $u_i \sim N(0, \sigma_u^2)$：患者 $i$ の変量切片
真の仕組みには患者ごとの傾き $u_{1i} \times \mathrm{visit}_{ij}$ もあるが、このモデルには入っていない。この「入れ忘れ」が、どの係数の標準誤差に効いてくるかを見ていく。
Rでの計算例
GLMMは lme4::glmer() で当てはめ、ロバスト分散は merDeriv パッケージ（Wang et al., 2022）と sandwich パッケージで計算する。
動作確認環境：R 4.6.1、lme4 2.0.6、merDeriv 0.2.6、sandwich 3.1.3、GLMMadaptive 0.9.7、geepack 1.3.13、dplyr 1.2.1（2026年10月時点）。lme4 のバージョンによって、出力の小数点以下3〜4桁目がわずかに変わることがある。
※注意：よくあるつまずき
sandwich::vcovCL()はglmerの結果（glmerMod）に対応しておらず、「no applicable method for ‘estfun’ applied to an object of class glmerMod」というエラーで止まる。また、clubSandwich::vcovCR()も glmerMod には対応していない（線形混合モデルのlmerには対応）。merDerivを読み込むと、glmerMod 用のestfun()（患者ごとのスコア）とbread()が追加され、sandwich::sandwich()で計算できるようになる。
Step 0：パッケージの準備
pkgs <- c("lme4", "merDeriv", "sandwich", "GLMMadaptive", "geepack", "dplyr")
for (p in pkgs) if (!requireNamespace(p, quietly = TRUE)) install.packages(p)
library(lme4)
library(merDeriv) # glmerMod 用の estfun() / bread() を追加する
library(sandwich)
library(dplyr)merDeriv を読み込むと、依存パッケージ（lavaan、nonnest2）の起動メッセージが表示されるが、問題ない。
Step 1：疑似データの作成
# 100人 × 4回（visit = 0, 1, 2, 3）の2値アウトカム（1 = 改善）
# 真のモデルでは「改善のスピード」が患者ごとに違う（変量傾き u1）
# ------------------------------------------------------------
set.seed(123)
n_patients <- 100
n_visits <- 4
patients <- data.frame(
patient_id = 1:n_patients,
treatment = factor(sample(c("A", "B"), n_patients, replace = TRUE), levels = c("A", "B")),
sex = factor(sample(c("Female", "Male"), n_patients, replace = TRUE), levels = c("Female", "Male")),
u0 = rnorm(n_patients, 0, 0.8), # 患者ごとの「改善しやすさ」（変量切片）
u1 = rnorm(n_patients, 0, 0.6) # 患者ごとの「改善スピード」（変量傾き）
)
dat <- expand.grid(patient_id = 1:n_patients, visit = 0:(n_visits - 1)) %>%
left_join(patients, by = "patient_id") %>%
mutate(
logit_p = -1 + 0.5 * (treatment == "B") - 0.4 * (sex == "Male") +
0.3 * (treatment == "B") * (sex == "Male") +
0.3 * visit + u0 + u1 * visit,
improvement = rbinom(n(), 1, plogis(logit_p)),
patient_id = factor(patient_id)
) %>%
arrange(patient_id, visit) # GEE用に患者ごとに並べておく性別は Female を基準カテゴリに明示している（levels を指定しないとアルファベット順で決まるため、意図しない基準になりやすい）。
Step 2：変量切片だけのGLMMを当てはめる
model_glmm <- glmer(
improvement ~ treatment * sex + visit + (1 | patient_id),
data = dat, family = binomial(link = "logit")
)
summary(model_glmm)$coefficients Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.46115234 0.28185590 -1.6361280 0.10181280
treatmentB -0.26490874 0.36644416 -0.7229171 0.46973083
sexMale -0.41883313 0.32899818 -1.2730561 0.20299814
visit 0.22209649 0.09740672 2.2800941 0.02260211
treatmentB:sexMale 0.06461201 0.50454949 0.1280588 0.89810242真の値（treatmentB = 0.5 など）と推定値がずれているのは、患者100人という規模での偶然誤差である（標準誤差は約0.37あり、この程度のずれは偶然でも起こりうる。別のシードで実行すると推定値は真の値に近づいたり離れたりする）。
Step 3：患者単位のロバスト分散を計算する
p <- length(fixef(model_glmm))
vcov_robust_all <- sandwich(model_glmm,
bread = bread(model_glmm, full = TRUE),
meat = meat(model_glmm))
vcov_robust <- vcov_robust_all[1:p, 1:p] # 固定効果の部分だけ取り出す
se_model <- sqrt(diag(as.matrix(vcov(model_glmm))))
se_robust <- sqrt(diag(vcov_robust))bread(..., full = TRUE) は固定効果に加えて変量効果の分散も含めた行列を返すため、最後に固定効果の部分（左上の $p \times p$）だけを取り出している。
Step 4：通常SEとロバストSEを比較する
est <- fixef(model_glmm)
z <- qnorm(0.975)
comparison <- data.frame(
term = names(est),
estimate = est,
SE_model = se_model,
SE_robust = se_robust,
SE_ratio = se_robust / se_model,
p_model = 2 * pnorm(-abs(est / se_model)),
p_robust = 2 * pnorm(-abs(est / se_robust)),
row.names = NULL
)
print(comparison, digits = 3)
comparison %>%
filter(term != "(Intercept)") %>%
transmute(
term,
OR_model_CI = sprintf("%.2f (%.2f-%.2f)", exp(estimate),
exp(estimate - z * SE_model), exp(estimate + z * SE_model)),
OR_robust_CI = sprintf("%.2f (%.2f-%.2f)", exp(estimate),
exp(estimate - z * SE_robust), exp(estimate + z * SE_robust))
) term estimate SE_model SE_robust SE_ratio p_model p_robust
1 (Intercept) -0.4612 0.2819 0.256 0.909 0.1018 0.0720
2 treatmentB -0.2649 0.3664 0.378 1.031 0.4697 0.4831
3 sexMale -0.4188 0.3290 0.319 0.970 0.2030 0.1892
4 visit 0.2221 0.0974 0.115 1.176 0.0226 0.0525
5 treatmentB:sexMale 0.0646 0.5045 0.513 1.016 0.8981 0.8997
term OR_model_CI OR_robust_CI
1 treatmentB 0.77 (0.37-1.57) 0.77 (0.37-1.61)
2 sexMale 0.66 (0.35-1.25) 0.66 (0.35-1.23)
3 visit 1.25 (1.03-1.51) 1.25 (1.00-1.56)
4 treatmentB:sexMale 1.07 (0.40-2.87) 1.07 (0.39-2.91)Step 5：別ルートで確認する
計算が正しいかを、別のパッケージでも確かめておく。GLMMadaptive は vcov(..., sandwich = TRUE) だけでロバスト分散を出せる。あわせて、変量傾きを入れたGLMMと、GEEも当てはめてみる。
library(GLMMadaptive)
model_ga <- mixed_model(improvement ~ treatment * sex + visit,
random = ~ 1 | patient_id,
data = dat, family = binomial())
cbind(estimate = fixef(model_ga),
SE_model = sqrt(diag(vcov(model_ga)))[1:p],
SE_robust = sqrt(diag(vcov(model_ga, sandwich = TRUE)))[1:p])
model_slope <- glmer(
improvement ~ treatment * sex + visit + (1 + visit | patient_id),
data = dat, family = binomial(link = "logit")
)
cbind(estimate = fixef(model_slope),
SE_model = sqrt(diag(as.matrix(vcov(model_slope)))))
library(geepack)
model_gee <- geeglm(improvement ~ treatment * sex + visit,
id = patient_id, data = dat,
family = binomial, corstr = "exchangeable")
summary(model_gee)$coefficients# GLMMadaptive（変量切片のみ）
estimate SE_model SE_robust
(Intercept) -0.46712474 0.29004902 0.2623063
treatmentB -0.26804277 0.38019809 0.3892545
sexMale -0.42419455 0.34123394 0.3285430
visit 0.22477266 0.09812843 0.1156556
treatmentB:sexMale 0.06515193 0.52339361 0.5276670
# 変量傾きも入れたGLMM
estimate SE_model
(Intercept) -0.5272414 0.3247425
treatmentB -0.1164794 0.4320546
sexMale -0.4322640 0.3835457
visit 0.2178499 0.1506007
treatmentB:sexMale -0.1250080 0.5893599
# GEE（交換可能相関、ロバストSE）
Estimate Std.err Wald Pr(>|W|)
(Intercept) -0.4216020 0.2334879 3.26044423 0.07096990
treatmentB -0.2457650 0.3530362 0.48462015 0.48633700
sexMale -0.3826376 0.2984716 1.64349818 0.19984618
visit 0.2034903 0.1032775 3.88218036 0.04880114
treatmentB:sexMale 0.0602403 0.4781236 0.01587427 0.89973735結果の解釈
- 推定値（オッズ比）は通常SEでもロバストSEでも同じである。変わるのは標準誤差・信頼区間・P値だけである
- 患者単位の共変量（治療法・性別・交互作用）では、ロバストSEと通常SEの比は0.97〜1.03で、ほとんど差がない
- 測定ごとに変わる共変量（visit）では、ロバストSEが通常SEの約1.18倍になった。変量傾きの入れ忘れは「時間とともに患者ごとに違う動き」を無視することなので、その影響は時間の係数の標準誤差に表れる
visit のオッズ比は1.25で、95%信頼区間は通常SEで1.03〜1.51、ロバストSEで1.00〜1.56となった。P値は0.023から0.053に変わり、ちょうど0.05をまたいだ。ここで「有意でなくなったから効果はない」と読むのは誤りである。点推定（1回の評価ごとに改善のオッズが約1.25倍）は同じで、その不確かさが当初の見積もりより大きかったというのが正しい読み方である。
筆者が同じ設定で300回データを生成し直して確認したところ、visit の係数の実際のばらつき（推定値の標準偏差）は0.109で、通常SEの平均（0.100）はやや小さく、ロバストSEの平均（0.112）はほぼ一致した。今回の例のずれは「1割程度の過小評価」で、ロバスト分散が補正しているのはまさにこの部分である。
別ルート (Step 5) の結果の読み方
- GLMMadaptive：通常SE・ロバストSEとも
merDerivとほぼ同じ値になった（visit のロバストSE：0.115 vs 0.116）。lme4はラプラス近似、GLMMadaptiveは適応的ガウス求積で尤度を計算しているため小数点以下にわずかな差が出るが、計算手順そのものは正しいと確認できる - 変量傾きを入れたGLMM：visit の標準誤差は0.151と、ロバストSE（0.115）よりさらに大きくなった。これは「患者ごとに傾きが違う中での平均的な傾き」という少し意味の異なる量を推定しているためで、単純な大小比較はできない。それでも、ロバストSEとモデルベースSEの食い違いは変量効果の構造を見直すサインとして使える
- GEE：推定値はGLMMより全体的に0に近い（visit：0.222 → 0.203、性別：−0.419 → −0.383）。これは前述の「条件付き効果と周辺効果の違い」によるもので、どちらかが誤っているわけではない。GEEの標準誤差（visit：0.103）は最初からロバスト分散で計算されている
論文でそのまま使える英文テンプレート
以下の数値は本記事の疑似データによる記述例である。実際の論文では、自分のデータの結果に置き換えてほしい。
Methods
To account for within-patient correlation of repeated binary outcomes, we fitted a logistic generalized linear mixed model with a patient-specific random intercept using the lme4 package in R. Fixed effects included treatment, sex, their interaction, and visit. As a sensitivity analysis for potential misspecification of the random-effects structure, we additionally calculated cluster-robust (sandwich) standard errors clustered at the patient level using the merDeriv and sandwich packages, and compared the resulting 95% confidence intervals with the model-based ones.
Results
Robust standard errors were similar to model-based standard errors for patient-level covariates (ratio, 0.97–1.03). For visit, the robust standard error was approximately 18% larger than the model-based one; the odds ratio per visit was 1.25 (95% CI, 1.03–1.51 with model-based standard errors; 1.00–1.56 with robust standard errors).
クラスター数が少ない研究では、Methodsに「小標本補正を行った」「クラスター数が少ないため、ロバスト分散の結果は参考として扱った」などの一文を加えるとよい。
まとめ
- ロバスト分散（サンドイッチ分散）は、患者単位のスコアの実際のばらつきから標準誤差を計算し直す方法である
- 守ってくれるのは標準誤差だけ。推定値は変わらず、平均構造の誤りや交絡による偏りは直せない
- 変量効果の構造の誤指定では、測定ごとに変わる共変量（時間など）の標準誤差に差が出やすい
- クラスター数が少ないと、ロバスト分散はかえって小さく出やすい
- Rでは
sandwich::vcovCL()やclubSandwichは glmer に対応していない。merDeriv+sandwich::sandwich()、またはGLMMadaptiveのvcov(..., sandwich = TRUE)を使う - 主解析のGLMMに対する感度分析（答え合わせ）として使い、大きな食い違いがあれば変量効果の構造を見直す
一括コピペ用Rスクリプト
# ============================================================
# ロバスト分散（サンドイッチ分散）を用いたGLMM
# lme4::glmer + merDeriv + sandwich で計算し、
# GLMMadaptive と geepack（GEE）で結果を確認する
# ============================================================
pkgs <- c("lme4", "merDeriv", "sandwich", "GLMMadaptive", "geepack", "dplyr")
for (p in pkgs) if (!requireNamespace(p, quietly = TRUE)) install.packages(p)
library(lme4)
library(merDeriv) # glmerMod 用の estfun() / bread() を追加する
library(sandwich)
library(dplyr)
# ------------------------------------------------------------
# Step 1. 疑似データの作成
# 100人 × 4回（visit = 0, 1, 2, 3）の2値アウトカム（1 = 改善）
# 真のモデルでは「改善のスピード」が患者ごとに違う（変量傾き u1）
# ------------------------------------------------------------
set.seed(123)
n_patients <- 100
n_visits <- 4
patients <- data.frame(
patient_id = 1:n_patients,
treatment = factor(sample(c("A", "B"), n_patients, replace = TRUE), levels = c("A", "B")),
sex = factor(sample(c("Female", "Male"), n_patients, replace = TRUE), levels = c("Female", "Male")),
u0 = rnorm(n_patients, 0, 0.8), # 患者ごとの「改善しやすさ」（変量切片）
u1 = rnorm(n_patients, 0, 0.6) # 患者ごとの「改善スピード」（変量傾き）
)
dat <- expand.grid(patient_id = 1:n_patients, visit = 0:(n_visits - 1)) %>%
left_join(patients, by = "patient_id") %>%
mutate(
logit_p = -1 + 0.5 * (treatment == "B") - 0.4 * (sex == "Male") +
0.3 * (treatment == "B") * (sex == "Male") +
0.3 * visit + u0 + u1 * visit,
improvement = rbinom(n(), 1, plogis(logit_p)),
patient_id = factor(patient_id)
) %>%
arrange(patient_id, visit) # GEE用に患者ごとに並べておく
# ------------------------------------------------------------
# Step 2. よく使われる「変量切片だけ」のGLMMを当てはめる
# （真の構造にある変量傾きは、あえて入れていない＝誤指定）
# ------------------------------------------------------------
model_glmm <- glmer(
improvement ~ treatment * sex + visit + (1 | patient_id),
data = dat, family = binomial(link = "logit")
)
summary(model_glmm)$coefficients
# ------------------------------------------------------------
# Step 3. 患者（クラスター）単位のロバスト分散を計算する
# ※ sandwich::vcovCL() は glmerMod に未対応でエラーになる
# ※ clubSandwich::vcovCR() も glmerMod には未対応
# ------------------------------------------------------------
p <- length(fixef(model_glmm))
vcov_robust_all <- sandwich(model_glmm,
bread = bread(model_glmm, full = TRUE),
meat = meat(model_glmm))
vcov_robust <- vcov_robust_all[1:p, 1:p] # 固定効果の部分だけ取り出す
se_model <- sqrt(diag(as.matrix(vcov(model_glmm))))
se_robust <- sqrt(diag(vcov_robust))
# ------------------------------------------------------------
# Step 4. 通常SEとロバストSE、オッズ比と95%信頼区間を比較する
# ------------------------------------------------------------
est <- fixef(model_glmm)
z <- qnorm(0.975)
comparison <- data.frame(
term = names(est),
estimate = est,
SE_model = se_model,
SE_robust = se_robust,
SE_ratio = se_robust / se_model,
p_model = 2 * pnorm(-abs(est / se_model)),
p_robust = 2 * pnorm(-abs(est / se_robust)),
row.names = NULL
)
print(comparison, digits = 3)
comparison %>%
filter(term != "(Intercept)") %>%
transmute(
term,
OR_model_CI = sprintf("%.2f (%.2f-%.2f)", exp(estimate),
exp(estimate - z * SE_model), exp(estimate + z * SE_model)),
OR_robust_CI = sprintf("%.2f (%.2f-%.2f)", exp(estimate),
exp(estimate - z * SE_robust), exp(estimate + z * SE_robust))
)
# ------------------------------------------------------------
# Step 5. 別ルートで確認する
# (a) GLMMadaptive：vcov(..., sandwich = TRUE) で直接計算できる
# (b) 変量傾きも入れたGLMM（構造を正しく指定した場合）
# (c) GEE：周辺（集団平均）モデル。ロバスト分散が標準
# ------------------------------------------------------------
library(GLMMadaptive)
model_ga <- mixed_model(improvement ~ treatment * sex + visit,
random = ~ 1 | patient_id,
data = dat, family = binomial())
cbind(estimate = fixef(model_ga),
SE_model = sqrt(diag(vcov(model_ga)))[1:p],
SE_robust = sqrt(diag(vcov(model_ga, sandwich = TRUE)))[1:p])
model_slope <- glmer(
improvement ~ treatment * sex + visit + (1 + visit | patient_id),
data = dat, family = binomial(link = "logit")
)
cbind(estimate = fixef(model_slope),
SE_model = sqrt(diag(as.matrix(vcov(model_slope)))))
library(geepack)
model_gee <- geeglm(improvement ~ treatment * sex + visit,
id = patient_id, data = dat,
family = binomial, corstr = "exchangeable")
summary(model_gee)$coefficientsおすすめ書籍
- 久保拓弥『データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC』岩波書店, 2012. GLMからGLMMへの考え方を、Rの例とともに直感的に学べる定番書。 Amazonで見る
- 丹後俊郎『経時的繰り返し測定デザイン―治療効果を評価する混合効果モデルとその周辺―』（医学統計学シリーズ）朝倉書店, 2015. 臨床試験の反復測定データを混合効果モデルで解析する考え方を医学統計の立場から解説。 Amazonで見る
- Fitzmaurice GM, Laird NM, Ware JH. Applied Longitudinal Analysis. 2nd ed. Wiley; 2011. GLMMとGEE（条件付き効果と周辺効果）の違い、ロバスト分散の使いどころを丁寧に扱う縦断データ解析の標準的教科書。 Amazonで見る
- Diggle P, Heagerty P, Liang KY, Zeger S. Analysis of Longitudinal Data. 2nd ed. Oxford University Press; 2002（ペーパーバック版 2013）. GEEの提唱者らによる古典。サンドイッチ分散の理論的背景を深く知りたい人向け。 Amazonで見る
参考文献
- White H. A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity. Econometrica. 1980;48(4):817–838. https://doi.org/10.2307/1912934
- White H. Maximum likelihood estimation of misspecified models. Econometrica. 1982;50(1):1–25. https://doi.org/10.2307/1912526
- Liang KY, Zeger SL. Longitudinal data analysis using generalized linear models. Biometrika. 1986;73(1):13–22. https://doi.org/10.2307/2336267
- Mancl LA, DeRouen TA. A covariance estimator for GEE with improved small-sample properties. Biometrics. 2001;57(1):126–134. 10.1111/j.0006-341x.2001.00126.x
- McCulloch CE, Neuhaus JM. Misspecifying the shape of a random effects distribution: why getting it wrong may not matter. Statistical Science. 2011;26(3):388–402. https://doi.org/10.1214/11-STS361
- Wang T, Graves B, Rosseel Y, Merkle EC. Computation and application of generalized linear mixed model derivatives using lme4. Psychometrika. 2022;87(3):1173–1193. https://doi.org/10.1007/s11336-022-09840-2
- 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. https://ideas.repec.org/a/jss/jstsof/v067i01.html
コメント