【学習内容】
マッチド・ケースコントロール研究(疫学), パネルデータの固定効果モデル(計量経済学)のような問題設定で, アウトカムが二値の場合, 通常のロジスティック回帰では, マッチングによって生じる層（ペア）ごとのnuisance parameter（局外母数）が存在するため, 漸近性の条件が満足されず, 推定量にバイアスが生じることが知られている. このあたりがどうなっているのかを学んでいきたい.
この記事では, まず, サンプルデータのinfertデータセットで, 通常のロジスティック回帰と条件付きロジスティック回帰を比較する. その後に, シミュレーションデータを生成して, そのデータセットで比較を行う.
使用するパッケージの読み込み
library(tidyverse)
library(survival)
library(broom)
infertデータセット
infertデータセットは, 過去の流産や人工妊娠中絶が, その後の不妊に影響するかどうかを調査したマッチド・ケースコントロール研究のデータ.
1人のケース（不妊患者）に対して, 年齢（age）, 教育歴（education）, 出産回数（parity）をマッチングさせたコントロールを原則2人選出（※74番の層のみコントロールは1人）.
変数の定義
-
age: 年齢 (ケースの年齢にマッチ)
-
parity: 出産回数
-
induced: 過去の人工妊娠中絶の回数 (0 = 0回, 1 = 1回, 2 = 2回以上)
-
case: ケース/コントロールのステータス (1 = ケース, 0 = コントロール)
-
spontaneous: 過去の自然流産の回数 (0 = 0回, 1 = 1回, 2 = 2回以上)
-
stratum: 層（マッチングされたセット）の番号 (1-63)
-
pooled.stratum: プールされた層の番号 (1-83)
データの確認
head(infert)
## education age parity induced case spontaneous stratum pooled.stratum
## 1 0-5yrs 26 6 1 1 2 1 3
## 2 0-5yrs 42 1 1 1 0 2 1
## 3 0-5yrs 39 6 2 1 0 3 4
## 4 0-5yrs 34 4 2 1 0 4 2
## 5 6-11yrs 35 3 1 1 1 5 32
## 6 6-11yrs 36 4 2 1 1 6 36
summary(infert)
## education age parity induced
## 0-5yrs : 12 Min. :21.00 Min. :1.000 Min. :0.0000
## 6-11yrs:120 1st Qu.:28.00 1st Qu.:1.000 1st Qu.:0.0000
## 12+ yrs:116 Median :31.00 Median :2.000 Median :0.0000
## Mean :31.50 Mean :2.093 Mean :0.5726
## 3rd Qu.:35.25 3rd Qu.:3.000 3rd Qu.:1.0000
## Max. :44.00 Max. :6.000 Max. :2.0000
## case spontaneous stratum pooled.stratum
## Min. :0.0000 Min. :0.0000 Min. : 1.00 Min. : 1.00
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:21.00 1st Qu.:19.00
## Median :0.0000 Median :0.0000 Median :42.00 Median :36.00
## Mean :0.3347 Mean :0.5766 Mean :41.87 Mean :33.58
## 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:62.25 3rd Qu.:48.25
## Max. :1.0000 Max. :2.0000 Max. :83.00 Max. :63.00
table(infert$parity)
##
## 1 2 3 4 5 6
## 99 81 36 18 6 8
table(infert$induced)
##
## 0 1 2
## 143 68 37
table(infert$case)
##
## 0 1
## 165 83
table(infert$spontaneous)
##
## 0 1 2
## 141 71 36
自然流産（spontaneous）と人工妊娠中絶（induced）が不妊（case）に与える影響を, 3つの異なるモデルで比較.
Model 1: 通常のロジスティック回帰（未調整） マッチングの構造（層）を完全に無視して解析するモデル.
model1 <- glm(case ~ spontaneous + induced, data = infert, family = binomial())
#summary(model1)
Model 2: 層(ペアID)をダミー変数として投入した通常のロジスティック回帰
model2 <- glm(case ~ spontaneous + induced + factor(stratum),
data = infert, family = binomial())
#summary(model2)
Model 3: 条件付きロジスティック回帰 マッチングの構造をモデル化. survival パッケージの clogit() 関数を使用し, strata() で層の変数を指定.
model3 <- clogit(case ~ spontaneous + induced + strata(stratum), data = infert)
#summary(model3)
spontaneous と induced の推定値（対数オッズ比）と標準誤差を比較.
res <- bind_rows(
tidy(model1) %>% filter(term %in% c("spontaneous","induced")) %>% mutate(model = "1. Unadjusted"),
tidy(model2) %>% filter(term %in% c("spontaneous","induced")) %>% mutate(model = "2. Dummy Adjusted"),
tidy(model3) %>% filter(term %in% c("spontaneous","induced")) %>% mutate(model = "3. Conditional")
) %>%
select(model, term, estimate, std.error, p.value) %>%
mutate(
estimate = round(estimate, 2),
std.error = round(std.error, 2),
p.value = round(p.value, 4)
)
res
## # A tibble: 6 × 5
## model term estimate std.error p.value
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 1. Unadjusted spontaneous 1.2 0.21 0
## 2 1. Unadjusted induced 0.42 0.21 0.042
## 3 2. Dummy Adjusted spontaneous 3.23 0.46 0
## 4 2. Dummy Adjusted induced 2.19 0.46 0
## 5 3. Conditional spontaneous 1.99 0.35 0
## 6 3. Conditional induced 1.41 0.36 0.0001
3つ目の条件付きロジスティック回帰の結果と他の結果は異なっている.
シミュレーションデータ
サンプルデータだと, 上手く推定できているのか分からないので, シミュレーションデータを生成して動きを確認する.
ここでは, 1対1マッチングで500ペア作って, 上と同様, 3つのモデルをあてはめる. 曝露Xの係数パラメータの真の値は1としている.
set.seed(2026)
n_pairs <- 500 # 1対1マッチングのペア数（層の数）
true_beta <- 1.0 # 曝露Xの真の対数オッズ比
# データを格納する空のデータフレーム
sim_data <- tibble()
# 1. マッチド・ケースコントロールデータの生成
for (i in 1:n_pairs) {
# 層(ペア)ごとに固有のベースラインリスク(alpha_i)を設定
alpha_i <- rnorm(1, mean = 0, sd = 1.5)
valid_pair <- FALSE
while(!valid_pair) {
# 曝露Xの生成: 層の特性(alpha_i)と相関
X <- alpha_i * 0.5 + rnorm(2, mean = 0, sd = 1)
# 真のロジスティックモデルに従ってアウトカムYを生成
logit_p <- alpha_i + true_beta * X
p <- exp(logit_p) / (1 + exp(logit_p))
Y <- rbinom(2, size = 1, prob = p)
# ここでアウトカムを2人分生成している
# 1人がケース(Y=1), 1人がコントロール(Y=0)になったらペアとして採用
if (sum(Y) == 1) {
sim_data <- bind_rows(sim_data, tibble(
stratum = i,
alpha = alpha_i, # 観測されない背景因子
X = X,
Y = Y
))
valid_pair <- TRUE
}
}
}
# 2. 3つのモデルによる推定の比較
# Model 1: 通常のロジスティック回帰（未調整）
model1 <- glm(Y ~ X, data = sim_data, family = binomial())
# Model 2: 層(ペアID)をダミー変数として投入した通常のロジスティック回帰
model2 <- glm(Y ~ X + factor(stratum), data = sim_data, family = binomial())
# Model 3: 条件付きロジスティック回帰
model3 <- clogit(Y ~ X + strata(stratum), data = sim_data)
# 3. 結果の抽出と表示 (曝露Xの係数のみを比較)
res <- bind_rows(
tidy(model1) %>% filter(term == "X") %>% mutate(model = "1. Unadjusted"),
tidy(model2) %>% filter(term == "X") %>% mutate(model = "2. Dummy Adjusted"),
tidy(model3) %>% filter(term == "X") %>% mutate(model = "3. Conditional")
) %>%
select(model, estimate, std.error) %>%
mutate(
estimate = round(estimate, 2),
std.error = round(std.error, 2),
true_beta = true_beta)
res
## # A tibble: 3 × 4
## model estimate std.error true_beta
## <chr> <dbl> <dbl> <dbl>
## 1 1. Unadjusted 0.81 0.07 1
## 2 2. Dummy Adjusted 1.97 0.13 1
## 3 3. Conditional 0.98 0.09 1
ダミー変数で調整したものは, 真の値の2倍近い値になってしまっている. これはIncidence Parameter Problemで, 2倍になることは, Chamberlain (1980)のSection 2に載っている. 数理を理解したわけではないので, このあたりは次回以降学びたい.
1:Mのマッチングだと(M+1)/M倍となる?(パネルデータでT時点のデータなら, T/(T-1)倍になる?). このあたりについても次回以降で確認したい.
参考文献(関連していると思われるからメモしているだけのものも含む)
-
藤井良宜『カテゴリカルデータ解析 (Rで学ぶデータサイエンス 1)』
https://amzn.to/4dfxJ4U
-
丹後俊郎・松井茂之『医学統計学ハンドブック』
https://amzn.to/41jOPrf
-
Neyman and Scott (1948) “Consistent Estimates Based on Partially Consistent Observations”
-
Andersen (1970) “Asymptotic Properties of Conditional Maximum-Likelihood Estimators”
-
Chamberlain (1980) “Analysis of covariance with qualitative data.”
-
Breslow (1981) “Odds ratio estimators when the data are sparse.”
-
Lancaster (2000) “The incidental parameter problem since 1948”
