一般化線形混合モデル(GLMM)の基本と応用

臨床研究や生物統計学の分野では、患者ごとのばらつきや測定の反復性など、データが持つ複雑な構造を考慮することが不可欠である。しかし、基礎的な統計モデルでは、このような複雑性を十分に捉えきれない。そこで必要となってくるのが、「一般化線形混合モデル（Generalized Linear Mixed Model: GLMM）」である。GLMMは、線形混合モデルの柔軟性と一般化線形モデルの多様性を組み合わせることで、さまざまなタイプの応答変数と複雑なデータ構造を同時に分析できる強力なフレームワークである。本記事では、GLMMの概要から線形混合モデルとの違い、臨床研究における具体的な応用例、Rでの計算例まで、わかりやすく解説する。

一般化線形混合モデルの概要

一般化線形混合モデル（GLMM）は、固定効果と変量効果の両方を含む統計モデルであり、応答変数が正規分布に従わない場合でも適用可能である。線形混合モデルが正規分布に従う応答変数に限定されるのに対し、GLMMは、ロジスティック回帰（二項データ）、ポアソン回帰（カウントデータ）、順序ロジスティック回帰（順序カテゴリデータ）など、一般化線形モデルのフレームワークが持つ様々な分布とリンク関数を、変量効果の概念と統合したモデルである。

GLMMの主な特徴は以下の通りである。

• 非正規分布の応答変数に対応: 二項、ポアソン、ガンマ、負の二項分布など、様々な分布に対応。
• リンク関数: 応答変数の期待値と線形予測子を関連付ける関数（例：ロジット、対数、逆数）。
• 固定効果: 研究デザインによって決定される、グループ間の平均的な効果。
• 変量効果: 個体差や施設差など、ランダムなばらつきを表現する効果。これにより、同じ個体からの反復測定や、クラスター化されたデータなど、データの相関構造を適切にモデル化できる。

線形混合モデルとの違い

一般化線形混合モデルと線形混合モデルの主な違いは、応答変数の分布である。

特徴	線形混合モデル (LMM)	一般化線形混合モデル (GLMM)
応答変数の分布	正規分布に限定される	二項、ポアソン、ガンマ、負の二項分布など、多様な分布に対応
リンク関数	恒等関数（Identity link）が一般的	ロジット、対数、逆数など、応答変数の分布に応じた関数を使用
主な用途	連続データで、残差が正規分布に従うと仮定できる場合	二値、カウント、順序データなど、非正規分布に従うデータ

LMMは、例えば血圧や体重といった連続的な測定値で、これらの測定値のばらつきが正規分布に従うと仮定できる場合に適している。一方、GLMMは、疾患の有無（二値）、イベントの発生回数（カウント）、疼痛の程度（順序）など、正規分布を仮定できない多様なタイプのデータに対応できるため、より広範な臨床研究に応用可能である。

臨床研究の具体例

GLMMは、複雑な臨床データを分析するための強力なツールである。以下に具体的な例を挙げる。

二項イベントの例：新薬の投与による疾患の再発有無

ある疾患の患者を対象に、新薬とプラセボの有効性を比較する臨床試験を考える。各患者は複数回の診察を受け、それぞれの診察時に疾患の再発有無（再発あり/なし）が記録される。このデータでは、同じ患者からの再発有無の測定は独立ではない（時間的な相関がある可能性が高い）。

• 応答変数: 疾患の再発有無（二値データ：0=なし, 1=あり）
• 固定効果: 治療群（新薬 vs プラセボ）、診察時期
• 変量効果: 各患者のベースラインの再発傾向のばらつき（患者ごとの切片）

この場合、ロジットリンク関数と二項分布を仮定したGLMM（ロジスティック混合モデル）を用いることで、治療群と診察時期が再発に与える影響を評価しつつ、患者ごとのばらつきを適切に考慮した解析が可能である。

順序カテゴリの例：薬剤投与による症状改善度の評価

