LoginSignup
12
13

More than 3 years have passed since last update.

@saltcooky(saltcooky)

調査観測データの統計科学（二重にロバストな推定/操作変数法）

Last updated at Posted at 2019-03-24

岩波書店の確率と情報の科学シリーズの星野崇宏 著「調査観測データの統計科学」を読んでいき、まとめていってます。
(と言っても、前回からは何ヶ月も経ってますが...)
今回は、前回の続きで第3章「セミパラメトリック解析」で「二重にロバストな推定量」と「操作変数法」についてです。

二重にロバストな推定

傾向スコアに用いた推定では主に二つの改善点が指摘されています。
1. 傾向スコア推定の際には、z=0のグループの共変量の情報を利用していうるが、結果変量の周辺分布の母数推定の際にはz=0のグループの共変量の情報を利用していない。
2. 傾向スコアを計算する際の「割り当て$z$と共変量$x$のモデル」が正しくない場合には、傾向スコアを用いた解析法は間違った結果を与える可能性がある。

この点を解消するために共変量で結果変量を説明する回帰関数を推定方程式に加えることにより、２点は解消できます。
この推定方程式による推定量が一致推定量になるには,次のどちらかの条件が成立していることが必要です。

A. 傾向スコアを算出する「割り当て$z$と共変量$x$のモデル」が正しく指定されている
B. 「共変量で結果変量を説明する回帰関数」が正しく指定されている

どちらか一方でも成り立っていれば、一致推定量がえられることから、二重にロバストな推定量(doubly robust estimator)と呼ばれています。
二重にロバストな推定量はIPW推定量よりも漸近分散が小さい推定量です。
結果変量$y_1,y_0$の期待値$E(y_1), E(y_0)$の二重にロバストな推定量は次に様に定義されます。

\hat E^{DR}(y_1) = \frac{1}{N}
\sum^N 
\left[ 
\frac{z_{i}y_{i1}}{e(x_i,\hat \alpha)} +
\left(1-\frac{z_i}{e(x_i,\hat \alpha)}\right)g(x_i,\hat \beta_1)
\right] \\

= \frac{1}{N}
\sum^N 
\left[ y_{i1}+
\left(\frac{z_i-e(x_i,\hat \alpha)}{e(x_i,\hat \alpha)}\right)(y_{i1}-g(x_i,\hat \beta_1))
\right] 


\hat E^{DR}(y_0)= \frac{1}{N}
\sum^N 
\left[ 
\frac{z_{i}y_{i0}}{e(x_i,\hat \alpha)} +
\left(1-\frac{z_i}{e(x_i,\hat \alpha)}\right)g(x_i,\hat \beta_0)
\right]

ここで、$e(x_i, \hat \alpha)$は一致推定量の母数$\hat \alpha$を用いた際の傾向スコアをを算出するためのモデル、$g(x_i,\hat \beta)$は母数$\beta$の一致推定量$\hat \beta$を用いた、結果変量を説明するための回帰モデル$E(y|x)$のです。

ここで$\hat E^{DR}(y_0)$の期待値は、次に様に計算できます。

E \left[
\hat E^{DR}(y_0)
\right]
= E(y_1) +
E \left[
\frac{z-e(x,\alpha^*)}{e(x,\alpha^*)}(y_{1}-g(x,\beta_1^*))
\right]

第二項は条件Aを満たしているとき、次の様に0となります。

E_{y,x} \left[
E_{z|y,x}
\Bigl\{
\frac{z-e(x,\alpha^0)}{e(x,\alpha^0)}
\Bigr\}
(y_{1}-g(x,\hat \beta_1^*))
\right] \\
= E_{y,x} \left[
\frac{e(x,\alpha^0)-e(x,\alpha^0)}{e(x,\alpha^0)}
(y_{1}-g(x,\hat \beta_1^*))
\right]=0

また、条件Bを満たしているときも、次の様に０となります。

E_{z,x} \left[
E_{y|z,x}
\Bigl\{
\frac{z-e(x,\alpha^*)}{e(x,\alpha^*)}
\Bigr\}
(y_{1}-g(x,\hat \beta_1^0))
\right] \\
= E_{z,x} \left[
\frac{z-e(x,\alpha^*)}{e(x,\alpha^*)}
(E_{y|z,x}(y_{1})-g(x,\hat \beta_1^0))
\right]=0

そのため、２つの条件のどちらかが満たせている場合、結果変量$y_1,y_0$の期待値$E(y_1), E(y_0)$は一致推定量になります。

Rで二重にロバストな推定

Rで二重にロバストな推定を行いたいと思います。
用いたデータはMatchingパッケージに含まれるデータセットlalondeを利用しました。
このデータは1976年に実施された米国職業訓練プログラムを受講した群としていない群で、その後の年収を追跡したものです。
今回は1978年の年収への職業訓練プログラムを受講の効果を推定します。
共変量には年齢、教育年数74年の年収、75年の年収等の10個があります。
はじめに、データを整形します。

library(tidyverse)
library(Matching)
data(lalonde)

