準備

いつも通り今回使うパッケージを読み込む。

library('foreign')
library('MASS')
library('dplyr')
library('ggplot2')
library('rethinking')

欠落変数バイアスと処置後変数バイアス

重回帰分析では複数の予測変数を使う。複数の予測変数を使う理由の1つは、ある結果に影響を与える原因が複数あると考えられるからである。そのようなとき、原因と考えられる複数の予測変数を回帰分析に含めるというの自然な発想である。しかし、結果変数の原因の中には、必ず回帰分析に含める必要があるものもあれば、回帰分析にいれてもいれなくてもよいものや、回帰分析にいれてはいけないものもある。回帰分析では主な予測変数以外の変数を統制変数 (control variables) と呼ぶことがあるが、回帰分析で統制 (control) すべき変数はどのようなものだろうか。


欠落変数バイアス

2つの変数 \(y\) と\(x\) があり、この2変数間の相関が高いとする。このとき、\(x\)が\(y\)の原因であるとは限らない。1つの可能性は、\(y\)が\(x\)の原因であるというものである。因果の向きが逆の場合は比較的見抜きやすいので、ここではその可能性はとりあえず考えない（実際の研究では、フィードバック効果などもあり、注意すべき問題である）。

もう1つの可能性は、第3の変数 \(z\) が存在し、この\(z\)が\(x\)の原因でもあり、\(y\)の原因でもあるという場合である。 \(x\)と\(y\)の相関が\(z\)によって説明されてしまうとき、\(x\)と\(y\)の相関は、見せかけの相関 (spurious correlation) と呼ばれる。また、実際に\(x\)が\(y\)の原因だとしても、 \(z\)のように\(x\)と\(y\)の両者に影響する変数があるかもしれない。このような\(z\) は、交絡変数 (confouding variable or confounder) と呼ばれる。

交絡変数は、必ず統制する必要がある（厳密には、バックドアパスをブロックすればよいのだが、詳しい話は統計的因果推論の教科書を参照されたい）。交絡変数を統制しないと、推定にバイアスが生じる。このバイアスを欠落変数バイアス (omitted variable bias) と呼ぶ。

シミュレーションで確認してみよう。

まず、データを作る。ここでは、x, z1, z2 という3つの変数がyの原因となっている（muを計算する行を参照）。また、z1はxの原因でもあるので、z1は交絡変数である。したがって、z1を統制しないと、xの係数が正しく推定できないはずである。z4は結果変数とは無関係の変数である。

set.seed(777)
n <- 100
z1 <- runif(n, -5, 5)
z2 <- runif(n, -5, 5) 
z3 <- runif(n, -5, 5)
x <- 0.2*z1 + rnorm(n)
mu <- 0.5 + 0.4*x + 0.2*z1 - 0.4*z2 
y <- rnorm(n, mu, sd = 1)
df <- data.frame(y, x, z1, z2, z3)

まず、正しいモデルを作り、パラメタを推定する。

true_model <- alist(
    y ~ dnorm(mu, sigma),
    mu <- alpha + beta*x + gamma1*z1 + gamma2*z2,
    alpha ~ dnorm(0, 10),
    beta ~ dnorm(0, 10),
    gamma1 ~ dnorm(0, 10),
    gamma2 ~ dnorm(0, 10),
    sigma ~ dunif(0, 50)
)
true_fit <- map(true_model, data = df)
precis(true_fit)
##         Mean StdDev  5.5% 94.5%
## alpha   0.58   0.11  0.40  0.75
## beta    0.42   0.10  0.26  0.58
## gamma1  0.21   0.04  0.14  0.27
## gamma2 -0.36   0.04 -0.42 -0.29
## sigma   1.11   0.08  0.99  1.24

xの係数であるbeta に注目すると、パラメタとして設定した 0.4 にある程度近い値が得られる。

次に、交絡変数であるz1を除外したモデルでパラメタを推定する。

omit_model1 <- alist(
    y ~ dnorm(mu, sigma),
    mu <- alpha + beta*x + gamma2*z2,
    alpha ~ dnorm(0, 10),
    beta ~ dnorm(0, 10),
    gamma2 ~ dnorm(0, 10),
    sigma ~ dunif(0, 50)
)
fit_omit1 <- map(omit_model1, data = df)
precis(fit_omit1)
##         Mean StdDev  5.5% 94.5%
## alpha   0.57   0.12  0.37  0.76
## beta    0.65   0.10  0.49  0.81
## gamma2 -0.36   0.05 -0.43 -0.29
## sigma   1.24   0.09  1.10  1.38