関節炎患者を対象に、新薬が症状改善に与える影響を評価する試験を考える。患者は一定期間、薬剤を服用し、定期的に症状の改善度を「全く改善なし」「少し改善」「かなり改善」「完全に改善」の4段階で評価する。

• 応答変数: 症状改善度（順序カテゴリデータ：1=全く改善なし, 2=少し改善, 3=かなり改善, 4=完全に改善）
• 固定効果: 治療群、投与期間
• 変量効果: 各患者の症状改善傾向のばらつき（患者ごとの切片）

この場合、累積ロジットリンク関数と多項分布（または順序ロジスティック回帰）を仮定したGLMM（順序ロジスティック混合モデル）を用いることで、薬剤が症状改善度に与える影響を評価し、患者間のばらつきを考慮した解析が可能である。

カウントデータの例：喘息患者の年間入院回数

喘息患者を対象に、新しい治療法の効果を評価する研究を考える。各患者について、治療開始後の年間入院回数がアウトカムである。入院回数は非負の整数であり、しばしば右に歪んだ（裾を引いた）分布を示す。

• 応答変数: 年間入院回数（カウントデータ）
• 固定効果: 治療法（新治療 vs 標準治療）、年齢、喫煙歴
• 変量効果: 各患者のベースラインの入院リスクのばらつき（患者ごとの切片）

この場合、対数リンク関数とポアソン分布を仮定したGLMM（ポアソン混合モデル）を用いることで、治療法が年間入院回数に与える影響を評価し、患者ごとの入院回数のばらつきを適切に考慮した解析が可能である。もしカウントデータに過分散がある場合（分散が平均より大きい場合）は、負の二項分布を仮定したGLMMを用いることも有効である。

R での計算例

R では、lme4パッケージやglmmTMBパッケージなどを用いてGLMMを簡単に実装できる。

二項イベントの計算例

ここでは、lme4パッケージを使った二項イベントの簡単な例を示す。

# 必要に応じてパッケージをインストール
# install.packages("lme4")
# install.packages("dplyr") # データ操作用
# install.packages("ggplot2") # 可視化用

library(lme4)
library(dplyr)
library(ggplot2)

# 例として、架空のデータを生成
# 20人の患者が3回ずつ診察を受け、再発の有無を記録
set.seed(123)
n_patients <- 20
n_visits <- 3

patient_id <- rep(1:n_patients, each = n_visits)
visit <- rep(1:n_visits, times = n_patients)
treatment <- sample(c("NewDrug", "Placebo"), n_patients, replace = TRUE) %>% rep(each = n_visits)

# 再発の確率を治療群と診察時期に応じて設定
# 新薬の方が再発確率が低い、診察回数が増えるごとに再発確率が変化
true_prob <- case_when(
  treatment == "NewDrug" & visit == 1 ~ 0.2,
  treatment == "NewDrug" & visit == 2 ~ 0.15,
  treatment == "NewDrug" & visit == 3 ~ 0.1,
  treatment == "Placebo" & visit == 1 ~ 0.4,
  treatment == "Placebo" & visit == 2 ~ 0.35,
  treatment == "Placebo" & visit == 3 ~ 0.3,
  TRUE ~ NA_real_ # ありえない値はNA
)

# 患者ごとのランダムな効果（切片）を追加
random_effect_patient <- rnorm(n_patients, mean = 0, sd = 0.5) %>% rep(each = n_visits)

# ロジット変換した真の確率にランダム効果を加える
logit_prob_with_re <- log(true_prob / (1 - true_prob)) + random_effect_patient

# 逆ロジット変換して、個別の観察における確率を計算
prob_outcome <- exp(logit_prob_with_re) / (1 + exp(logit_prob_with_re))

# 結果（再発有無）を生成
recurrence <- rbinom(n_patients * n_visits, 1, prob_outcome)

