「受動喫煙防止法について論点整理①：受動喫煙による健康リスク・死亡者数の推定はどのくらい信用できるか？」と言うブログのエントリーで、受動喫煙の相対リスク1.22倍が未測定交絡因子によって引き起こされたバイアスだとして、未測定交絡因子の相対リスクは何倍ぐらい無いといけないかと言う議論をしていて、ノンパラメトリックな分析方法を使って最低でも1.74としているのですが、参照している分析手法の論文の記述をよく見ると、もっと低い値でも済みそうなので*1数値例を作ってみました。
set.seed(5301929)
### データ生成 ###
# 観測数
n <- 300# 観察不能な1か0の未測定交絡因子
X <- c(rep(0, n/2), rep(1, n/2))# 90%の確率で未測定交絡因子と一致する観察可能因子
CV <- 0.9
E <- c(ifelse(runif(n/2) < CV, 0, 1), ifelse(runif(n/2) < CV, 1, 0))# ロジット関数
plogit <- function(p){
exp(p)/(1 + exp(p))
}# ロジットモデルにおける定数項と未測定交絡因子Xの係数
CONS <- -1
COEF_X <- 0.5# イベント発生の有無を示す1か0かの被説明変数
D <- 0 + (runif(n) < plogit(CONS + COEF_X*X))*1sprintf("未測定交絡因子Xの相対リスク: %.4f", plogit(CONS + COEF_X)/plogit(CONS + 0))
sprintf("未測定交絡因子Xと観察可能因子Eの相関係数: %.4f", cor(X, E))
# 観察可能因子と未測定交絡因子のリスク比RR_EXを計算する（注意：論文ではRR_{EU}）
r_xtabs <- xtabs(~X + E, data=data.frame(D = D, X = X, E = E))
RR_EX <- max(r_xtabs[4]/r_xtabs[3], r_xtabs[2]/r_xtabs[1]) # = max(P(X=1|E=1)/P(X=1|E=0), P(X=0|E=1)/P(X=0|E=0))
sprintf("観察可能因子Eと未測定交絡因子Xのリスク比: %.4f", RR_EX)
### 説明変数をEで、ロジット関数の推定をする ###
r_E <- glm(D ~ E, family = binomial("logit"))
sprintf("観察可能因子Eの観察された相対リスク: %.4f", plogit(coef(r_E)[1] + coef(r_E)[2])/plogit(coef(r_E)[1]))
（追記：nを30000、CONSをlog(0.0038/(1-0.0038))、COEF_Xを0.32ぐらいにすると、実際の分析データに近い設定になります。）
未測定交絡因子の相対リスクは1.40倍で、観察可能因子Eの観察された相対リスクは1.35倍になりました。未測定交絡因子と観察可能因子の強い相関を仮定したら必要となる相対リスクは低いものの、そうでないと高くなると言う話のようです。
*1: the observed relative risk of cigarette smoking on lung cancer was ... The work of Cornfield et al. showed that for a binary unmeasured confounder to completely explain away the observed relative risk, both the exposure-confounder relative risk and the confounder-outcome relative risk would have to be at least 10.73
