はじめに

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

1章「調査観測データの統計科学（１）」
2章「調査観測データの統計科学（2）」
3章A 「調査観測データの統計科学（3）」
3章B 「調査観測データの統計科学（二重にロバストな推定/操作変数法）」
3章C 「調査観測データの統計科学（回帰分断デザイン/DID推定））」

今回は、第5章「選択バイアスとその除去」についてです。

選択バイアス

選択バイアスとは、情報を得る際に偏ったデータの一部のみを取得することにより、分析した際の結果に歪み(バイアス)のことをさします。
有名なお話としては、第二次世界大戦の米軍における、帰還した戦闘機の被弾場所を集計についてがあげられます。

はじめに、2つの変数$y_1$,$y_2$が、２変量正規分布に従うとします。

  \left(
    \begin{array}{r}
      y_1  \\
      y_2 
    \end{array}
  \right)

  \sim
  N \left(
  \left(
    \begin{array}{r}
      \mu_1  \\
      \mu_2 
    \end{array}
  \right)
  ,
  \left(
    \begin{array}{rr}
      \sigma_1^2 & \rho \sigma_1 \sigma_2  \\
      \rho \sigma_1 \sigma_2 & \sigma_2^2 
    \end{array}
  \right)
\right)

ここで$\rho$は、$y_1$と$y_2$の相関係数です。
ここで、$y_2$がある閾値$c$を超えた場合のみに$y_1$が観測されるとします。
例えば、$y_2$がある個人の仕事における能力、$y_1$を賃金とするとき、一定の能力以下であれば失業してしまう場合での、観測される賃金の分布などが考えられます。

観測された$y_1$が従う分布は次のようになります。

p(y_1|y_1 > c) = \int_c^{\infty} \frac{p(y_1,y_2)}{p(y_2>c)}dy_2

また、期待値を分散は次のようになります。

E(y_1|y_2>c) = \mu_1 + \rho_1 g(\alpha) \\
V(y_1|y_2>c) = \sigma_1^2(1-\rho^2h(\alpha))

ここで、$g()$と$h()$はそれぞれ標準正規分布の確率密度関数と累積分布関数です、
$y_1$と$y_2$が正の相関を持つ場合に$E(y_1)$より$E(y_1|y_2>c)$の方が大きくなり、分散についても$V(y_1)$より$V(y_1|y_2>c)$の方が小さくなってしまいます。
このように、特定の条件のみの対象のデータだけが観測されている場合に、単純に平均値や分散を計算しても正しい値が得られないこと
があります。
特定の被験者の変数の値のみが観測されるメカニズムを、選択メカニズムと呼びます。
また、選択メカニズムを無視して単純な推定を行うことで生じるバイアスを、選択バイアスと呼びます。

ヘックマンのプロビットモデル

選択バイアスを適切に扱う方法として、選択メカニズムをモデリングすることが多いです。
上記の２変量正規分布に従うモデルを拡張して2つの変数が回帰モデルにしたがっているとします。
具体的には、次の２つの回帰モデルを考えます。
賃金の回帰モデル:賃金はその女性の学歴や卒業後の年数、職種、雇用形態などに依存する。
$y_{i1}$を調査対象$i$における$y_1$の値、$x_1$を$y_1$を説明する独立変数ベクトルとします。

y_{i1}=x_{i1}^t \beta_1 + \epsilon_{i1} \tag{1}

ここで、$\beta_1$は偏回帰係数ベクトルです。

ただし、結果変数は全ての調査対象に関して観測されるのではなく、特定の調査対象者のみ観測されます。
観測されるかどうかは、共変量$x_2$に依存すると考えます。

就業のプロビット回帰モデル(選択方程式)：
就業するかの就業により得られる賃金が就業することにより発生るるコストよりうわまるかで決まると考えられます。
そして、就業によるコストは、結婚の有無や子供の年齢や人数、夫の就業形態などの要因に依存します。

ここで、$y_{i2}$を調査対象者$i$の潜在的な上程変数と考え、共変量の値がこの潜在変数に影響を与える回帰モデルを考えます。