# データフレームの作成
df_glmm <- data.frame(
  patient_id = factor(patient_id),
  visit = factor(visit),
  treatment = factor(treatment),
  recurrence = recurrence
)

# データの確認
head(df_glmm)
str(df_glmm)

# GLMMの実行（ロジスティック混合モデル）
# recurrence ~ treatment + visit: 固定効果
# (1 | patient_id): patient_idごとの切片を変量効果として考慮
model_glmm <- glmer(recurrence ~ treatment + visit + (1 | patient_id), 
                    data = df_glmm, 
                    family = binomial(link = "logit"))

# モデル結果の表示
summary(model_glmm)

# オッズ比の計算（解釈を容易にするため）
exp(fixef(model_glmm))

# 固定効果の95%信頼区間を計算
confint_fixed <- confint(model_glmm, method="Wald", oldNames=FALSE)[-1,]
print("固定効果の95%信頼区間:")
print(confint_fixed)

# オッズ比とその95%信頼区間
odds_ratios <- data.frame(
  OR = round(exp(fixef(model_glmm))[-1], 2),
  CI_lower = round(exp(confint_fixed[-1,1]), 2),
  CI_upper = round(exp(confint_fixed[-1,2]), 2),
  p_value = round(summary(model_glmm)$coefficients[-1,4], 3)
)
print("オッズ比とその95%信頼区間:")
print(odds_ratios)

実行結果：

> # モデル結果の表示
> summary(model_glmm)
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: recurrence ~ treatment + visit + (1 | patient_id)
   Data: df_glmm

      AIC       BIC    logLik -2*log(L)  df.resid 
     71.7      82.2     -30.8      61.7        55 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-0.8745 -0.5908 -0.4233 -0.3546  2.8200 

Random effects:
 Groups     Name        Variance Std.Dev.
 patient_id (Intercept) 0        0       
Number of obs: 60, groups:  patient_id, 20

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)
(Intercept)       -0.9354     0.5696  -1.642    0.101
treatmentPlacebo   0.6671     0.6318   1.056    0.291
visit2            -0.7841     0.7379  -1.063    0.288
visit3            -1.1381     0.7906  -1.440    0.150

Correlation of Fixed Effects:
            (Intr) trtmnP visit2
tretmntPlcb -0.553              
visit2      -0.512 -0.044       
visit3      -0.471 -0.054  0.389
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

treatmentPlacebo の estimate がプラセボ群の再発オッズが新薬群と比較してどれくらい高いかを表している。ただし、対数オッズ比になっているので、のちほど、真数に変換して、確認する。

> print("オッズ比とその95%信頼区間:")
[1] "オッズ比とその95%信頼区間:"
> print(odds_ratios)
                   OR CI_lower CI_upper p_value
treatmentPlacebo 1.95     0.56     6.72   0.291
visit2           0.46     0.11     1.94   0.288
visit3           0.32     0.07     1.51   0.150

この結果が、真数にした結果である。treatmentPlacebo の 1.95 がプラセボ群の再発オッズ比の推定値である。ただし、95％信頼区間を見ると、1より小さい（つまり、関連しない方向の）可能性もあるため、名目上であっても、統計学的有意な結果ではないときは、結果の解釈・位置づけを慎重にしないといけない。

順序カテゴリの計算例

順序カテゴリデータを一般化線形混合モデルで解析する場合、ordinalパッケージのclmm()関数が便利である。この関数は、累積ロジットリンク関数を用いた順序ロジスティック混合モデルを実行できる。

# 必要に応じてパッケージをインストール
# install.packages("ordinal")
# install.packages("dplyr") # データ操作用

library(ordinal)
library(dplyr)

# 例として、架空のデータを生成
set.seed(456)
n_patients <- 60 
n_visits <- 3

patient_id <- rep(1:n_patients, each = n_visits)
treatment <- sample(c("DrugA", "DrugB"), n_patients, replace = TRUE) %>% rep(each = n_visits)
visit <- rep(c("Baseline", "Post_Treatment_1", "Post_Treatment_2"), times = n_patients)

