はじめに

@a_bicky さんから指摘された内容について、僕が考えていたことと僕が間違えていたことと、僕が調べたことをメモとして記載する。

僕が間違えたやりとりはこちら。

僕がこのやりとりで考えていた事

複数の変数(treat, x1,x2)がy(outcome)と関連しているとする。
自分がoutcomeとの因果効果を知りたいのは複数の変数の内treatのみであるとする。
x1,x2 のうちx1 はTreatment と関連しており、x2はTreatment に対して独立している。
outcomeとtreatmentの関連を推計したいときに調整すべき変数は何か？

疫学系の人だとここでDAGを頭に浮かべるだろう。今回の状況をDAGに表すと、以下のリンク通り である。
Treat とOutcome の間にバイアスパス（赤色のライン）が存在するのがわかる。
x1 をブロックすれば(x1のノードを選択してaをタイプすればブロックした事になる)、バイアスパスは消えるためTreat のOutcome に対する効果を正しく推計出来る事になる。

よって、僕は「相関 0 の 2 つの説明変数の片方を除くいても係数は変わらない。」と考えた。

問題点

しかし、この考えは誤っていた。
outcomeが連続変数で回帰を一般線形回帰で回帰する場合は正しいが、outcomeがカテゴリカル変数でlogistic modelを用いて一般的なodds ratio を報告する場合は正しくはない。
以下に具体例とその理由、対処を示す。

シミュレーション

パラメーターセッティング
x1(treat),x2,x3のy(Outcome)に対する効果を推計する。
いずれの変数もyと関連しておりそれぞれの係数は、beta1,beta2,beta3である。
treatとx2は相関しており、相関係数はbeta1_beta2である。
treatとx3は独立している。

モデル
モデルは以下のmodel1~model4 を用いた。

一般線形回帰で係数を報告する。

連続量のアウトカムを上記計数から設定し、100回ほど乱数を生成した。
一般線形モデルで回帰した結果をストックし、1-4のモデルにおけるtreatの係数を箱ひげ図にて示した。

boot <- 100
beta1 <- 2
beta2 <- -2
beta3 <- 1.5
beta1_beta2 <- 1
res_mat <- matrix(,boot,4)
res_mat_2 <- matrix(,boot,2)
res_mat_3 <- matrix(,boot,3)
res_mat_4 <- matrix(,boot,3)

for(i in 1:boot){
  n <- 10000
  x1 <- rnorm(n,mean = 0,sd=10)
  e2 <- rnorm(n,mean = 0,sd=10)
  x2 <- beta1_beta2*x1 + e2
  x3 <- rnorm(n,mean = 0,sd=10)
  e <- rnorm(n,mean = 0,sd=10)
  y <- beta1*x1+beta2*x2+ beta3*x3 +e
  
  res <- lm(y~x1+x2+x3)
  res_mat[i,] <- res$coefficients

  res_2 <- lm(y~x1)
  res_mat_2[i,] <- res_2$coefficients

  
  res_3 <- lm(y~x1+x3)
  res_mat_3[i,] <- res_3$coefficients
  
  res_4 <- lm(y~x1+x2)
  res_mat_4[i,] <- res_4$coefficients

}
treat_coef <- cbind(res_mat[,2],res_mat_2[,2],res_mat_3[,2],res_mat_4[,2])
colnames(treat_coef) <- c("model1","model2","model3","model4")
boxplot(treat_coef, xlab = "model", ylab = "treat_coef")