n      <- nrow(lalonde)
Y      <- lalonde %>% dplyr::select(re78)
z1     <- lalonde %>% dplyr::select(treat) %>% unlist()
data1  <- lalonde %>% filter(treat==1)
data0  <- lalonde %>% filter(treat==0)

介入の割り当てモデルを作成し、傾向スコアを算出します。

# 割当モデルを作成
logit.model <- glm(treat~age+educ+black+hisp+married+nodegr+re74+re75+u74+u75,
                   data=lalonde, 
                   family=binomial(link = "logit"))
ps     <- logit.model$fitted.values

次に$E(y|x)$のモデルの推定を行います。

# 説明変数への回帰モデルを作成
fomula <- re78~age+educ+black+hisp+married+nodegr+re74+re75+u74+u75
model1 <- lm(fomula, data = data1)
model0 <- lm(fomula, data = data0)

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

二重にロバストな結果変数の周辺期待値を推定します。
そして、最後に因果効果を算出します。

# 二重にロバストな周辺期待値を推定
dre1    <- (1/n)*sum(Y+((z1-ps)/ps)*(Y-fitted1))
dre0    <- (1/n)*sum(((1-z1)*Y)/(1-ps)+(1-(1-z1)/(1-ps))*fitted0)

# 因果効果を推定
dre <- dre1-dre0

# 効果は以下の通りです
# > dre
# [1] 1571.777

通常の回帰分析では1727、IPW推定では1641であったため、これらより因果効果は少し小さめに推定した様です。

操作変数による推定と分析

回帰分析では、説明変数と誤差が無相関であることが、一致推定量であることの条件となります。
しかし、結果変数と誤差を分離できるほどの説明変数を観測することができない場合あり、上記の条件が満たされていない場合も多いです。

そこで、特立変数と相関があり、誤差との相関がゼロである変数を用いることで、条件を満たした推定ができるとされています。
割り当て$d(=1,0)$を説明する操作変数$z=(1,0)$が存在し、$z=1$であれば得られる割り当て変数を$d_1$、$z=0$であれば観測される変数を$d_0$とします。

d=zd_1+(1-z)d_0

また、結果変数についても$d=1$であれば$y_1$、$d=0$であれば$y_0$が得られるとします。

y=dy_1+(1-d)y_0

そして、操作変数の仮定は次に様に表すことができます。

(y_0,y_1) \perp z|d, 　d_z \perp　z

操作変数は直接結果変数に影響を与えず、割り当て変数を通じてもに結果変量に影響を与えるという状態を表現しています。
この状態を除外制約が成立していると言います。

操作変数法を回帰分析に適用する場合、二段階最小二乗法を用いて推定を行います。
$x_1$(内生変数)の$y$への影響を分析する場合、回帰分析の形で表現すると次の様になります。

y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 +\epsilon

そこで、$x_1$と$\epsilon$が無相関となる操作変数$d$を導入します。
$d$と内生変数以外の変数を用いて、内生変数を表現する回帰モデルを考えます。これが１段階目の推定です。

x_1 = \gamma + \gamma x_2 + \delta d + \omega

このモデルにより内生変数の予測値を用います。

\hat x_1 = \gamma + \gamma x_2 + \delta d

この予測値を用いて結果変量を表現する回帰モデルに用いて推定を行います。
これが、二段階目の推定です。

y = \beta_0 + \beta_{iv} \hat x_1 + \beta_2 x_2 +\epsilon

ここで、$\beta_{iv}$が操作変数法による$x_1$の推定された因果効果となります。

次に、操作変数を用いない場合のバイアスとなぜ操作変数法により一致推定量になるのかについてです。
目的変数をy, 観測されていない説明変数を$x_a$、効果を知りたい説明変数$x_1$とした時の回帰モデルを次の様に定義します。

y=\alpha + \beta x_1 + \gamma x_a + \epsilon

この時の回帰分析による偏回帰係数$\hat \beta_s$の期待値は次の様になります。

\begin{align}
E[\hat \beta_s]
&= \frac{Cov(Y,x_1)}{Var(x_1)} \\
&= \frac{Cov(\alpha + \beta x_1 + \gamma x_a + \epsilon,x_1)}{Var(x_1)} \\
&= \frac{Cov(\beta x_1,x_1)}{Var(x_1)} + \frac{\gamma Cov(x_a,x_1)}{Var(x_1)} \\
&= \beta + \frac{\gamma Cov(x_a,x_1)}{Var(x_1)} \\
\end{align}

$x_a$は観測できていない場合のバイアスは第二項部分であり、$x_1$で$x_a$を説明する場合の回帰係数です。
操作変数法による編回帰係数$\beta_{iv}$は次の様に表現でき、変形することができます。

\begin{align}
\beta_{iv}
&= \frac{Cov(y,x_1)}{Cov(x_1,d)} \\
&= \frac{Cov(\alpha + \beta x_1 + \gamma x_a+\epsilon,x_1)}{Cov(x_1,d)} \\
&= \frac{1}{Cov(x_1,d)} \{\beta Cov(x_1,d) + \gamma Cov(x_a,d) + Cov(\epsilon,d)\} \\
\end{align}