# 患者ごとのランダムな効果（切片）を生成
random_effect_patient <- rnorm(n_patients, mean = 0, sd = 1.2) %>% rep(each = n_visits)

# 症状改善度の真の確率的な傾向を設定
linear_predictor <- case_when(
    treatment == "DrugA" & visit == "Baseline" ~ 0.2,
    treatment == "DrugA" & visit == "Post_Treatment_1" ~ 0.8,
    treatment == "DrugA" & visit == "Post_Treatment_2" ~ 1.2,
    treatment == "DrugB" & visit == "Baseline" ~ 0.1,
    treatment == "DrugB" & visit == "Post_Treatment_1" ~ 0.4,
    treatment == "DrugB" & visit == "Post_Treatment_2" ~ 0.7,
    TRUE ~ NA_real_
)

# ランダム効果を線形予測子に加える
linear_predictor_with_re <- linear_predictor + random_effect_patient

# 累積確率から症状改善度を生成
thresholds <- c(-1.5, 0.5, 2.0) 

# 症状カテゴリを割り当てる関数
get_ordinal_category <- function(lp_val, thresholds) {
    prob_cum1 <- plogis(thresholds[1] - lp_val) # P(Y <= 1)
    prob_cum2 <- plogis(thresholds[2] - lp_val) # P(Y <= 2)
    prob_cum3 <- plogis(thresholds[3] - lp_val) # P(Y <= 3)

    u <- runif(1) # 一様乱数を生成

    if (u < prob_cum1) {
        return(1) # 全く改善なし
    } else if (u < prob_cum2) {
        return(2) # 少し改善
    } else if (u < prob_cum3) {
        return(3) # かなり改善
    } else {
        return(4) # 完全に改善
    }
}

# 各行について症状改善度を生成
symptom_improvement <- mapply(get_ordinal_category, linear_predictor_with_re, MoreArgs = list(thresholds))

# データフレームの作成
df_ordinal_glmm <- data.frame(
    patient_id = factor(patient_id),
    treatment = factor(treatment),
    visit = factor(visit, levels = c("Baseline", "Post_Treatment_1", "Post_Treatment_2")), # 順序を考慮
    symptom_improvement = ordered(symptom_improvement,
        levels = 1:4
    )
)

# データの確認
head(df_ordinal_glmm)
str(df_ordinal_glmm)

# 順序ロジスティック混合モデルの実行
cat("\n=== 順序ロジスティック混合モデルの実行 ===\n")
model_ordinal_glmm <- clmm(symptom_improvement ~ treatment + visit + (1 | patient_id),
    data = df_ordinal_glmm)

# モデル結果の表示
cat("\n=== モデル結果 ===\n")
print(summary(model_ordinal_glmm))

# オッズ比の計算（解釈を容易にするため）
exp(coef(model_ordinal_glmm))

# 固定効果の95%信頼区間を計算
confint_fixed <- confint(model_ordinal_glmm, method = "Wald", oldNames = FALSE)[-1, ]
print("固定効果の95%信頼区間:")
print(exp(confint_fixed))

# オッズ比とその95%信頼区間
odds_ratios <- data.frame(
    OR = round(exp(coef(model_ordinal_glmm))[c(-1,-2,-3)], 2),
    CI_lower = round(exp(confint_fixed[c(-1,-2), 1]), 2),
    CI_upper = round(exp(confint_fixed[c(-1,-2), 2]), 2),
    p_value = round(summary(model_ordinal_glmm)$coefficients[c(-1,-2,-3), 4], 3)
)
print("オッズ比とその95%信頼区間:")
print(odds_ratios)

実行結果：

=== モデル結果 ===
> print(summary(model_ordinal_glmm))
Cumulative Link Mixed Model fitted with the Laplace approximation

