はじめに

数年前に岩波書店の確率と情報の科学シリーズの星野崇宏著「調査観測データの統計科学」を読んでいきまとめていました。

調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)

作者:星野崇宏
岩波書店

Amazon

しかし、感度分析の部分について放置していたので、今回まとめました。

隠れた共変量によるバイアス

適切に因果効果を推定するにあたり、共変量は全て観測されている状態であることが求められます。しかし、実際の観測データでは全ての共変量は観測されていない可能性が高いです。

観測されていない共変量が存在するとき、重回帰分析や傾向スコアを用いた分析により推定した因果効果にバイアスがのってしまいます。

未観測の共変量の影響を分析する

隠れた/観測されてない共変量が存在している場合に因果効果の推定値にバイアスが乗る事は避けされないことがわかっています。そこで、未観測の共変量をモデリングし、その影響力を変化させることで、得られる因果効果の変動を調べることができます。この分析方法は感度分析(Sensitivity analysis)と呼ばれています。

この感度分析では、処置変数の割当と結果変数を説明する二つのモデルを考えていきます

一つ目は、未観測の共変量 $u$ を含んだ処置変数を説明するロジスティック回帰モデルであり、次のように表現します。

$\displaystyle{ P(W=1|\textbf{X}, U) = \frac{exp (\gamma^t \textbf{X} + \alpha U )}{1+exp (\gamma^t \textbf{X} + \alpha U )} }$

二つ目は結果変数を説明する回帰モデルで、次のように表現します。

$\displaystyle{ Y \sim N(\tau w +\beta^t \textbf{X}+\delta U+\epsilon_i, \sigma^2) }$

(本来は $\beta$ はベクトルなのでボールド体です)

感度分析の基本的な方法は、ロジスティック回帰モデルにおける隠れた共変量 $U_{i} (0 \le U \le 1 )$ とその係数 $\gamma$ を変化させることで、隠れた共変量の影響を調べる方法です。

しかし、この方法では「隠れた共変量が結果変量に与える影響」を考慮できないため、隠れた共変量が因果効果の推定に影響を及ぼす効果を過小評価している可能性があります。

そのため、割り当てを説明するモデルと結果変数を説明するモデルの両方の隠れた共変量の影響を分析する方法が提案されている。

具体的には、 $\gamma$ と $\delta$ を固定すれば、他の母数は最尤推定することができることから、 $\gamma$ と $\delta$ を様々変化させ両方の回帰モデルを推定を行うことで、「隠れた共変量の割り当てや結果変数の影響力」、「因果効果の推定値への影響力」を評価することができるようになります。

$\gamma$ と $\delta$ を固定して介入効果を推定するには、以下の尤度関数を用います。

$\displaystyle{ \begin{align} L(\tau, \beta, \sigma^2 , \gamma,\alpha, \delta) = \sum_{i=1} ^n \rm{ln} \left[ \frac{1}{2} \left( \frac{1}{\sqrt{2 \pi \sigma^2}} \right) \times exp \left( -\frac{1}{2 \sigma^2} (Y_i -\tau W_i-\beta' \textbf{X}_1)^2\right) \\ \times \frac{(exp(\gamma ' \textbf{X}_i))^{W_i}}{1-exp(\gamma ' \textbf{X}_i)} + \frac{1}{2} \left( \frac{1}{\sqrt{2 \pi \sigma^2}} \right) \\ \times exp \left( \frac{1}{2 \sigma^2} (Y_i -\tau W_i- \beta ' \textbf{X}_i )^2 \right) \\ \times \frac{(exp(\gamma ' \textbf{X}_i + \alpha ))^{W_i}}{1-exp(\gamma ' \textbf{X}_i + \alpha )} \right] \end{align} }$

Imbens(2003)では、この方法による感度分析だと結果が解釈しにくいとしており、未観測共変量から結果変数と処置変数への影響を偏決定係数(partial R²)を用いた感度分析の方法が用いられています。

partial R²(偏決定係数)とは、多変量で目的変数を表現する場合にある一つの説明変数が目的変数の変動量をどれほど表現しているかを表した指標です。

まず、partial R²を使って未観測交絡変数から結果変量への影響を表現する。

未観測共変量から結果変数へ影響の強さを表現する偏決定係数 $R^2$ は次のように表現される。

$\displaystyle{ R^2_{Y,par}(\alpha, \delta) = \frac{R^2(\alpha, \delta) -R^2(0,0)}{1-R^2_Y(0,0)} }$

ここで、 $R^2 _Y ( \alpha , \delta ) = 1 − \sigma ^2 ( \alpha , \delta ) / \sum Y$ と $\sum _Y = \sum _i ( Y _i − \bar{Y} ) ^2/N.$ である。

未観測共変量から処置変数へ影響の強さを表現するモデルはロジスティック回帰であり、上記の係数は利用できない。そのため、ここではロジスティック回帰のための次のような決定係数 $R^2$ を考える。