除外制約により、$Cov(x_a,d)=0$、かつ$Cov(\epsilon,d)=0$、また$Cov(x_a,d) \neq 0$であることから上記式は次に様になります。

\beta_{iv} = \frac{1}{Cov(x_1,d)} \beta Cov(x_1,d) = \beta

この様に、操作変数法は除外制約を満たしていれば、バイアスのない因果効果を推定することができます。

Rで操作変数法

Rで操作変数法を行いたいと思います。
用いたデータは、矢内さんの記事に習って「計量経済学の第一歩 -- 実証分析のススメ」のサポートページからダウンロードできる修学年数と収入のデータを利用します。

> income_df <- readr::read_csv("./data/csv/8_income.csv")
> str(income_df)
Classes ‘spec_tbl_df’, ‘tbl_df’, ‘tbl’ and 'data.frame':    734 obs. of  8 variables:
 $ exper  : num  19 19 19 19 23 23 23 23 9 9 ...
 $ exper2 : num  361 361 361 361 529 529 529 529 81 81 ...
 $ mbirth : num  6 4 1 9 8 11 3 6 10 11 ...
 $ moyeduc: num  9 12 9 9 9 12 9 9 13 12 ...
 $ payeduc: num  9 12 9 9 9 13 9 9 13 12 ...
 $ sibs   : num  1 1 1 1 3 2 2 6 1 2 ...
 $ yeduc  : num  9 9 9 9 9 9 9 9 11 11 ...
 $ lincome: num  5.7 5.99 6.21 6.21 6.31 ...

はじめに、比較するまでに通常の修学年数で収入を説明する回帰分析を行います。

> # 普通に修学年数と収入の関係を分析してみる
> lm_model <- lm(lincome ~ yeduc, data = income_df)
> summary(lm_model)

Call:
lm(formula = lincome ~ yeduc, data = income_df)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.26331 -0.19445  0.00148  0.20215  1.24404 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 5.387691   0.087018  61.915   <2e-16 ***
yeduc       0.055391   0.006091   9.094   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3377 on 732 degrees of freedom
Multiple R-squared:  0.1015,    Adjusted R-squared:  0.1003 
F-statistic:  82.7 on 1 and 732 DF,  p-value: < 2.2e-16

それでは、操作変数法による二段階推定を行なっていきます。
１段階目の推定は,親の修学年数で対処者の修学年数を説明する回帰モデルを作成します。

> # 1段階目の推定
> lm_model_1 <- lm(yeduc ~ payeduc, data = income_df)
> summary(lm_model_1)

Call:
lm(formula = yeduc ~ payeduc, data = income_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.3640 -1.1819 -0.0685  1.9315  4.8181 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 10.52203    0.35015   30.05   <2e-16 ***
payeduc      0.29554    0.02803   10.54   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.909 on 732 degrees of freedom
Multiple R-squared:  0.1319,    Adjusted R-squared:  0.1307 
F-statistic: 111.2 on 1 and 732 DF,  p-value: < 2.2e-16

次に、二段階目の推定で、先ほどのモデルの予測値$\hat x_1$を用いて、収入を説明する回帰モデルを作成します。

> income_df <- income_df %>% 
+   mutate(pred = predict(lm_model_1))
> 
> lm_model_2 <- lm(lincome ~ pred, data = income_df)
> summary(lm_model_2)

Call:
lm(formula = lincome ~ pred, data = income_df)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.15814 -0.17731  0.01089  0.22816  1.13571 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.75290    0.25005  23.007   <2e-16 ***
pred         0.02956    0.01766   1.674   0.0946 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3556 on 732 degrees of freedom
Multiple R-squared:  0.003813,  Adjusted R-squared:  0.002452 
F-statistic: 2.802 on 1 and 732 DF,  p-value: 0.09459

結果として、修学年数の影響は0.0296と通常の回帰分析の結果より小さく、操作変数法を用いない場合は修学年数の影響を過大評価してしまっていると考えられます。

Rでは、2段階推定による操作変数法を行うことができるAREパッケージがあります。
ivreg関数を用いいれば、２段階推定を一つの関数で行うことができます。
$|$の後の１段階推定目で用いる結果変数と説明変数、$|$の2段階推定目で用いる説明変数を記述します。

# AERパッケージを使った推定
library(AER)
fit_iv_1 <- ivreg(lincome ~ yeduc | yeduc + payeduc, data = income_df)

推定された結果を見てみましょう。

> summary(fit_iv_1)

Call:
ivreg(formula = lincome ~ yeduc | yeduc + payeduc, data = income_df)

Residuals:
      Min        1Q    Median        3Q       Max 
-1.263305 -0.194445  0.001479  0.202150  1.244035 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 5.387691   0.087018  61.915   <2e-16 ***
yeduc       0.055391   0.006091   9.094   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3377 on 732 degrees of freedom
Multiple R-Squared: 0.1015, Adjusted R-squared: 0.1003 
Wald test:  82.7 on 1 and 732 DF,  p-value: < 2.2e-16

0.0553と通常の回帰分析の結果と同程度となりました。
共変量を入れると変わってくるかもしれません。

以上で、二重にロバストな推定量と操作変数法の部分を終了します。

参考

12
13
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
12
13