formula: symptom_improvement ~ treatment + visit + (1 | patient_id)
data:    df_ordinal_glmm

 link  threshold nobs logLik  AIC    niter     max.grad cond.H 
 logit flexible  180  -241.99 497.98 342(1028) 3.59e-05 3.1e+01

Random effects:
 Groups     Name        Variance Std.Dev.
 patient_id (Intercept) 0.5177   0.7195  
Number of groups:  patient_id 60 

Coefficients:
                      Estimate Std. Error z value Pr(>|z|)
treatmentDrugB        -0.02135    0.33177  -0.064    0.949
visitPost_Treatment_1  0.40715    0.33501   1.215    0.224
visitPost_Treatment_2  0.40299    0.33866   1.190    0.234

Threshold coefficients:
    Estimate Std. Error z value
1|2  -1.4436     0.3464  -4.167
2|3   0.3255     0.3174   1.026
3|4   1.6851     0.3532   4.771

treatmentDrugA に対する DrugB の対数オッズ比や visitBaseline に対する Post_Treatment1, Post_Treatment2 のオッズ比が計算されている。Threshold も計算されている。Threshold は、最初の想定と概ね一致しているのがわかる。

[1] "オッズ比とその95%信頼区間:"
> print(odds_ratios)
                        OR CI_lower CI_upper p_value
treatmentDrugB        0.98     0.51     1.88   0.949
visitPost_Treatment_1 1.50     0.78     2.90   0.224
visitPost_Treatment_2 1.50     0.77     2.91   0.234

最後に、Treatment と Visit に関するとオッズ比と95％信頼区間が出力される。DrugB のオッズ比はほぼ1であり、Drug A と B は異なるとは言えなかったという結果になっている。

カウントデータの計算例

カウントデータ（イベント発生回数など、非負の整数値）を一般化線形混合モデルで解析する場合、応答変数がポアソン分布に従うと仮定するポアソン混合モデルが一般的である。ただし、カウントデータには過分散（overdispersion）と呼ばれる現象がしばしば見られる。これは、データの分散が平均よりはるかに大きい場合に起こり、ポアソン分布の仮定（平均と分散が等しい）が破綻していることを示唆している。このような場合、負の二項混合モデルを使用することで、過分散を適切にモデル化できる。

Rでは、lme4パッケージのglmer()関数がポアソン混合モデルを、glmmTMBパッケージのglmmTMB()関数がポアソン混合モデルと負の二項混合モデルの両方をサポートしており、より柔軟なモデル構築が可能である。ここでは、glmmTMBパッケージを使用した例を示す。

# 必要に応じてパッケージをインストール
# install.packages("glmmTMB")
# install.packages("dplyr") # データ操作用
# install.packages("ggplot2") # 可視化用

library(glmmTMB)
library(dplyr)
library(ggplot2)

# 例として、架空のデータを生成
# 50人の喘息患者が2つの治療群に割り付けられ、
# ベースライン、1年目、2年目の3時点で年間入院回数を記録
set.seed(1234) # 再現性のためのシード設定
n_patients <- 50
n_time_points <- 3 # ベースライン(0), 1年目(1), 2年目(2)

patient_id_vec <- rep(1:n_patients, each = n_time_points)
time_vec <- rep(0:(n_time_points - 1), times = n_patients) # 時点: 0, 1, 2

# 患者に割り当てられる固定効果は時間を通じて不変
treatment_vec <- sample(c("NewTreatment", "StandardCare"), n_patients, replace = TRUE) %>% 
  rep(each = n_time_points)
age_vec <- round(rnorm(n_patients, mean = 45, sd = 10)) %>% 
  rep(each = n_time_points)
smoking_history_vec <- sample(c("Non-Smoker", "Ex-Smoker", "Current-Smoker"), n_patients, replace = TRUE, prob = c(0.5, 0.3, 0.2)) %>% 
  rep(each = n_time_points)

# 患者ごとのランダムな効果（切片）を生成
random_effect_patient <- rnorm(n_patients, mean = 0, sd = 0.7) %>% # 標準偏差を少し大きくして患者差を強調
  rep(each = n_time_points)