y_{i2}=x_{i2}^t \beta_2 + \epsilon_{i2} \tag{2}

$y_{i2}>0$なら$y_{i1}$が観測される。
ここで、$\epsilon_{i1}$と$\epsilon_{i2}$は多変量正規分布を仮定する。ただし、母数空いて時の識別性から$\epsilon_{i2}$の分散は１とします。
このときのモデルはHeckitモデルやType Ⅱ Tobitモデルと呼ばれ、ヘックマンのプロビット選択モデルとも呼ばれている。

ヘックマンの二段階推定量

ここで、$y_2>0$という条件のもとでの$y_1$の期待値は次式のように表現できます。

\begin{eqnarray}
E(y_{i1}|y_{i2}>0) &=& E(x_{i1}^t \beta_1| y_{i2} > 0) + E(\epsilon_{i1} | y_{i2}>0) \\
  &=& x_{i1}^t \beta_1 + E(\epsilon_{i1} | \epsilon_{i2} > -x_{i2}^t \beta_2) \\
  &=& x_{i1}^t \beta_1 + (\rho \sigma_1) \frac{\phi (x_{i2}^t \beta_2)}{\Phi (x_{i2}^t \beta_2)}
\end{eqnarray}

最後の第２項は、選択バイアスの修正のために付加された項であり、コントロール関数と呼ばれます。

ヘックマンの二段階推定法とは、式(2)をプロビット回帰モデルと考えて推定量$\hat \beta_2$を計算し、その推定量を用いて式(1)を推定する方法です。
具体的には、次のような手順を踏みます。
第一段階(1):$y_1$が観測されるかの目的変数を$x_2$で説明するプロビット回帰分析モデルから$\beta_2$の推定値$\hat\beta_2$と得る。
第一段階(2)：推定値$\beta_2$を用いてコントロール関数を計算する。
第二段階：$y_1$が観測されている対象者について、$y_i1$を$x_i1$と$\hat \lambda_i$に回帰する最小二乗法を行う

ヘックマンの二段階推定における欠点として次にようなものが挙げられています。
1.漸近有効性ンボない推定量である
2.$V(\epsilon_{i1}|y_{i2}>0)$は不均質であるため、推定値の分散の式を修正する必要がある.
3.$\beta_2$の定数項に対応する部分以外の値がある程度の大きさの値になる必要がある.
4.共変量に関して制約があり、$x_2$と$y_1$と相関がない、$x_1$と$x_2$の説明変数が重複する場合は二段階目の推定で多重共線性が生じてしまう.

このような欠点が多いことから実際には、最尤推定法により推定することが多いそうです。

Rでプロビット選択モデル

プロビット選択モデルは、sampleSelectionパッケージを用いると実装できます。
今回用いたデータは、sampleSelectionパッケージに含まれる1987年にMrozが実施した既婚女性の就業に関するパネルデータです。
このデータを用いて、選択バイアスを考慮した女性の収入に影響を与える要因について分析していきます。
はじめに、選択バイアスがある状態の通常の最小二乗法による重回帰モデルを作成します。

library(tidyverse)
library(sampleSelection)

data(Mroz87)
Mroz87 <- Mroz87 %>% 
  mutate(kids = if_else(kids5+kids618 > 0, 1, 0))

# 通常の最小二乗法によるモデル作成
basic_model <- lm(wage ~ exper + I(exper^2)  + educ + city + kids, data = Mroz87)

結果を確認してみます。
収入には、経験と教育が影響していることが考えられます。

> summary(basic_model)

Call:
lm(formula = wage ~ exper + I(exper^2) + educ + city + kids, 
    data = Mroz87)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.1008 -1.7628 -0.4442  1.1372 23.7810 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -4.150929   0.636362  -6.523 1.27e-10 ***