このモデルでは、beta の推定値が0.65 になり、かなり過大に推定されている。

続いて、\(y\) の原因ではあるが、交絡変数ではないz2 を除外してみる。

omit_model2 <- alist(
    y ~ dnorm(mu, sigma),
    mu <- alpha + beta*x + gamma1*z1,
    alpha ~ dnorm(0, 10),
    beta ~ dnorm(0, 10),
    gamma1 ~ dnorm(0, 10),
    sigma ~ dunif(0, 50)
)
fit_omit2 <- map(omit_model2, data = df)
precis(fit_omit2)
##        Mean StdDev 5.5% 94.5%
## alpha  0.67   0.15 0.43  0.91
## beta   0.34   0.13 0.12  0.56
## gamma1 0.21   0.06 0.11  0.30
## sigma  1.48   0.10 1.32  1.65

ここでは、beta が過小推定されたが、それほど大きく外れたわけではない。

最後に、yの原因ではない（関係のない）z3 を加えて回帰分析をしてみよう。

extra_model <- alist(
    y ~ dnorm(mu, sigma),
    mu <- alpha + beta*x + gamma1*z1 + gamma2*z2 + gamma3*z3,
    alpha ~ dnorm(0, 10),
    beta ~ dnorm(0, 10),
    gamma1 ~ dnorm(0, 10),
    gamma2 ~ dnorm(0, 10),
    gamma3 ~ dnorm(0, 10),
    sigma ~ dunif(0, 50)
)
fit_extra <- map(extra_model, data = df)
precis(fit_extra)
##         Mean StdDev  5.5% 94.5%
## alpha   0.57   0.11  0.40  0.75
## beta    0.42   0.10  0.26  0.59
## gamma1  0.21   0.04  0.14  0.28
## gamma2 -0.36   0.04 -0.42 -0.29
## gamma3 -0.01   0.04 -0.07  0.05
## sigma   1.11   0.08  0.99  1.24

betaについてはほぼ正しい値に近い推定値が得られた。また、gamma3の係数が0に近く、影響がないという事実に合致する結果が得られた。

ここまでのシミュレーションは、データセットを1つ生成し、それを分析しただけである。1つのデータだけでは偶然そうなっただけかもしれないので、上のシミュレーションを繰り返し行い、推定結果を調べてみよう。

まず、繰り返しシミュレーションを行うための関数を作る。

sim_regression <- function(trials = 200, n = 50, beta = .4) {
    ## Function to simulate multiple regression
    ## Args: trials = num. of simulation trials
    ##       n = sample size of each data set
    ##       beta = effect of the main predictor on the outcome
    ## Return: res = matrix contaning the mean estimates of beta's
    
    z1 <- matrix(runif(trials*n, -5, 5), nrow = n)
    z2 <- matrix(runif(trials*n, -5, 5), nrow = n)
    z3 <- matrix(runif(trials*n, -5, 5), nrow = n)
    x <- 0.2*z1 + rnorm(trials*n)
    mu <- 0.5 + beta*x + 0.2*z1 - 0.4*z2 
    y <- mu + rnorm(trials*n, mean = 0, sd = 1)
    
    res <- matrix(NA, nrow = trials, ncol = 4)
    colnames(res) <- c('true', 'omit1', 'omit2', 'extra')
   
    for (i in 1:trials) {
        df <- data.frame(y = y[,i], x = x[,i], 
                         z1 = z1[,i], z2 = z2[,i], z3 = z3[,i])
        true_fit <- map(true_model, data = df,
                        start = list(
                            alpha = .2,
                            beta = beta,
                            gamma1 = .2,
                            gamma2 = -.4,
                            sigma = 1))
        fit_omit1 <- map(omit_model1, data = df,
                        start = list(
                            alpha = .2,
                            beta = beta,
                            gamma2 = -.4,
                            sigma = 1))
        fit_omit2 <- map(omit_model2, data = df,
                        start = list(
                            alpha = .2,
                            beta = beta,
                            gamma1 = .2,
                            sigma = 1))
        fit_extra <- map(extra_model, data = df,
                         start = list(
                             alpha = .2,
                             beta = beta,
                             gamma1 = .2,
                             gamma2 = -.4,
                             gamma3 = 0,
                             sigma = 1))
 
        res[i, ] <- c(coef(true_fit)[2], coef(fit_omit1)[2],
                      coef(fit_omit2)[2], coef(fit_extra)[2])
    }
    return(res)
}

作った関数を使い、シミュレーションを実行する。xの係数は0.4に設定する。