# 入院回数の真の平均（ポアソン分布のλ）を設定
# 対数リンク関数を使用するため、線形予測子にランダム効果を加える
# 新治療の効果は1年目、2年目により顕著に出るようにする
# 時間経過による入院回数の変化 (例: 徐々に減少)
log_lambda <- 
  -1.0 +                                      # ベースラインの対数平均
  # 固定効果の影響
  ifelse(treatment_vec == "NewTreatment" & time_vec > 0, -0.6, 0) + # 新治療の効果 (治療後のみ)
  0.02 * (age_vec - mean(age_vec)) +          # 年齢の効果 (中心化)
  ifelse(smoking_history_vec == "Ex-Smoker", 0.4, 0) + # 喫煙歴の効果
  ifelse(smoking_history_vec == "Current-Smoker", 0.8, 0) +
  # 時間の効果 (例: 時間が経つと少しずつ入院が減る)
  ifelse(time_vec == 1, -0.2, 0) +            # 1年目の効果
  ifelse(time_vec == 2, -0.4, 0) +            # 2年目の効果
  # 患者ごとのランダム効果
  random_effect_patient

# 平均入院回数
lambda <- exp(log_lambda)

# 入院回数（ポアソン分布に従うと仮定）を生成
# 過分散をシミュレートするため、負の二項分布で生成してみる (ポアソンでやりたければ rpois に変更)
annual_hospitalizations <- rnbinom(n_patients * n_time_points, size = 1.5, mu = lambda) # sizeを1.5に設定して少し過分散を持たせる
# annual_hospitalizations <- rpois(n_patients * n_time_points, lambda)

# データフレームの作成
df_repeated_count <- data.frame(
  patient_id = factor(patient_id_vec),
  time = factor(time_vec, levels = 0:2, labels = c("Baseline", "Year1", "Year2")), # 時点を因子として指定
  treatment = factor(treatment_vec),
  age = age_vec,
  smoking_history = factor(smoking_history_vec, levels = c("Non-Smoker", "Ex-Smoker", "Current-Smoker")),
  annual_hospitalizations = annual_hospitalizations
)

# データの確認
head(df_repeated_count)
str(df_repeated_count)
summary(df_repeated_count$annual_hospitalizations)
table(df_repeated_count$time, df_repeated_count$treatment)

# ポアソン混合モデルの実行
# timeを固定効果として含める
# (1 | patient_id): patient_idごとの切片を変量効果として考慮
model_poisson_repeated <- glmmTMB(annual_hospitalizations ~ treatment + age + smoking_history + time + (1 | patient_id),
                                  data = df_repeated_count,
                                  family = poisson(link = "log"))

# モデル結果の表示
summary(model_poisson_repeated)

# 残差の確認（過分散の兆候がないか）
# モデル診断プロットなどを行うことが推奨される
# plot(model_poisson_repeated) # glmmTMBのplotメソッドはより詳細な診断をサポート

# もし過分散が疑われる場合は、負の二項混合モデルを使用
model_negbin_repeated <- glmmTMB(annual_hospitalizations ~ treatment + age + smoking_history + time + (1 | patient_id),
                                 data = df_repeated_count,
                                 family = nbinom2(link = "log"))

summary(model_negbin_repeated)

# モデルの比較 (AIC値)
AIC(model_poisson_repeated, model_negbin_repeated)

実行結果：

> # モデル結果の表示
> summary(model_poisson_repeated)
 Family: poisson  ( log )
Formula:          
annual_hospitalizations ~ treatment + age + smoking_history +  
    time + (1 | patient_id)
Data: df_repeated_count

      AIC       BIC    logLik -2*log(L)  df.resid 
    299.8     323.9    -141.9     283.8       142 

Random effects:

Conditional model:
 Groups     Name        Variance Std.Dev.
 patient_id (Intercept) 0.3296   0.5741  