$\displaystyle{ R^2_{W}(\alpha, \delta) = \frac{\hat{\gamma}(\alpha, \delta)' \Sigma_X \hat{\gamma}(\alpha, \delta) + \alpha^2/4}{\hat{\gamma}(\alpha, \delta)'\Sigma_X \hat{\gamma}(\alpha, \delta) + \alpha^2/4 +\pi^2/3} }$

(本来は $\Sigma_X$ はベクトルなのでボールド体) ここから、ロジスティック回帰モデルにおける偏決定係数は次のように定義される。

$\displaystyle{ R^2_{W,par}(\alpha, \delta) = \frac{R^2_W(\alpha, \delta) -R^2_W(0,0)}{1-R^2_W(0,0)} }$

これら $\hat \tau , R ^2 _{Y,par}, R^2 _{W,par}$ のペアを $\alpha, \delta$ を変化させて取得していき、同時に変化する因果効果 $\tau$ の関係性を見ていきます。具体的には、どれほど処置変数と結果変数への影響が強い未観測の変数があれば、現在観測されている因果効果を小さくできるのかや0にすることができるのかを議論します。

例えば、因果効果を0にするためには $\alpha=0.5, \delta=-0.8$ に対応するペアの影響が未観測変数から必要であるとなった場面を考えます。

現状の観測済みの共変量と比較し、ある共変量の値が $\alpha=0.5, \delta=-0.8$ に対応するペアに一致する場合、その共変量と同じような未観測共変量があれば、結果は覆ることになります。

共変量と同じような未観測共変量があるかどうかは、ドメイン知識や過去の取り組みから判断することになります。

また、 $\Gamma=exp(\gamma)$ 持つとすると、上記のロジスティック回帰において、全く等しい共変量の値を持つ対象 $i$ と $i'$ のオッズ比は次のようになる。

$\displaystyle{ \frac{1}{\Gamma} \le \frac{{p(z_i=1|x_i)}{(1-p(z_i'=1|x_i'))}}{{p(z_i'=1|x_i')}{(1-p(z_i=1|x_i))}} \le \Gamma }$

この性質を使用して二項検定やWilcoxonの符号順位検定等の検定を用いた感度分析方法も提案されている。(なお、感度分析はこれ以外にも様々な方法が提案されています。)

ここら辺の理解については、こちらのローゼンバウムの訳書が役に立つかもしれません。

ローゼンバウム統計的因果推論入門: 観察研究とランダム化実験

作者:Paul R. Rosenbaum
共立出版

Amazon

Rで感度分析

Rで感度分析を行なってみます。

用いたデータはMatchingパッケージに含まれるlalondeです。このデータは1976年に実施された米国職業訓練プログラムを受講した群としていない群で、その後の年収を追跡したものです。今回は1978年の年収への職業訓練プログラム受講の効果を推定する場面を考え、2値（バイナリ）の未観測共変量があることを仮定した感度分析を行います。

今回の感度分析では、回帰係数による方法と偏決定係数による方法の両方を実行してみました。観測済み共変量には年齢、教育年数等の7個を用います。

なお、こちらの資料に記載されているコードを~~パクリ~~ベースに計算を行なっています。

xmatt <- lalonde %>%
  dplyr::select("age", "educ", "black", "hisp", "married", "nodegr", "u74") %>% 
  mutate(const = 1) %>% 
  relocate(const, 1)
W <- lalonde$treat
Y <- lalonde$re78

# 処置モデルと介入モデルを作成する
lm.model <- lm(Y ~ .-1, data = cbind(xmatt, W))
glm.model <- glm(W ~ .-1, data = xmatt, family=binomial(link = "logit"))


# 未観測変数がない場合の結果を取得
Izero <- imbens(0, 
                0, 
                as.matrix(xmatt), 
                stvalnt)

S2zero <- Izero$Sigma2
R2wzero <- (t(Izero$Gamma) %*% Izero$VC %*% Izero$Gamma) /
           (t(Izero$Gamma) %*% Izero$VC %*% Izero$Gamma + ((pi^2)/3))

次に最尤推定によるパラメータ推定を行う関数の定義です。

imbens <- function(alpha, delta, xmat, stvaln) {
  logl <- function(b, alpha, delta, X, W, Y) {
    tau <- b[length(b)]
    lns2 <- b[length(b)-1]
    s2 <- exp(lns2)
    gamma <- b[1:ncol(X)]
    beta <- b[(ncol(X)+1):(length(b)-2)]
    llik <- log(
      0.5*(1/sqrt(2*pi*s2)
           *exp(-(1/(2*s2))*(Y-tau*W-X %*% beta)^2)
           *((exp(X%*%gamma))^W)/(1+exp(X %*% gamma)))
      + 0.5*((1/sqrt(2*pi*s2))
             *exp(-(1/(2*s2))*(Y-tau*W-X %*% beta-delta)^2)
             *((exp(X %*% gamma+alpha))^W)/(1+exp(X %*% gamma+alpha)))
    )
    sum(llik)
  }
  imbens.mle <- optim(stvaln, logl, hessian = F, method = "BFGS",
                      control = list(fnscale = -1, trace = 1, maxit = 2500, reltol = 1e-17),
                      alpha = alpha, delta = delta, X = xmat, W = W, Y = Y )
  sigma2 <- exp(imbens.mle$par[length(stvaln)-1])
  tau <- imbens.mle$par[length(stvaln)]
  gamma <- imbens.mle$par[2:ncol(xmat)]
  vc <- as.matrix(var(xmat[,-1]))
  stvaln <- imbens.mle$par
  result <- list( Sigma2 = sigma2, Tau = tau, Gamma = gamma, VC = vc, mle=imbens.mle)
  result
}