ご覧の通り、モデル1の結果と合致したのはモデル4の結果であった。
つまり、調整すべき変数はx2のみであった。
JD Angrist[1]らによれば、単純な一般線形回帰の場合、変数をOmit した際のTreatmentの係数は\(\rho + \gamma'\sigma_{As}\) で表される。
ここで\(\rho\) はTreat の真の係数、\(\gamma\)を上記モデルでいうbeta2、\(\sigma_{As}\) をx2をTreatに回帰した際の回帰係数とする。
つまり、TreatとOmitted Variable が関連しない場合、係数は正しく推計される。

一般化線形回帰（ロジスティックモデル）でOdds ratioを報告する。

次にロジスティックモデルの、シミュレーション結果を示す。
用いた分布とリンク関数、そして使用モデルが一般化線形モデルでロジスティック回帰を使用している点を除き設定は上記と同様である。 縦軸にはtreatの係数をodds ratioに変換した結果(exp(coef))を示した。

boot <- 100
beta1 <- 2
beta2 <- -2
beta3 <- 1.5
beta1_beta2 <- 1
res_mat <- matrix(,boot,4)
res_mat_2 <- matrix(,boot,2)
res_mat_3 <- matrix(,boot,3)
res_mat_4 <- matrix(,boot,3)
mf_mat <- matrix(,boot,4)
mf_mat_2 <- matrix(,boot,2)
mf_mat_3 <- matrix(,boot,3)
mf_mat_4 <- matrix(,boot,3)

for(i in 1:boot){
  n <- 10000
  x1 <- rnorm(n,mean = 0,sd=10)
  e2 <- rlogis(n,0,1)
  x2 <- beta1_beta2*x1 + e2
  x3 <- rnorm(n,mean = 0,sd=10)
  e <- rlogis(n,0,1)
  z <- beta1*x1+beta2*x2+ beta3*x3 +e
  y_pred <- 1/(1+exp(-z))
  y <- ifelse(y_pred<0.5,0,1)
  
  res <- glm(y~x1+x2+x3, family = binomial(logit))
  res_mat[i,] <- res$coefficients
  mf_mat[i,] <- coef(res) * mean((dlogis(predict(res,type = "link"))))
  
  res_2 <- glm(y~x1, family = binomial(logit))
  res_mat_2[i,] <- res_2$coefficients
  mf_mat_2[i,] <- coef(res_2) * mean((dlogis(predict(res_2,type = "link"))))
  
  res_3 <- glm(y~x1+x3, family = binomial(logit))
  res_mat_3[i,] <- res_3$coefficients
  mf_mat_3[i,] <- coef(res_3) * mean((dlogis(predict(res_3,type = "link"))))

  
  res_4 <- glm(y~x1+x2, family = binomial(logit))
  res_mat_4[i,] <- res_4$coefficients
  mf_mat_4[i,] <- coef(res_4) * mean((dlogis(predict(res_4,type = "link"))))
}
treat_coef <- cbind(res_mat[,2],res_mat_2[,2],res_mat_3[,2],res_mat_4[,2])
colnames(treat_coef) <- c("model1","model2","model3","model4")
boxplot(treat_coef, xlab = "model", ylab = "treat_coef")

ご覧の通り、model1以外全てのモデルがbiasされている。

JM Wooldrige[1] によれば、ロジスティックモデルではOutcome に関連しTreatment には完全独立な変数であってもomitすることで係数を0にする方向のbiasがかかる。
具体的には係数を\(\frac{1}{\sqrt{1+\beta^2\tau^2}}\)倍にする。
ここで\(\beta\) は上記モデルのbeta3 をさし、\(\tau\) はx3 の分散である。
係数の正負(と多分p値)は変わらない。よって、outcomeに関連する変数を全て含まないとodds ratioが正しくないことになる。
これは、ロジスティックモデルを使って因果効果を推計したい場合にそのハードルをかなり高くする。

対応

では、この事態をどう回避すればよいのか？
Average Partial Effect(APE) を推計すれば回避できる。
APE とはロジスティックモデルの結果を利用した仮想的なRisk Differnce にあたる。
2値変数x のAPEを求める場合、すべての症例がx=1 であった場合の推計イベント割合から、すべての症例がx=0 であった場合の推計イベント割合を引いて算出される。
具体的には以下の式による[2]。

\(g(z) = \frac{1}{1+exp(-z)}\)

\(APE = \hat{\beta_K}(N^{-1}\sum_{i=1}^{N}{g(x_i\hat{\beta})}) \cdots x_K がcontinuous\)

\(APE = N^{-1}\sum_{i=1}^{N}{[g(\hat{\beta_1}+\hat{\beta_2}x_{i2}+\cdots+\hat{\beta_{K-1}}x_{i,K-1}+\hat{\beta_{K}})-g(\hat{\beta_1}+\hat{\beta_2}x_{i2}+\cdots+\hat{\beta_{K-1}}x_{i,K-1})]} \cdots x_K がbinary\)

先ほどの推計からAPEを推計した結果を下記に示す。

ご覧の通り、モデル1とモデル4の結果が合致している。 結果に大きく影響しTreatと独立した変数が無視されている場合、APEを報告した方が何かと無難っぽい。
ちなみに、Treat に関連した変数が無視しされた場合は内生性の問題になるので対応は別問題になる。
今まであんまり意識してなかったなぁ。あと、各種Reporting Guideもodds ratioの報告を推奨しているけどどうなんこれ？

参考資料

[1] Angrist JD, Pischke J-S. Mostly harmless econometrics : an empiricists companion. Cram101; 2013. Section 3.2.2.
[2] Wooldridge JM. Econometric analysis of cross section and panel data. 2nd ed. MIT Press; 2010. Section 15.7.