beta <- .4
set.seed(2016-07-19)
sim1 <- sim_regression(trials = 200, n = 50, beta = beta)

各モデルの係数のMAP推定値の平均値と標準偏差を確認する。

apply(sim1, 2, mean)
##      true     omit1     omit2     extra 
## 0.4070648 0.6591522 0.4156521 0.4027188
apply(sim1, 2, sd)
##      true     omit1     omit2     extra 
## 0.1618894 0.1499448 0.2371405 0.1662764

この結果をみると問題がある（推定値の平均0.4から外れている）のはomit1だけである。 それぞれのモデルから得られたbetaのMAP推定値の分布を図示してみよう。

図のベースを設定。

sim1 <- as.data.frame(sim1)
plt_base <- ggplot(sim1, aes(y = ..density..)) +
    xlab(expression(paste('estimated ', beta)))

正しいモデル。

plt_true <- plt_base + 
    geom_histogram(aes(x = true), color = 'black') +
    geom_vline(xintercept = beta, color = 'red') +
    ggtitle("True Model's Estimates")
print(plt_true)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

交絡変数を除外したモデル。

plt_omit1 <- plt_base + 
    geom_histogram(aes(x = omit1), color = 'black') +
    geom_vline(xintercept = beta, color = 'red') +
    ggtitle("Estimates When We Omit z1")
print(plt_omit1)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

交絡ではない原因を除外したモデル。

plt_omit2 <- plt_base + 
    geom_histogram(aes(x = omit2), color = 'black') +
    geom_vline(xintercept = beta, color = 'red') +
    ggtitle("Estimates When We Omit z2")
print(plt_omit2)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

結果変数の原因ではない余分な変数を加えたモデル。

plt_extra <- plt_base + 
    geom_histogram(aes(x = extra), color = 'black') +
    geom_vline(xintercept = beta, color = 'red') +
    ggtitle("Estimates When We Include z3")
print(plt_extra)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

このシミュレーションから、交絡変数ではない原因を入れ損ねたり、原因ではない変数を入れてしまうのは問題ないが、交絡変数を予測変数に加え忘れると、平均して誤った分析結果を出してしまうことがわかる。したがって、交絡変数は必ず回帰分析に加える必要がある。

交絡を入れ損ねるとバイアスが生じ、関係ない変数を入れても問題がないのであれば、できるだけ多くの変数を統制した方がよさそうである。実際、欠落変数バイアスを防ぐためには、できるだけたくさんの要因を統制した方がよい。ただし、手当たり次第に変数を投入すると起きる問題が、（少なくとも）2つある。

1つ目の問題は、モデルが現実から乖離する確率が大きくなるということである。 この問題が起きるのは、モデルに含む予測変数が増えるに連れて、変数同士の関係のあり方のパタン（例えば、2変数以上の相互作用があるかどうか）が増えるのに対し、実際に正しいモデル（実際にデータが生成される過程）は1つしかないはずだからである。この問題は、ノンパラメトリックな方法を使えば、回避することができる（今回は考えない）。

もう1つの問題は、処置後変数バイアスという異なる種類のバイアスが生じる可能性である。 この問題は、以下のシミュレーションで理解しよう。


処置後変数バイアス

処置後変数バイアス（post treatment variable bias）とは、\(y\) の原因の1つであるが、主な予測変数（因果推論の文脈ではこれを処置変数 [treatment variable] と呼ぶ）である\(x\)の結果でもあるような変数、つまり、\(x\)と\(y\)の関係を仲介するような変数である\(z\)を予測変数として使うことによって生じるバイアスである。 処置後変数バイアスがあると、\(x\)が\(y\)に及ぼす効果を正しく推定することができない。

以下のシミュレーションで、処置後変数バイアスを確認してみよう。

まず、\(x\)が\(y\)に与える効果を設定する。

true_effect <- .25

次に、処置後変数を含むデータ生成過程の関数を作る。

ptb_dgp <- function(n = 100) {
    # Function to generate data to simulate post-treatment bias
    x <- rbinom(n, size = 1, prob = .5)
    z <- .3 + .2*x + rnorm(n, mean = 0, sd = .1)
    y <- 0.2 + (true_effect/0.2)*z  + rnorm(n, mean = 0, sd = .1)
    return(data.frame(y, x, z))
}

試しにデータセットを1つ作る。

set.seed(2016-07-20)
df <- ptb_dgp()
glimpse(df)
## Observations: 100
## Variables: 3
## $ y (dbl) 0.8688067, 0.6859681, 0.6691589, 1.0062167, 0.5793728, 0.766...
## $ x (int) 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, ...
## $ z (dbl) 0.5073898, 0.2416466, 0.4415609, 0.5798076, 0.4432278, 0.454...