回帰係数を最尤推定で求める方法

# 結果変量と処置変数への回帰係数をランダムに生成
N <- 3000
alpha <- rnorm(N, 0, 1) 
delta <- rnorm(N, 0, 1) 

## 様々なパラメータで隠れた共変量を用いて影響度を算出
mase_results_df <- data.frame()
for (i in 1:N) {
  Izero <- imbens(matrix(rep(alpha[i], length(Y)), ncol = 1), 
                  matrix(rep(delta[i], length(Y)), ncol = 1), 
                  as.matrix(xmatt), 
                  stvalnt)
  estALPHA <- alpha[i]
  estDELTA <- delta[i]
  estTAU <- Izero$Tau
  
  res <- data.frame(alpha = estALPHA,  delta = estDELTA, tau = estTAU)
  mase_results_df <- rbind(mase_results_df, res)
}

## 共変量の回帰係数を取得
coef_df <- data_frame(col = names(xmatt[-1),
                      alpha = coef(lm.model)[c(-1,-9)], 
                      delta = coef(glm.model)[-1])

## 想定する影響の大きさの定義
bais_rate <- (1671-500)/1671 
bais_rate2 <- (1671-1000)/1671
bais_rate3 <- (1671-1500)/1671

## 基準線のためのデータ生成
pair_target_df <- 
  mase_results_df %>% 
  mutate(color = ifelse(round(tau,2) == round(bais_rate2 * coef(lm.model)["W"], 2), "p","o")) %>% 
  filter(color=="p") %>% 
  rowwise() %>% 
  mutate(inv_delta = 1/delta)

$1671の結果が-$1000($671)となる場合の $\alpha, \delta$ のペアの曲線と、共変量の回帰係数を描写する。

あまり良い方法ではないですが、ペアの曲線は反比例の関数を推定して引いています。

f:id:saltcooky:20210728222943p:plain

ここで、観測されている共変量の中では黒人かどうか(black)が比較的大きな影響があります。現在推定されている職業訓練の効果を$1000低下させるには、黒人であることの影響よりも3~4倍程度影響が大きい未観測変数がないといけないということが言えそうです。

正しいドメイン知識を有していないため正確なことはわかりませんが、感覚としてはかなり大きい影響を持つ変数がないと結果は覆りにくいのではないかと思います。

偏決定係数による方法

Imbens(2003)のfigure1を再現することに対応しています。

aldelval <- cbind(runif(N, -3, 3), runif(N, -3, 3)) # alphaとdeltaをランダム生成

rsc <- map(1:N, Rsquare, xmatt, aldelval, S2zero, R2wzero) %>% 
  bind_rows()

rsc_all <- cbind(aldelval, rsc)
colnames(rsc_all) <- c("alpha","delta","Tau","R2p","R2wp","S2")

totalR2w <- R2wzero

# 共変量の偏重相関係数の算出
partials <- function(i, df){
  xmatp <- df[,-i]
  stvalnp <- setting(xmatp)
  imb <- imbens(0, 0, as.matrix(xmatp), stvalnp)
  gammas <- as.matrix(imb$Gamma)
  vxb <- var(as.matrix(xmatp[,-1]) %*% gammas)
  exclR2wp <- vxb/(vxb + pi^2/3)
  parR2w <- (totalR2w - exclR2wp)/(1 - totalR2w)
  s2p <- imb$Sigma2
  parR2 <- (s2p - S2zero)/s2p
  par <- data.frame(parR2w = parR2w, parR2=parR2)
  return(par)
}

parR2_df <- 
  map(2:length(xmatt), partials, df=xmatt) %>% 
  bind_rows() %>% 
  mutate(col = names(xmatt[,-1])) %>% 
  filter(col %in% c("age", "black", "educ")) #描画の都合で3つに絞る

こちらにおいても得られている因果効果が$1000減少する場合の曲線を引きます。

## 影響度が大きいと考えられる部分色を変える 
results <- 
  rsc_all %>% 
  mutate(color = ifelse(round(Tau,2) == round(bais_rate * coef(lm.model)["W"],2), "500", NA)) %>% 
  mutate(color = ifelse(round(Tau,2) == round(bais_rate2 * coef(lm.model)["W"],2), "1000", color)) %>% 
  mutate(color = ifelse(round(Tau,2) == round(bais_rate3 * coef(lm.model)["W"],2), "1500", color)) 

pair_target_df <- 
  results %>% 
  filter(color=="1000", R2wp>0) %>% 
  filter(R2wp < 0.25, R2p <0.25) %>% 
  rowwise() %>% 
  mutate(inv_R2wp = 1/R2wp)