欠落変数バイアス

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

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

交絡変数は、必ず統制する必要がある（厳密には、バックドアパスをブロックすればよい）。交絡変数を統制しないと、推定にバイアスが生じる。このバイアスを欠落変数バイアス (omitted variable bias) と呼ぶ。経済学では、選択バイアス（selection bias） とも呼ばれる。

\(y\)と\(x\)の両者に影響を与える\(z\)という変数があるとき、\(z\)を無視して、 \[y = \beta_0 + \beta_1 x + u\] という式を考えると、\(z\) は誤差項\(u\)に含まれることになる。 そうすると、当然ながら、説明変数\(x\)と誤差項\(u\)の間に相関があるので、最小二乗推定量が不偏性もつための仮定（教科書 p.122の (3)）が満たされず、\(\beta_1\) の推定が偏ったもの（つまり、バイアスをもったもの）になってしまうのである。

このことを、シミュレーションで確認してみよう。 まず、データを作る。ここでは、x, z1, z2 という3つの変数がyの原因となっている（muを計算する行を参照）。また、z1はxの原因でもあるので、z1は交絡変数である。したがって、z1を統制(（コントロール）しないと、xの係数が正しく推定できないはずである。z3は結果変数とは無関係の変数である。

# シミュレーションを再度実行したときに同じ結果が得られるように、乱数の種を指定する
# 異なる結果を得たいときは、set.seed() の中身を変える
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 <- 1.5 + 0.4*x + 0.5*z1 - 0.2*z2 
y <- rnorm(n, mean = mu, sd = 1)
df <- tibble(y, x, z1, z2, z3)

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

true_model <- lm(y ~ x + z1 + z2, data = df)
tidy(true_model)

## # A tibble: 4 x 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)    1.58     0.114      13.8  1.55e-24
## 2 x              0.419    0.104       4.04 1.07e- 4
## 3 z1             0.507    0.0432     11.7  3.19e-20
## 4 z2            -0.159    0.0415     -3.83 2.26e- 4

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

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

omitted_1 <- lm(y ~ x + z2, data = df)
tidy(omitted_1)

## # A tibble: 3 x 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)    1.55     0.177       8.77 6.16e-14
## 2 x              0.983    0.142       6.91 5.14e-10
## 3 z2            -0.157    0.0644     -2.44 1.64e- 2

このモデルでは、xの係数の推定値が 0.98 になり、xのyに対する影響がかなり過大に推定されている。

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

omitted_2 <- lm(y ~ x + z1, data = df)
tidy(omitted_2)

## # A tibble: 3 x 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)    1.62     0.121      13.3  1.18e-23
## 2 x              0.384    0.110       3.48 7.48e- 4
## 3 z1             0.506    0.0462     11.0  1.16e-18

ここでは、xの係数は、正しい値である0.4に近い。

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

extra_model <- lm(y ~ x + z1 + z2 + z3, data = df)
tidy(extra_model)

## # A tibble: 5 x 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)  1.57       0.115     13.7   2.82e-24
## 2 x            0.422      0.105      4.02  1.17e- 4
## 3 z1           0.507      0.0435    11.7   5.21e-20
## 4 z2          -0.159      0.0417    -3.81  2.46e- 4
## 5 z3          -0.00954    0.0395    -0.241 8.10e- 1

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

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

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

sim_regression <- function(trials = 200, n = 50, beta = .4) {
    ## 重回帰分析をシミュレートするための関数
    ## 引数：trials = シミュレーションの繰り返し回数（既定値は200)
    ##       n = 標本サイズ（既定値は50）
    ##       beta = x 係数 (beta) の母数（パラメタ）
    ## 返り値：res = 係数の推定値を要素にもつ行列
    
    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 <- 1.5 + beta*x + 0.5*z1 - 0.2*z2 
    y <- mu + rnorm(trials*n, mean = 0, sd = 1)
    
    beta_hat <- matrix(NA, nrow = trials, ncol = 4)

    colnames(beta_hat) <- c('true', 'omit1', 'omit2', 'extra')
   
    for (i in 1:trials) {
        df <- tibble(y = y[,i], x = x[,i], 
                     z1 = z1[,i], z2 = z2[,i], z3 = z3[,i])
        true_fit <- lm(y ~ x + z1 + z2, data =df)
        fit_omit1 <- lm(y ~ x + z2, data = df)
        fit_omit2 <- lm(y ~ x + z1, data = df)
        fit_extra <- lm(y ~ x + z1 + z2 + z3, data = df)  
        beta_hat[i, ] <- c(coef(true_fit)[2], coef(fit_omit1)[2],
                           coef(fit_omit2)[2], coef(fit_extra)[2])
    }
    return(beta_hat)
}

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

beta <- 0.4
set.seed(2019-05-01)
sim1 <- sim_regression(trials = 200, n = 50, beta = beta)

各モデルの係数の最小二乗推定値の平均値を確認する。

apply(sim1, 2, mean)

##      true     omit1     omit2     extra 
## 0.4011777 1.0318576 0.4069522 0.4006314

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

sim1_beta <- as_tibble(sim1)  # ggplot2 を使うためにデータフレームに変換
# 全モデルに共通する図のベースを作る
plt_base <- ggplot(sim1_beta, aes(y = ..density..)) +
  xlab(expression(paste(beta, "の推定値", sep = ""))) +
  ylab("確率密度")

正しいモデル。

plt_true <- plt_base + 
    geom_histogram(aes(x = true), color = 'black', bins = 12) +
    geom_vline(xintercept = beta, color = 'red') +
    ggtitle("正しいモデルの推定値の分布")
plot(plt_true)

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

plt_omit1 <- plt_base + 
    geom_histogram(aes(x = omit1), color = 'black', bins = 12) +
    geom_vline(xintercept = beta, color = 'red') +
    ggtitle("交絡変数を除外したモデルの推定値の分布")
plot(plt_omit1)

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

plt_omit2 <- plt_base + 
    geom_histogram(aes(x = omit2), color = 'black', bins = 12) +
    geom_vline(xintercept = beta, color = 'red') +
    ggtitle("交絡ではない原因を除外したモデルの推定値の分布")
plot(plt_omit2)

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

plt_extra <- plt_base + 
    geom_histogram(aes(x = extra), color = 'black', bins = 12) +
    geom_vline(xintercept = beta, color = 'red') +
    ggtitle("余分な変数を追加したモデルの推定値の分布")
plot(plt_extra)

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

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

まず、モデルが現実（真実）から乖離する確率が大きくなる。 この問題が起きるのは、モデルに含む説明変数が増えるにつれて、変数同士の関係のあり方のパタン（例えば、2変数以上の相互作用があるかどうか）が増えるのに対し、実際に正しいモデル（実際にデータが生成される過程）は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) {
    # 処置後変数バイアス問題をシミュレートする関数
    x <- rbinom(n, size = 1, prob = .5)
    z <- 0.3 + 0.2 * x + rnorm(n, mean = 0, sd = 0.1)
    y <- 0.2 + (true_effect / 0.2) * z  + rnorm(n, mean = 0, sd = 0.1)
    return(tibble(y, x, z))
}

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