exper        0.188135   0.039068   4.816 1.78e-06 ***
I(exper^2)  -0.003316   0.001276  -2.598  0.00956 ** 
educ         0.415469   0.048707   8.530  < 2e-16 ***
city         0.070307   0.229956   0.306  0.75989    
kids        -0.049009   0.253644  -0.193  0.84684    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.977 on 747 degrees of freedom
Multiple R-squared:  0.1621,	Adjusted R-squared:  0.1565 
F-statistic:  28.9 on 5 and 747 DF,  p-value: < 2.2e-16```

次にヘックマンの二段階推定によるプロビット選択モデルを検討します。

# ヘックマンの二段階推定法によるモデル作成
heckman_model <- selection(lfp ~ educ + age + I(age^2) + kids + faminc,
                          wage ~ exper + I(exper^2)  + educ + city + kids, data = Mroz87 ,method = "2step")

結果を確認してみます。
教育は収入に影響していることは変わりませんでした。
一方で、経験は影響する可能性がないことが示唆されました。
選択バイアスにより適切な推定ができていなかった可能性があります。

> summary(heckman_model)
--------------------------------------------
Tobit 2 model (sample selection model)
2-step Heckman / heckit estimation
753 observations (325 censored and 428 observed)
14 free parameters (df = 740)
Probit selection equation:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -4.157e+00  1.402e+00  -2.965 0.003127 ** 
educ         9.818e-02  2.298e-02   4.272 2.19e-05 ***
age          1.854e-01  6.597e-02   2.810 0.005078 ** 
I(age^2)    -2.426e-03  7.735e-04  -3.136 0.001780 ** 
kids        -4.490e-01  1.309e-01  -3.430 0.000638 ***
faminc       4.580e-06  4.206e-06   1.089 0.276544    
Outcome equation:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.9712003  2.0593505  -0.472    0.637    
educ         0.4170174  0.1002497   4.160 3.56e-05 ***
exper        0.0210610  0.0624646   0.337    0.736    
I(exper^2)   0.0001371  0.0018782   0.073    0.942    
city         0.4438379  0.3158984   1.405    0.160    
Multiple R-Squared:0.1264,	Adjusted R-Squared:0.116
   Error terms:
              Estimate Std. Error t value Pr(>|t|)
invMillsRatio   -1.098      1.266  -0.867    0.386
sigma            3.200         NA      NA       NA
rho             -0.343         NA      NA       NA
--------------------------------------------

次に、最尤推定を行いました。

# 最尤法によるモデル作成
heckman_model_ml <- selection(lfp ~ educ + age + I(age^2) + kids + faminc,
                           wage ~ exper + I(exper^2)  + educ + city + kids, data = Mroz87 ,method = "ml")

得られたモデルを確認します。

> summary(heckman_model_ml)
--------------------------------------------
Tobit 2 model (sample selection model)
Maximum Likelihood estimation
Newton-Raphson maximisation, 5 iterations
Return code 2: successive function values within tolerance limit
Log-Likelihood: -1581.258 
753 observations (325 censored and 428 observed)
13 free parameters (df = 740)
Probit selection equation:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -4.120e+00  1.401e+00  -2.942 0.003368 ** 
educ         9.528e-02  2.315e-02   4.115  4.3e-05 ***
age          1.840e-01  6.587e-02   2.794 0.005345 ** 
I(age^2)    -2.409e-03  7.723e-04  -3.119 0.001886 ** 
kids        -4.506e-01  1.302e-01  -3.461 0.000568 ***
faminc       5.680e-06  4.416e-06   1.286 0.198782    
Outcome equation:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.9630242  1.1982209  -1.638    0.102    
educ         0.4570051  0.0732299   6.241 7.33e-10 ***
exper        0.0278683  0.0615514   0.453    0.651    
I(exper^2)  -0.0001039  0.0018388  -0.056    0.955    
city         0.4465290  0.3159209   1.413    0.158    
   Error terms:
      Estimate Std. Error t value Pr(>|t|)    
sigma   3.1084     0.1138  27.307   <2e-16 ***
rho    -0.1320     0.1651  -0.799    0.424    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
--------------------------------------------

こちらも、教育の影響はありますが、経験による影響はないようです。

プロビット選択モデルが有する問題点

プロビット選択モデルには次にような問題点が挙げられています。
1.回帰関数の指定を誤ると推定値に大きなバイアスが生じる。
2.誤差の分布が２変量正規分布を仮定しているため、誤差分布へのロバスト性がなく、分布仮定のチェックもできない。
3.選択バイアスの補正における母数$\rho$の推定が不安定。

因果効果推定にプロビット選択モデルを適用する

プロビット選択モデルに潜在的な結果変数の考え方を導入すれば、因果効果の推定に利用できるようになります。
$y_{i2}>0$なら「介入を加えた場合の結果変数$y_{i1}$が観測されるだけではなく、介入を与えられないばいいの結果変数$y_{i0}$は観測されない」、一方で、$y_{i2} \leq 0$ならば「$y_{i2}>0$は観測され、$y_{i1}>0$は観測されない」とする。さらに結果変数$y_{i1}$,$y_{i0}$についての回帰モデルを仮定する。

y_{i1}=x_{i1}^t\beta_1+\epsilon_{i1}\\
y_{i0}=x_{i1}^t\beta_0+\epsilon_{i0}

  \left(
    \begin{array}{r}
      \epsilon_0  \\
      \epsilon_1  \\
      \epsilon_2
    \end{array}
  \right)

  \sim
  N \left(
  \left(
    \begin{array}{r}
      0  \\
      0  \\
      0 
    \end{array}
  \right)
  ,
  \left(
    \begin{array}{rr}
      \sigma_0^2 & \sigma_{01} & \rho \sigma_0 \\
      \sigma_{02} & \sigma_1^2 & \rho \sigma_1 \\
      \rho \sigma_0& \rho \sigma_1 & 1 
    \end{array}
  \right)
\right)

ここで関心のある因果効果や処置軍での因果効果TETは、次のようになる。

E(y_1-y_0)={E(x_1)}^t(\beta_1-\beta_0)

TET=E(y_1-y_0|z=1)={E(x_1)}^t(\beta_1-\beta_0)+(\rho \sigma_1-\rho \sigma_0)E \left( \left. \frac{\phi (x_{i2}^t \beta_2)}{\Phi (x_{i2}^t \beta_2)}  \right| z=1 \right)

ルービンの因果モデルを選択バイアスのモデリングの選択的な結果変数に対して適用したものであると考えることができます。

強く無視できる割り当て条件

欠測データの枠組みを利用して、欠測インディケータをモデリングすることを考えます。
例えば、市場調査での無回答のバイアスを扱います。
調査対象者N人分において、調査に回答してくれた対象者は$z=1$、拒否または不在で調査できないなら$z=0$となる欠測インディケーター$z$を導入します。
$z$が目的変数$y$と独立な変数であれば、母集団における$y$の期待値と調査に回答していくれた場合に限定した場合の$y$の期待値は等しくなります。

E(y)=E(y|z=1)=E(y|z=0)

しかし、$y$と$z$が独立でなければ、$\hat y_1$の推定量にはバイアスが存在してしまいます。
調査の実施内容により回答する人に偏りがある場合は、一般的に存在することが考えられます。
そのため、目的変数と観測されるかどうかが独立でないとなり、バイアスが存在してしまします。
このような場合、因果推論と同様に共変量を利用して調整を行うことを考えます。

はじめに、「$y$が観測されるかどうかは$x$の値に依存する=強く無視できる割り当て条件」という次の条件が成立すると仮定します。

p(z|x,y)=p(z|x)

期待値に関しては次のように表現できます。

E(y)=E_x(E(y|X))=E_x(E(y|x,z=1))

ここで、$y$の$x$への回帰モデル$y=g(x|\theta)+e$がわかっている場合に、推定する期待値は次にようなる。

\hat E(y)=\frac{1}{N} \sum_i^N \{z_iy_i+(1-z_i)g(x|\theta)\}

傾向スコアや二重にロバストな推定方、IPW推定量などを利用することでバイアスのない推定を行うことができます。
これらの推定方法は、ヘックマンの選択モデルと異なりモデル仮定が非常に弱いことやヘックマンの選択モデルでは少数の共変量を扱う状況を考えているのに対して、傾向スコアを用いる方法では多数の共変量を利用することが出来るといった利点があげられます。

Rで傾向スコアによる選択バイアス補正

上記で扱った女性の就業と収入のデータにおける、収入の期待値を推定します。
傾向スコアの推定にはロジスティック回帰分析、期待値の推定には二重にロバストな推定法を用います。
はじめに、傾向スコアを算出して共変量の分布の乖離がない状態になっていることを確認します。

library(tidyverse)
library(sampleSelection)
library(WeightIt)
library(cobalt)

data(Mroz87)
Mroz87 <- Mroz87 %>% 
  mutate(kids = if_else(kids5+kids618 > 0, 1, 0))

# 傾向スコアの推定
ps_model <- weightit(lfp ~ exper + I(exper^2)  + educ + city + kids, 
                      data = Mroz87,
                      method = "ps",
                      estimand = "ATE")

# 共変量のバランス確認
love.plot(ps_model, threshold = 0.1)

共変量のバランスは取れているようです。

続いて、収入の期待値の二重にロバストな推定量を求めます。

n      <- nrow(Mroz87)
Y      <- Mroz87 %>% pull(wage)
z1     <- Mroz87 %>% pull(lfp) %>% unlist()

data1  <- Mroz87 %>% filter(lfp==1)
data0  <- Mroz87 %>% filter(lfp==0)

# 説明変数への回帰モデルを作成
fomula <- wage ~ exper + I(exper^2)  + educ + city + kids
model1 <- lm(fomula, data = data1)
model0 <- lm(fomula, data = data0)

fitted1 <- predict(model1, Mroz87)
fitted0 <- predict(model0, Mroz87)

# 二重にロバストな推定
dre1    <- (1/n)*sum(Y+((z1-ps_model$ps)/ps_model$ps)*(Y-fitted1))

> dre1
[1] 3.963792

共変量シフト

共変量シフトとは、機械学習の分野でよく議論されている内容です。
機械学習では入力$\omega$から出力$y$へ精度よく予測することが、主な観点となっている。
これは、予測の問題は条件付き分布$p(y|\omega)$の選択と推定の問題に帰着する。
しかし、データを「条件付き分布のモデル選択、および母数推定のために利用する訓練データA」と「これらから得られた条件付き分布の性能をテストするテストデータB」の２つのデータに分けて考える。
この時、入力を与えられた時の出力の条件付き分布$p(y|\omega)$は訓練データAとテストデータBで共通であるとが前提となっています。
つまり、訓練データAとテストデータBは同じ生成規則により得られていることが前提である。
しかし、必ずしも訓練データAとテストデータBは同じ生成規則により得られているとは限らず、このような状態を共変量シフトと呼びます。
共変量シフトの状態での予測問題では、テストデータBでの予測精度をよくすることを考えます。

これまでと同様に$z=1$を訓練データへの割り当て、$z=0$をテストデータへの割り当てであると考えれます。
この時の入力分布は$p(\omega|z=1)=p(\omega|z=0)$です。
訓練データのペア$(y,\omega)$は、同時分布$p(y|\omega)p(\omega|z=0)$からの無作為抽出標本であり、そのうち訓練データは次にような確率で$(y,\omega)$が観測されていると考える。

p(z=1|\omega)=\frac{p(\omega|z=1)p(z=1)}{p(\omega)}=\frac{p(\omega|z=1)p(z=1)}{p(\omega|z=0)p(z=0)+p(\omega|z=1)p(z=1)}

$\omega$を条件付けることで、$y$はランダムな欠測または観測値による選択が行われていると考えることができる。
条件付き分布$p(y|\omega)$の母数推定だけに関心があれば訓練データに対して最尤推定や最小二乗推定などを行えば良い。

本来関心があるのは$z=0$での同時分布$p(y|\omega)p(\omega|z=0)$とすると

予測値$\hat y$と実データ$y$の損失関数を$l(\hat y,y)$と表現すると、共変量シフト下で学習を行う際の期待損失は次にように重みを付けた形となり、テストデータの期待損失と等しくなる。

\begin{eqnarray}
E_{p(y|\omega)p(\omega|z=1)} \left[ \frac{p(\omega|z=0)}{p(\omega|z=1)}l(\hat y,y)\right]\\
&=&\int p(y|\omega)p(\omega|z=1)\frac{p(\omega|z=0)}{p(\omega|z=1)}l(\hat y,y)dyd\omega\\
&=&\int p(y|\omega)p(\omega|z=0)l(\hat y,y)dyd\omega\\
&=&E_{p(y|\omega)p(\omega|z=0)}[l(\hat y,y)]
\end{eqnarray}

この損失関数は重要度重み付きERM()と呼ばれる。
テキストでは、$y$の$\omega$への回帰関数が２次以上の高次であり、線形回帰モデルを用いて予測を行う場合には、訓練データから単純に最小二乗法によって得られた回帰直線を用いてテイストデータのyの予測を行うよりも
、上記の重みつけを行なった最小二乗法によて計算された回帰直線でテストデータを予測した方が精度が良くなると紹介しています。

実際に推定を行う際には、${p(\omega|z=0)}$と${p(\omega|z=1)}$をそれぞれ推定せず、$\frac{p(\omega|z=0)}{p(\omega|z=1)}$で推定することが推奨されている。
$\frac{p(\omega|z=0)}{p(\omega|z=1)}$は密度比と呼ばれています。
密度比の推定には,KMM (Kernel Mean Matching),Logistic Regression (LogReg), Least-Square Importance Fitting(LSIF), uLSIF (unconstrained LSIF)といったアルゴリズムが存在する。

Rで共変量シフト

今回、真の関数は次のような二次関数であると仮定します。

y = 0.1 (x-5)^2+5

この関数の出力に対して、$N(0,2.4)$のバラツキがあります。
そして、この関数の入力になる$x$の値は、学習とテストそれぞれ次のような正規分布に従い発生しているとしました。

x_{train} \sim N(0,4),  x_{test} \sim N(9,2)

library(tidyverse)

# 二次関数
houteishiki <- function(x){
  y = 0.1* (x-5)^2 + 5 + rnorm(length(x), 0, 2.4)
  return(y)
}

# データ生成
test_point <- rnorm(100, 9, 2)
train_point <- rnorm(200, 0, 4)
base_line <- data.frame(x = seq(-15,15,0.1)) %>% 
  mutate(y = houteishiki(x))

test_df <- data.frame(x = test_point, 
                      y = test_point %>% houteishiki(),
                      type = "test") 

train_df <- data.frame(x = train_point, 
                       y = train_point %>% houteishiki(), 
                       type = "train") 
test_df <- test_df 
train_df <- train_df 
all_df <- bind_rows(test_df,train_df)

データを確認してみます.
学習とテストのデータが明確に異なることが確認できます。

all_df %>% 
  ggplot(aes(x = x, y = y, alpha = 1))+
  geom_point(aes(col = type))+
  stat_smooth(method = loess, formula = y ~ x,linetype = "dashed", se = FALSE, size=0.6)+
  theme_bw()

密度比推定にはdensratioパッケージを用います。
推定アルゴリズムは、デフォルトのuLSIFです。

library(densratio)

# 密度比推定
dens_rate <- densratio(test_df$x, train_df$x)
train_dens_rate <- dens_rate$compute_density_ratio(train_df$x)

# 重み付き回帰分析
train_model <- lm(y ~ x, data = train_df) #普通の回帰分析
cs_model <- lm(y ~ x, data = train_df, weights = (train_dens_rate))

base_line <- base_line %>% 
  mutate(model1 = predict(train_model, base_line),
         model2 = predict(cs_model, base_line))

結果を確認してみます。

all_df %>% 
  ggplot(aes(x = x, y = y, alpha = 1))+
  geom_point(aes(col = type))+
  ylim(-5,45)+
  stat_smooth(method = loess, formula = y ~ x,linetype = "dashed", se = FALSE, size=0.6)+
  geom_line(base_line, mapping=aes(x = x, y = model1))+
  geom_line(base_line, mapping=aes(x = x, y = model2))+
  theme_bw()

学習データをそのまま用いたモデルでは、学習データに対する誤差は小さいですが、テストデータに対しての誤差は大きいです.
テストデータに対して比較的誤差のない推定ができていると考えられます。

今回はここまでで終わりたいと思います。

参考