Number of obs: 150, groups:  patient_id, 50

Conditional model:
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                   -2.76770    0.70730  -3.913 9.11e-05 ***
treatmentStandardCare          0.56851    0.34141   1.665   0.0959 .  
age                            0.02433    0.01330   1.829   0.0673 .  
smoking_historyEx-Smoker       0.58910    0.36131   1.630   0.1030    
smoking_historyCurrent-Smoker  0.85493    0.38553   2.218   0.0266 *  
timeYear1                     -0.17435    0.29601  -0.589   0.5558    
timeYear2                      0.07696    0.27756   0.277   0.7816    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

treatmentStadardCare の 0.56851 は、新治療群と比較して、標準治療群は、年間入院回数が、対数平均で 0.56851 回、真数にすると exp(0.56851) = 1.765734、つまり、約 1.8 回多いと推定された。ただし、統計学的有意ではなかった。

> # もし過分散が疑われる場合は、負の二項混合モデルを使用
> model_negbin_repeated <- glmmTMB(annual_hospitalizations ~ treatment + age + smoking_history + time + (1 | patient_id),
+                                  data = df_repeated_count,
+                                  family = nbinom2(link = "log"))
> summary(model_negbin_repeated)
 Family: nbinom2  ( log )
Formula:          
annual_hospitalizations ~ treatment + age + smoking_history +  
    time + (1 | patient_id)
Data: df_repeated_count

      AIC       BIC    logLik -2*log(L)  df.resid 
    286.6     313.7    -134.3     268.6       141 

Random effects:

Conditional model:
 Groups     Name        Variance  Std.Dev. 
 patient_id (Intercept) 1.754e-08 0.0001324
Number of obs: 150, groups:  patient_id, 50

Dispersion parameter for nbinom2 family (): 0.728 

Conditional model:
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                   -2.51040    0.72986  -3.440 0.000583 ***
treatmentStandardCare          0.56779    0.35536   1.598 0.110085    
age                            0.02244    0.01398   1.604 0.108622    
smoking_historyEx-Smoker       0.61934    0.37778   1.639 0.101125    
smoking_historyCurrent-Smoker  0.84247    0.39962   2.108 0.035014 *  
timeYear1                     -0.23412    0.39391  -0.594 0.552284    
timeYear2                      0.08020    0.37520   0.214 0.830750    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

負の二項分布を想定して計算すると、対数平均で、0.56779 で、真数にすると、exp(0.56779) = 1.764363 とこちらも、約 1.8 回、標準治療群のほうが多いという計算結果であった。しかし、こちらでも統計学的有意ではなかった。

> # モデルの比較 (AIC値)
> AIC(model_poisson_repeated, model_negbin_repeated)
                       df      AIC
model_poisson_repeated  8 299.8395
model_negbin_repeated   9 286.5831

ポアソン分布を仮定したモデルと負の二項分布を仮定したモデルを当てはまりの観点から AIC で比較すると、負の二項分布を仮定したモデルのほうが、若干 AIC が低く、比較的当てはまりはよいとの結果になった。ただし、上述の通り、どちらのモデルでも、新薬治療群と標準治療群は、統計学的有意な差はなく、両群に差があるとは言えなかった。

まとめ

一般化線形混合モデル（GLMM）は、線形混合モデルの枠組みを変量効果に拡張しつつ、一般化線形モデルの柔軟性を利用して、様々な分布を持つ応答変数に対応できる強力な統計モデリングツールである。二項イベント、順序カテゴリ、カウントデータなど、多様な臨床研究データにおいて、個体間のばらつきや反復測定の相関構造を適切に考慮した解析を可能にする。GLMMを理解し、適切に活用することで、より正確でロバストな研究結果を得ることができ、臨床的意思決定に貢献することが期待される。Rなどの統計ソフトウェアを用いることで、GLMMは容易に実行できるので、結果の解釈が適切にできるようにしていきたい。