set.seed(2018-04-02)
df <- ptb_dgp()
glimpse(df)

## Observations: 100
## Variables: 3
## $ y <dbl> 0.4259625, 0.9550180, 0.3682632, 0.6271188, 0.9494557, 0.77470…
## $ x <int> 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1,…
## $ z <dbl> 0.2257827, 0.5996252, 0.1408864, 0.2506434, 0.7002162, 0.42286…

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

with(df, cor(x, z))

## [1] 0.6664344

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

with(df, cor(y, z))

## [1] 0.8400825

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

そこで、「欠落変数バイアスを避けるために」、\(z\)を予測変数に含む以下のモデルを考える。 \[y_i = \alpha + \beta x_i + \gamma z_i + u.\] このモデルのパラメタを推定してみよう。

fit_with <- lm(y ~ x + z, data = df)
tidy(fit_with)

## # A tibble: 3 x 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)   0.230     0.0336     6.84  7.07e-10
## 2 x             0.0261    0.0292     0.892 3.75e- 1
## 3 z             1.13      0.104     10.8   2.25e-18

推定値を見ると、xの係数\(\beta\)は過小推定（本当は 0.25）されている。

ここで、説明変数 \(z\)を除外して、以下のモデルを考えてみよう。 \[y_i = \alpha + \beta x_i + u.\] このモデルのパラメタを推定しよう。

fit_wo <- lm(y ~ x, data = df)
tidy(fit_wo)

## # A tibble: 2 x 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)    0.550    0.0234     23.5  5.22e-42
## 2 x              0.237    0.0322      7.35 5.89e-11

「誤って」コントロールすべきと考えていた \(z\) を除外したモデルなのに、\(\beta\)の推定値は設定した値である0.25 に近い。 これは偶然だろうか？シミュレーションで確かめてみよう。

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

sim_ptb <- function(true_effect = .25, n = 100) {
    df <- ptb_dgp(n = n)
    fit_with <- lm(y ~ x + z, data = df)
    fit_wo <- lm(y ~ x, data = df)
    betas <- c(coef(fit_with)[2], coef(fit_wo)[2])
    names(betas) <- c('with', 'without')
    return(betas)
}

この関数を、replicate() で複数回実行する。500回繰り返してみよう。

ptb_1 <- replicate(500, sim_ptb())

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

dd <- tibble(with = ptb_1[1, ],
             without = ptb_1[2, ])  # ggplot で使うためのデータセット
hist_with <- ggplot(dd, aes(x = with, y = ..density..)) +
  geom_histogram(color = "black", bins = 12) +
  geom_vline(xintercept = true_effect, color = "red") +
  labs(x = "処置後変数を含むモデルの推定値", y = "確率密度")
plot(hist_with)

hist_wo <- ggplot(dd, aes(x = without, y = ..density..)) +
  geom_histogram(color = "black", bins = 12) +
  geom_vline(xintercept = true_effect, color = "red") +
  labs(x = "処置後変数を含まない（正しい）モデルの推定値", y = "確率密度")
plot(hist_wo)

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