ここで、\(x\)と\(z\)の相関関係を確認してみよう。

with(df, cor(x, z))
## [1] 0.707455

\(z\)と\(x\)の相関が高いことがわかる。また、\(y\)と\(z\)については、

with(df, cor(y, z))
## [1] 0.8662613

となり、やはり高い相関を示している。

そこで、「欠落変数バイアスを避けるために」、\(z\)を予測変数に含む以下のモデルを考える。 \[y_i \sim \mbox{Normal}(\mu_i, \sigma)\] \[\mu_i = \alpha + \beta x_i + \gamma z_i\] \[\alpha \sim \mbox{Normal}(0, 10)\] \[\beta \sim \mbox{Normal}(0, 10)\] \[\gamma \sim \mbox{Normal}(0, 10)\] \[\sigma \sim \mbox{Uniform}(0, 10)\]

このモデルのパラメタを推定してみよう。

model_with <- alist(
    y ~ dnorm(mu, sigma),
    mu <- alpha + beta*x + gamma*z,
    alpha ~ dnorm(0, 10),
    beta ~ dnorm(0, 10),
    gamma ~ dnorm(0, 10),
    sigma ~ dunif(0, 10)
)
fit_with <- map(model_with, data = df)
precis(fit_with)
##       Mean StdDev  5.5% 94.5%
## alpha 0.19   0.03  0.13  0.24
## beta  0.04   0.03 -0.01  0.08
## gamma 1.24   0.11  1.07  1.42
## sigma 0.10   0.01  0.09  0.11

MAP推定値の平均値を見ると、\(\beta\)は過少推定、\(\gamma\)は過大推定されている。

ここで、試しに、予測変数 \(z\) を除外して、正しい以下モデルを考えてみよう。

\[y_i \sim \mbox{Normal}(\mu_i, \sigma)\] \[\mu_i = \alpha + \beta x_i\] \[\alpha \sim \mbox{Normal}(0, 10)\] \[\beta \sim \mbox{Normal}(0, 10)\] \[\sigma \sim \mbox{Uniform}(0, 10)\]

このモデルのパラメタを推定しよう。

model_wo <- alist(
    y ~ dnorm(mu, sigma),
    mu <- alpha + beta*x,
    alpha ~ dnorm(0, 10),
    beta ~ dnorm(0, 10),
    sigma ~ dunif(0, 10)
)
fit_wo <- map(model_wo, data = df)
precis(fit_wo)
##       Mean StdDev 5.5% 94.5%
## alpha 0.55   0.02 0.52  0.59
## beta  0.26   0.03 0.22  0.31
## sigma 0.15   0.01 0.13  0.17

正しくないと思っていたモデルなのに、\(\beta\)の推定値は設定した値である0.25 に近い。 これは偶然だろうか？シミュレーションで確かめてみよう。

1回シミュレーションを実行し、処置後変数を含むモデルと含まないモデルの係数のMAP推定値を返す関数を作る。

sim_ptb <- function(true_effect = .25, n = 100) {
    df <- ptb_dgp(n = n)
    fit_with <- map(model_with, data = df,
                    start = list(
                        alpha = 0.5,
                        beta = 0.2,
                        gamma = 1.2,
                        sigma = 0.1))
    fit_wo <- map(model_wo, data = df,
                  start = list(
                      alpha = 0.5,
                      beta = 0.25,
                      sigma = 0.15))
    betas <- c(coef(fit_with)[2], coef(fit_wo)[2])
    names(betas) <- c('with', 'without')
    return(betas)
}

この関数を、replicate() で複数回実行する。

ptb1 <- replicate(200, sim_ptb())

結果をヒストグラムで確認する。

hist(ptb1[1, ], col = 'gray', freq = FALSE,
     main = 'Model with a Post-treatment Variable')

#abline(v = true_effect, col = 'red')
hist(ptb1[2, ], col = 'gray', freq = FALSE,
     main = 'Model without Post-treatment Variable')
abline(v = true_effect, col = 'red')

このように、ある変数\(x\)の効果を推定したいときは、その変数の結果として出てくる変数を統制してはいけない。変数間に時間的な前後関係があれば、このバイアスを回避するのは比較的容易である。しかし、時間的な前後関係が不明なとき、ある変数が交絡変数か処置後変数かを見抜くのは難しい場合がある。統計モデルを作るときには、自分が統制する変数は交絡であり、処置後変数ではないことを理論的に示す必要がある。



授業の内容に戻る