具体例を用いた交絡の確認

# はじめに
本稿では"**A Crash Course in Good and Bad Controls**"というタイトルの文献(以下リンク)について, 簡単な具体例を用いてゆっくり確認していきます. 

https://ftp.cs.ucla.edu/pub/stat_ser/r493.pdf

なお, 本稿で頻出する**control**(statistical control)という用語は, 私が調べる限り, 関心下の説明変数が目的変数に与える効果を検証する際に, 関心のない変数の影響を固定するために敢えて説明変数に追加されるもののことを指すようです. 



#  関連ライブラリのインポート

はじめに, 以下のように関連ライブラリをインポートしておきます.

```python
import numpy as np
import graphviz
from scipy import stats 
import statsmodels.api as sm 

import matplotlib.pyplot as plt 
plt.style.use('fivethirtyeight')


# Define the sample size 
N_SAMPLES = 1000
```

# Good controls
### Model.1; good control(blocking back-door paths)

まずはじめに, Pythonのライブラリ`graphviz`を用いて, DAGを描画します.

以下のようなコードを実行すると,

```python
# Create a directed graph
g_1 = graphviz.Digraph(format='png')

## Add nodes
nodes_1 = ['Z', 'X', 'Y']
[g_1.node(n) for n in nodes_1] 

g_1.edges(['ZX', 'ZY'])

# Render for print 
g_1.render('img/graph_01')
```

下図のようなDAG

![graph_01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/62e8258d-c487-c797-033a-d448a406e565.png)

を作ることができました. ここで, Y, X, Zはそれぞれ被説明変数, 関心下の説明変数, コントロールしたい変数を表します.

つまり, **Good controls**, **Bad controls**のいずれに該当するのか判断したいのは変数Zです.

実際に**Good control**であることを確認したいので, このDAGに対して, 次のような**Data generating process**(DGP)を仮定しましょう.

```math
\begin{align} \\
&Z \sim \mathcal{N}(0, 1) \\
&X = 2z + 0.5\epsilon_x, \quad \epsilon_x \sim \mathcal{N}(0, 1) \\
&Y = -0.1z + 0.5\epsilon_y, \quad \epsilon_y \sim \mathcal{N}(0, 1) \\
\end{align}
```

単純化のため, 線形のDGPとしました.

このDGPに従って, 以下のコードのように変数をアサインします.

```python
# Set random seed for reproductivity 
np.random.seed(45)

# assign values to the variables: 
z = np.random.randn(N_SAMPLES)
x = 2 * z + 0.5 * np.random.randn(N_SAMPLES)
y = -0.1 * z + 0.5 * np.random.randn(N_SAMPLES)
```

次に, 線形重回帰を行います. 2種類のモデル

-   y \~ x

-   y \~ x + z

を作成します.

次に, 作成したモデルに使用したそれぞれの説明変数について以下のメトリクス

-   偏回帰係数

-   偏回帰係数を0とする帰無仮説のP-値

-   有意水準5%としたとき, 偏回帰係数が有意かどうか

を取得し, 表示させます.

それでは以下のとおりコードを実行して結果を出力させます.

```python
# Define four model variants 
variants = [
    [x], 
    [x, z]
]

# Fit models iteratively and store the results 
results = []
for variant in variants: 
    X = sm.add_constant(np.stack(variant).T)

    # Instantiate the model and fit it 
    model = sm.OLS(y, X)
    fitted_model = model.fit()

    results.append([fitted_model.params, fitted_model.pvalues])

    print(f'Params: {fitted_model.params}')
    print(f'p-vals: {fitted_model.pvalues}')
    print(f'Signif: {fitted_model.pvalues <= .05}\n')
```
```
## Params: [-0.01044718 -0.04884977]
## p-vals: [5.13973229e-01 9.49711998e-10]
## Signif: [False  True]
## 
## Params: [-0.01107309  0.01410072 -0.13272288]
## p-vals: [0.48847265 0.65651128 0.04055804]
## Signif: [False False  True]
```

上記出力結果のうち, 説明変数Xは2列目の要素に対応します.

1つ目のモデル,

-   y \~ x

では, Xの偏回帰係数は有意でした.

一方で, 2つ目のモデル,

-   y \~ x + z

においては帰無仮説を棄却することができませんでした.

また, 1つ目のモデルにおいて, Xの偏回帰係数の推定値は-0.04884977でした.

一方で, 以下に再掲するDGPにおいて, 

```math
\begin{align} \\
&Z \sim \mathcal{N}(0, 1) \\
&X = 2z + 0.5\epsilon_x, \quad \epsilon_x \sim \mathcal{N}(0, 1) \\
&Y = -0.1z + 0.5\epsilon_y, \quad \epsilon_y \sim \mathcal{N}(0, 1) \\
\end{align}
```

確かに, YをXだけを用いてで回帰した時の偏回帰係数は-0.05になりそうな気がします.

以上のことから, このケースでは回帰式のControlとしてZを追加することによって, 偏回帰係数上で生じていた交絡による擬相関を除去できたということができます.

### Model.2; good control(blocking back-door paths)

次に, 以下のようなDAG

```python
# Create a directed graph
g_2 = graphviz.Digraph(format='png')

## Add nodes
nodes_2 = ['Z', 'X', 'Y', 'U']
[g_2.node(n) for n in nodes_2] 

g_2.edges(['ZX', 'UZ', 'UY'])

# Render for print 
g_2.render('img/graph_02')
```

![graph_02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/8554f9fd-2493-d1b4-549e-84d859a8da2e.png)

を仮定します.

先ほどのDAGとは異なり, 観測できない(Unobservable)変数Uを追加しました.

更に, DGPとして,

```math
\begin{align} \\
&U \sim \mathcal{N}(0, 1) \\
&Z = u + \epsilon_z, \quad \epsilon_z \sim \mathcal{N}(0, 1) \\
&X = 2z + 0.5 \epsilon_x, \quad \epsilon_x \sim \mathcal{N}(0, 1) \\
&Y = -0.1u + 0.5 \epsilon_y, \quad \epsilon_y \sim \mathcal{N}(0, 1)
\end{align}
```

を仮定します.

そして, 先ほどと同様, YをXのみで回帰するモデル, YをXとZで回帰するモデルを作成し, メトリクスを出力します.

```python
# Set random seed for reproductivity 
np.random.seed(45)

# assign values to the variables: 
u = np.random.randn(N_SAMPLES)
z = 1 * u + np.random.randn(N_SAMPLES)
x = 2 * z + 0.5 * np.random.randn(N_SAMPLES)
y = -0.1 * u + 0.5 * np.random.randn(N_SAMPLES)

# Define four model variants 
variants = [
    [x], 
    [x, z]
]

# Fit models iteratively and store the results 
results = []
for variant in variants: 
    X = sm.add_constant(np.stack(variant).T)

    # Instantiate the model and fit it 
    model = sm.OLS(y, X)
    fitted_model = model.fit()

    results.append([fitted_model.params, fitted_model.pvalues])

    print(f'Params: {fitted_model.params}')
    print(f'p-vals: {fitted_model.pvalues}')
    print(f'Signif: {fitted_model.pvalues <= .05}\n')
```
```
## Params: [ 0.0260086  -0.01964127]
## p-vals: [0.09582058 0.00041798]
## Signif: [False  True]
## 
## Params: [ 0.02671423  0.05317639 -0.15057924]
## p-vals: [0.08646669 0.08515535 0.01663936]
## Signif: [False False  True]
```
モデル1と同様, 1つ目のモデルでは, Xの偏回帰係数は有意である一方で, 2つ目のモデルにおいては帰無仮説を棄却することができませんでした.

上図のDAGによりデータが生成される場合, 未観測の変数Uが存在したとしても, Zを観測できていれば, XがYに与える関係を正しく推論できるということになります.

### Model.3; good control(blocking back-door paths)

モデル1, モデル2と類似した下図のようなDAG(と当該DAGに準拠したDGP)についても同様のことが言えます.

```python
# Create a directed graph
g_3 = graphviz.Digraph(format='png')

## Add nodes
nodes_3 = ['Z', 'X', 'Y', 'U']
[g_3.node(n) for n in nodes_3] 

g_3.edges(['ZY', 'UZ', 'UX'])

# Render for print 
g_3.render('img/graph_03')
```

![graph_03.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/b3a6f49d-d60d-19d7-5321-8798cf32b8db.png)


```math
\begin{align} \\
&U \sim \mathcal{N}(0, 1) \\
&Z = 1 \cdot u + \mathcal{N}(0, 1) \\
&X = 2 \cdot u + 0.5 \cdot \mathcal{N}(0, 1) \\
&Y = -0.1 \cdot z + 0.5 \cdot \mathcal{N}(0, 1) \\
\end{align}
```

前述の手順を実施すると, 以下のようになり,

```python
# Set random seed for reproductivity 
np.random.seed(45)

# assign values to the variables: 
u = np.random.randn(N_SAMPLES)
z = 1 * u + np.random.randn(N_SAMPLES)
x = 2 * u + 0.5 * np.random.randn(N_SAMPLES)
y = -0.1 * z + 0.5 * np.random.randn(N_SAMPLES)

# Define four model variants 
variants = [
    [x], 
    [x, z]
]

# Fit models iteratively and store the results 
results = []
for variant in variants: 
    X = sm.add_constant(np.stack(variant).T)

    # Instantiate the model and fit it 
    model = sm.OLS(y, X)
    fitted_model = model.fit()

    results.append([fitted_model.params, fitted_model.pvalues])

    print(f'Params: {fitted_model.params}')
    print(f'p-vals: {fitted_model.pvalues}')
    print(f'Signif: {fitted_model.pvalues <= .05}\n')
```

```
## Params: [ 0.0226853  -0.04019455]
## p-vals: [1.50171604e-01 2.32586886e-07]
## Signif: [False  True]
## 
## Params: [ 0.02420901  0.00486855 -0.09988378]
## p-vals: [1.16764648e-01 6.30741837e-01 3.83785779e-11]
## Signif: [False False  True]
```

Xに関する偏回帰係数のP-値について, Zが**Goood control**であることと整合的です.

### Model.4; good control

モデル4, モデル5, モデル6ではXとYの間に**Mediator**(**Chain**)が存在するケースについて考察します.　 まずはじめに, 以下のようDAG

```python
# Create a directed graph
g_4 = graphviz.Digraph(format='png')

## Add nodes
nodes_4 = ['Z', 'X', 'Y', 'M']
[g_4.node(n) for n in nodes_4] 

g_4.edges(['ZX', 'ZM', 'XM', 'MY'])

# Render for print 
g_4.render('img/graph_04')
```

![graph_04.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/ea1dd1a9-90e4-724e-1015-a26fbef606bd.png)


を考えます.

これまでのDAGとは異なり, XはMを介してYと因果関係があります.

このDAGに基いて, XとYをいずれもバイナリ変数とした以下のDGP

```math
\begin{align} \\
&Z \sim \mathcal{N}(0, 1) \\
&X = 1\{z + \epsilon_x > 0.8\}, \quad \epsilon_x \sim \mathcal{N}(0, 1)\\
&M = 0.5x + 1\{z > 0\} +  0.5 \epsilon_m, \quad \epsilon_m \sim \mathcal{N}(0, 1) \\
&Y = 1\{m >0\}\\
\end{align}
```

を仮定します.

DAGにおいて, Yの直接的な親はMで, そのMに対してはZとXが親になっています.

ZとXはいずれもMの数値にプラスの因果関係を持ちますが, 相対的な影響度合いはZの方が大きいことが見て取れます.

これまでと同様に, $y \sim x$, $y \sim x + z$で回帰した場合のそれぞれについて所定のメトリクスを出力すると,

```python
# Set random seed for reproductivity 
np.random.seed(45)

# assign values to the variables: 
#u = np.random.randn(N_SAMPLES)
z = np.random.randn(N_SAMPLES)
x = np.where(z + np.random.randn(N_SAMPLES) >0.8, 1, 0)
m = 0.5 * x + np.where(z>0, 1, 0) + 0.5 * np.random.randn(N_SAMPLES)
y = np.where(m>0, 1, 0)

def return_reg_summary(): 
  # Define four model variants 
  variants = [
      [x], 
      [x, z]
  ]
  
  # Fit models iteratively and store the results 
  results = []
  for variant in variants: 
      X = sm.add_constant(np.stack(variant).T)
  
      # Instantiate the model and fit it 
      model = sm.OLS(y, X)
      fitted_model = model.fit()
  
      results.append([fitted_model.params, fitted_model.pvalues])
  
      print(f'Params: {fitted_model.params}')
      print(f'p-vals: {fitted_model.pvalues}')
      print(f'Signif: {fitted_model.pvalues <= .05}\n')

return_reg_summary()

```
```
## Params: [0.66666667 0.31167268]
## p-vals: [3.65104979e-235 8.82797016e-026]
## Signif: [ True  True]
## 
## Params: [0.72181478 0.1305511  0.16434868]
## p-vals: [1.18031147e-261 3.19386870e-005 1.82325130e-029]
## Signif: [ True  True  True]
```

のようになります.

今回は, XとYの間には因果関係のパスが存在します.

従って, Xの偏回帰係数が2つのモデル共に有意となることとは, DAGと整合的です.

また, 2つ目のモデルにおけるXの偏回帰係数の推定値0.1305511が, 1つ目のモデルにおけるXの偏回帰係数の推定値0.31167268よりも小さな値であることは, Zが**Good control**であることを示しています.

今回のケースでは, Yがバイナリ変数であるなど, DGPを考慮すると, DGPで言及されている説明変数をそのまま重回帰式に投入することが適切ということではありませんが, シンプルな重回帰を行うことで, 交絡Zの影響を緩和した偏回帰係数を得ることができました.

### Model5; good control(blocking back-door paths)

次に, モデル4と同様に, **Mediator**(**Chain**)が存在するようなDAGで, 未観測の変数Uが存在するケース

```python
# Create a directed graph
g_5 = graphviz.Digraph(format='png')

## Add nodes
nodes_5 = ['Z', 'X', 'Y', 'M', 'U']
[g_5.node(n) for n in nodes_5] 

g_5.edges(['UZ', 'ZX', 'UM', 'XM', 'MY'])

# Render for print 
g_5.render('img/graph_05')
```

![graph_05.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/294196c6-735c-54fc-922b-31b2ec2b4247.png)


について考えます.

更に, 先ほどのDGPを少し手直しした,

```math
\begin{align} \\　
&U \sim \mathcal{U}(0, 1) \\
&Z = u + 0.1\epsilon_z, \quad \epsilon_z \sim \mathcal{N}(0, 1) \\
&X = 1\{z^2 + 0.1\epsilon_x > 0.5\}, \quad \epsilon_x \sim \mathcal{N}(0, 1) \\
&M = 1\{u+x > 0.5\} + \epsilon_m, \quad \epsilon_m \sim \mathcal{N}(0, 1) \\
&Y = 1\{m + \epsilon_y >0\}, \quad \epsilon_y \sim \mathcal{N}(0, 1)\\
\end{align}
```

のようなDGPを仮定します.

ここで, $\mathcal{U}(0, 1)$は区間(0, 1)の一様分布です.

全体としてX固有の影響というよりも, Xの親であるUの影響を大きくしてみました.

これまでと同様に, YをXで回帰した場合, YをX+Zで回帰した場合のそれぞれについて所定のメトリクスを出力すると,

```python
# Set random seed for reproductivity 
np.random.seed(45)

# assign values to the variables: 
u = np.random.uniform(0, 1, N_SAMPLES)
z =  u + 0.1 * np.random.randn(N_SAMPLES)
x = np.where(z**2 + 0.1 * np.random.randn(N_SAMPLES)>0.5, 1, 0)
m = np.where(u+x>0.5, 1, 0)+ np.random.randn(N_SAMPLES)
y = np.where(m + + np.random.randn(N_SAMPLES)>0, 1, 0)

return_reg_summary()
```
```
## Params: [0.58873239 0.20092278]
## p-vals: [7.65628930e-165 1.19305871e-009]
## Signif: [ True  True]
## 
## Params: [0.49955193 0.06525494 0.26018718]
## p-vals: [6.60600052e-50 2.08294106e-01 7.96183431e-04]
## Signif: [ True False  True]
```

2つ目のモデルにおけるXの偏回帰係数の推定値0.06525494が, 1つ目のモデルにおけるXの偏回帰係数の推定値0.20092278よりも小さな値であることは, DAG及びDGPと整合的です.

一方で, 2つ目のモデルにおいて, Xの偏回帰係数の推定値のP-値は2.08294106e-01となり, 有意水準を5%とした場合に有意であるとは言えない結果になりました.

このことは, DAG及びDGPと整合的であるとはいえません.

つまり, シンプルな重回帰モデルにおける偏回帰係数の推定値のP-値は, DGPの内容(または有意水準の設定)次第で, DAGと不整合となるケースが存在することが分かりました.

今回のケースでは, Yがバイナリ変数であるなど, DGPを考慮すると, DGPで言及されている説明変数をそのまま重回帰式に投入することが適切ということではありませんが, シンプルな重回帰を行うことで, 交絡Zの影響を緩和した偏回帰係数を得ることができました.

### Model.6; good control(blocking back-door paths)

モデル5の変形として下図のようなDAG, DGPを考えます.

```python
# Create a directed graph
g_6 = graphviz.Digraph(format='png')

## Add nodes
nodes_6 = ['Z', 'X', 'Y', 'M', 'U']
[g_6.node(n) for n in nodes_5] 

g_6.edges(['UX', 'UZ', 'ZM', 'XM', 'MY'])

# Render for print 
g_6.render('img/graph_06')
```

![graph_06.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/2deb316c-424e-a83f-280c-a992a83c84ae.png)


```math
\begin{align} \\　
&U \sim \mathcal{U}(0, 1) \\
&Z = u + 0.1\epsilon_z, \quad \epsilon_z \sim \mathcal{N}(0, 1) \\
&X = 1\{u^2 + 0.1\epsilon_x > 0.5\}, \quad \epsilon_x \sim \mathcal{N}(0, 1) \\
&M = 1\{z+x > 0.5\} + \epsilon_m, \quad \epsilon_m \sim \mathcal{N}(0, 1) \\
&Y = 1\{m + \epsilon_y >0\}, \quad \epsilon_y \sim \mathcal{N}(0, 1)\\
\end{align}
```

このケースでもZは**Good control**となります.

モデル5と同様に, Yに与える影響の大きさは, Xよりも, Zの方が大きくしています.

YをXで回帰した場合, YをX+Zで回帰した場合のそれぞれについて所定のメトリクスを出力すると,

```python
# Set random seed for reproductivity 
np.random.seed(45)

# assign values to the variables: 
u = np.random.uniform(0, 1, N_SAMPLES)
z =  u + 0.1 * np.random.randn(N_SAMPLES)
x = np.where(u**2 + 0.1 * np.random.randn(N_SAMPLES)>0.5, 1, 0)
m = np.where(z+x>0.5, 1, 0)+ np.random.randn(N_SAMPLES)
y = np.where(m + np.random.randn(N_SAMPLES)>0, 1, 0)

return_reg_summary()

```
```
## Params: [0.58707865 0.21500468]
## p-vals: [1.93437317e-165 7.45351535e-011]
## Signif: [ True  True]
## 
## Params: [0.50370716 0.09836625 0.23678328]
## p-vals: [3.19371106e-54 3.99274283e-02 9.14759834e-04]
## Signif: [ True  True  True]
```


のようになります.

今回はモデル5とは異なり, 2つのモデルにおけるXの偏回帰係数の推定値は共に有意となり, DAG, DGPと整合的な結果となりました.

また, 2つ目のモデルにおけるXの偏回帰係数の推定値0.09836625が, 1つ目のモデルにおけるXの偏回帰係数の推定値0.21500468よりも小さな値であることは, DAG及びDGPと整合的です.

モデル5で設定したDGPと同様, Yがバイナリ変数であるなど, DGPで言及されている説明変数をそのまま重回帰式に投入することが適切ということではありませんが, シンプルな重回帰を行うことで, 交絡Zの影響を緩和されたことが分かるような偏回帰係数を得ることができました.

### Good controlsのまとめ 

このセクションで説明したDAG

![graph_01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/3cfbaf25-f37e-6327-8174-6ca54b9f3511.png)

![graph_02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/ec6c8bf4-c64b-5af7-60d4-62aa4cc2860b.png)

![graph_03.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/a50dd7e5-7376-500e-7be2-03100ad491ce.png)
![graph_05.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/18346727-f94a-e6f6-3d62-8e8be71f1956.png)
![graph_06.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/4275cacc-a2ab-8858-2053-75eca5acd367.png)
![graph_04.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/504e9917-f05f-04fa-5d26-4845d03f5383.png)




のいずれにおいても, Zが**Good control**であることが確認できました.

これらのZは**common cause**又は**fork**とも呼ばれるようです。

## Neutral controls

### Model.8; Neutral control

下図のようなDAG, DGPを考えます.

```python
# Create a directed graph
g_8 = graphviz.Digraph(format='png')

## Add nodes
nodes_8 = ['Z', 'X', 'Y']
[g_8.node(n) for n in nodes_8] 

g_8.edges(['XY', 'ZY'])

# Render for print 
g_8.render('img/graph_08')
```

![graph_08.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/7aa9dacc-c797-1da0-6ec1-b092fa5c01da.png)


```math
\begin{align} \\
&Z \sim \mathcal{N}(0, 1) \\
&X \sim \mathcal{U}(0, 1) \\
&Y = 10z - x \\
\end{align}
```

このケースでZは漸近的には**Neutral control**とされます.

ただし, 標本サイズが有限なレジームにおいては, 上図のZを説明変数に追加することにより, XがYに与える影響の推定精度が**向上**する可能性があります.

YをXで回帰した場合, YをX+Zで回帰した場合のそれぞれについて所定のメトリクスを出力すると,

```python
# Set random seed for reproductivity 
np.random.seed(45)

# assign values to the variables: 
z = np.random.randn(N_SAMPLES)
x = np.random.uniform(0, 1, N_SAMPLES)
y = 10 * z - x

return_reg_summary()

```
```
## Params: [ 0.72089323 -3.02869546]
## p-vals: [0.25268028 0.00526912]
## Signif: [False  True]
## 
## Params: [-7.09501902e-16 -1.00000000e+00  1.00000000e+01]
## p-vals: [1.02083218e-07 0.00000000e+00 0.00000000e+00]
## Signif: [ True  True  True]
```


2つのモデルにおけるXの偏回帰係数の推定値は共に有意となり, DAG, DGPと整合的な結果となりました.

一方で, 2つ目のモデルにおけるXの偏回帰係数の推定値-1.00000000e+00が, 1つ目のモデルにおけるXの偏回帰係数の推定値-3.02869546よりも大きな値となり, 約3倍変動していることが分かります.

漸近的に**Neutral control**であるとのことでしたので, 標本サイズを1,000ではなく1,000,000に変更してみます.

```python
# Set random seed for reproductivity 
np.random.seed(45)

# temporarily modify the number of samples: 
tmp_n_samples = 1000000

# assign values to the variables: 
## z = np.random.randn(N_SAMPLES)
## x = np.random.uniform(0, 1, N_SAMPLES)
z = np.random.randn(tmp_n_samples)
x = np.random.uniform(0, 1, tmp_n_samples)
y = 10 * z - x

return_reg_summary()

```
```
## Params: [-0.00827811 -0.99907185]
## p-vals: [6.78722182e-001 6.86106700e-183]
## Signif: [False  True]
## 
## Params: [-4.14555543e-15 -1.00000000e+00  1.00000000e+01]
## p-vals: [1.22062246e-16 0.00000000e+00 0.00000000e+00]
## Signif: [ True  True  True]
```


今度はZが漸近的に**Neutral control**であることと整合的な結果となりました.

このような**Neutral control**を説明変数から除外しないように注意する必要があります.

**Neutral control**であるようなZを説明変数から除外すると, 標本サイズ次第では, Xの偏回帰係数の推定値が真の値から少なくない乖離をする可能性について考慮する必要があると言えるでしょう.

### Model.9; Neutral control

下図のようなDAG, DGPを考えます.

```python
# Create a directed graph
g_9  = graphviz.Digraph(format='png')

## Add nodes
nodes_9 = ['Z', 'X', 'Y']
[g_9.node(n) for n in nodes_9] 

g_9.edges(['ZX', 'XY'])

# Render for print 
g_9.render('img/graph_09')
```

![graph_09.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/fcfcaa5a-5ee4-28bb-5d59-3fcfddaa9ef0.png)


```math
\begin{align} \\
&Z \sim \mathcal{N}(0, 1) \\
&X := 1\{z \cdot \epsilon_x < -1.96 \}, \quad \epsilon_x \sim \mathcal{U}(0, 1) \\
&Y := x +  \epsilon_y, \quad \epsilon_y \sim \mathcal{N}(0, 1) \\
\end{align}
```

モデル8と同様に, このケースでZは漸近的には**Neutral control**とされます.

一方で, 標本サイズが有限なレジームにおいては, モデル8とは異なり,　上図のZを説明変数に追加することにより, XがYに与える影響の推定精度が**悪化**する可能性があります.

YをXで回帰した場合, YをX+Zで回帰した場合のそれぞれについて所定のメトリクスを出力すると,

```python
# Set random seed for reproductivity 
np.random.seed(45)

# assign values to the variables: 
z = np.random.randn(N_SAMPLES)
x = np.where(
  z * np.random.uniform(0, 1, N_SAMPLES) < -1.96, 
  1, 
  0)
y = x + np.random.randn(N_SAMPLES)

return_reg_summary()

```
```
## Params: [-0.04003066  1.69575722]
## p-vals: [0.20687297 0.00074249]
## Signif: [False  True]
## 
## Params: [-0.039679    1.74396334  0.01797812]
## p-vals: [0.21113638 0.00063357 0.57933439]
## Signif: [False  True False]
```


2つのモデルにおけるXの偏回帰係数の推定値は共に有意となり, DAG, DGPと整合的な結果となりました.

一方で, 2つ目のモデルにおけるXの偏回帰係数の推定値1.74396334が, 1つ目のモデルにおけるXの偏回帰係数の推定値1.69575722よりも大きな値となりました.

ちなみに, 標本サイズを1,000ではなく1,000,000に変更すると, 結果は以下のようになります.

```python
# Set random seed for reproductivity 
np.random.seed(45)

# temporarily modify the number of samples: 
tmp_n_samples = 1000000

# assign values to the variables: 
z = np.random.randn(tmp_n_samples)
x = np.where(
  z * np.random.uniform(0, 1, tmp_n_samples) < -1.96, 
  1, 
  0)
y = x + np.random.randn(tmp_n_samples)

return_reg_summary()

```
```
## Params: [-6.14266061e-04  1.00326741e+00]
## p-vals: [0.53997471 0.        ]
## Signif: [False  True]
## 
## Params: [-6.06301700e-04  1.00087441e+00 -9.24223843e-04]
## p-vals: [0.54525726 0.         0.3616708 ]
## Signif: [False  True False]
```


上図のZが漸近的に**Neutral control**であることと整合的な結果となりました.

### Model.13; Neutral control

下図のようなDAG, DGPを考えます.

```python
# Create a directed graph
g_13  = graphviz.Digraph(format='png')

## Add nodes
nodes_13 = ['Z', 'X', 'M', 'Y']
[g_13.node(n) for n in nodes_13] 

g_13.edges(['XM', 'MY', 'ZM'])

# Render for print 
g_13.render('img/graph_13')
```

![graph_13.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/ab3ab39b-3836-c1ad-7a60-b1b6ac5f6d98.png)


このケースにおけるZはモデル8, モデル9と同様に,漸近的には**Neutral control**とされます.

さらに, モデル8と同様, 標本サイズが有限なレジームにおいては, 上図のZを説明変数に追加することにより, XがYに与える影響の推定精度が**向上**する可能性があります.

設例として, 次のようなDGPを考えます.

```math
\begin{align} \\
&Z := 1 + 0.01\epsilon_z, \quad \epsilon_z \sim  \mathcal{N}(0, 1) \\
&X \sim \mathcal{U}(0, 1) \\
& M:= 1.5z - x + 0.01\epsilon_m, \quad \epsilon_m \sim \mathcal{N}(0, 1)\\
&Y := 1\{ m > 1 \} \\
\end{align}
```

$y \sim x$, $y \sim x + z$のそれぞれで回帰した場合の所定のメトリクスを出力してみましょう.

```python
# Set random seed for reproductivity 
np.random.seed(45)

# assign values to the variables: 
z = 1 + 0.01 * np.random.randn(N_SAMPLES)
x = np.random.uniform(0, 1, N_SAMPLES)
m = 1.5 * z - x + 0.01 * np.random.randn(N_SAMPLES)
y = np.where(m > 1, 1, 0)

return_reg_summary()

```
```
## Params: [ 1.24618095 -1.49054326]
## p-vals: [0.00000000e+000 4.36423264e-296]
## Signif: [ True  True]
## 
## Params: [-0.63962897 -1.48672028  1.88445143]
## p-vals: [4.30757031e-001 2.85316663e-295 2.03068907e-002]
## Signif: [False  True  True]
```



2つのモデルにおけるXの偏回帰係数の推定値は共に有意となり, DAG, DGPと整合的な結果となりました.

次に, 2つ目のモデルにおけるXの偏回帰係数の推定値-1.48672028が, 1つ目のモデルにおけるXの偏回帰係数の推定値-1.49054326よりも大きな値となりました.

今回のDAGは, 標本サイズが有限なレジームにおいては, Zを説明変数に追加することにより, XがYに与える影響の推定精度が**向上**することが期待されるとのことでした.

この主張が正しければ, 2つ目のモデルにおけるXの偏回帰係数の推定値-1.48672028を算出するために用いた推定量は, 1つ目のモデルにおけるXの偏回帰係数の推定値-1.4905432を算出するために用いた推定量よりも優れているということになります.

数値の変動幅を考慮すると私なら誤差による変動だと考えますが果たして本当なのでしょうか?

このことを確認するために, 標本サイズを1,000ではなく1,000,000に変更すると, 結果は以下のようになります.

```python
# Set random seed for reproductivity 
np.random.seed(45)

# temporarily modify the number of samples: 
tmp_n_samples = 1000000

# assign values to the variables: 
z = 1 + 0.01 * np.random.randn(tmp_n_samples)
x = np.random.uniform(0, 1, tmp_n_samples)
m = 1.5 * z - x + 0.01 * np.random.randn(tmp_n_samples)
y = np.where(m > 1, 1, 0) 

return_reg_summary()

```
```
## Params: [ 1.24886375 -1.49801363]
## p-vals: [0. 0.]
## Signif: [ True  True]
## 
## Params: [-0.2759154  -1.49801505  1.52479177]
## p-vals: [3.75921905e-28 0.00000000e+00 0.00000000e+00]
## Signif: [ True  True  True]
```


2つのモデルにおけるXの偏回帰係数の推定値は約-1.49801となりました.

ここで, 先ほどの重回帰結果において, 2つ目のモデルにおけるXの偏回帰係数の推定値-1.48672028が, 1つ目のモデルにおけるXの偏回帰係数の推定値-1.49054326よりも大きな値となったことに留意すると, この設例では, Zを説明変数に追加することにより, XがYに与える影響の推定精度が**向上**した訳ではないようです.　

このDGPにシンプルな重回帰をあてはめ, 連続値であるXの偏回帰係数を分析するのには無理があるとも言えるので, 今回限定的な標本サイズ下で推定精度が特段向上しなかったのは, 何か別の原因があるかもしれません.

### Model.14; Neutral control

下図のようなDAG, DGPを考えます.

```python
# Create a directed graph
g_14  = graphviz.Digraph(format='png')

## Add nodes
nodes_14 = ['Z', 'X', 'Y']
[g_14.node(n) for n in nodes_14] 

g_14.edges(['XY', 'XZ'])

# Render for print 
g_14.render('img/graph_14')
```

![graph_14.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/8d91986a-5daa-2094-e4cf-109d651e7462.png)


このケースにおけるZは**Neutral control**とされます.

なお, このZはXの**Post treatment variable**と呼ばれます.

さらに, 標本サイズが有限なレジームにおいて, 上図のZをコントロールすることでXのバリエーションが制限されるような状況では, XがYに与える影響の推定精度が**悪化**する可能性があります.

設例として, 次のようなDGPを考えます.

```math
\begin{align} \\
&X \sim \mathcal{U}(0, 1) \\
&Z := 1\{ x + 0.1 \epsilon_z < 0.5\}, \quad \epsilon_z \sim \mathcal{N}(0, 1)\\
&Y := 2.5x +  \epsilon_y, \quad \epsilon_y \sim \mathcal{U}(0, 1)\\
\end{align}
```


ご覧のとおり, ZでコントロールするとXのバリエーションが小さくなります.

$y \sim x$, $y \sim x + z$のそれぞれで回帰した場合の所定のメトリクスを出力してみましょう.

```python
# Set random seed for reproductivity 
np.random.seed(45)

# assign values to the variables: 
x = np.random.uniform(0, 1, N_SAMPLES)
z = np.where(
  x + 0.1 * np.random.randn(N_SAMPLES) < 0.5, 
  1, 
  0
)
y = x * 2.5 +  np.random.uniform(0, 1, N_SAMPLES)

return_reg_summary()

```
```
## Params: [0.47759209 2.53056834]
## p-vals: [1.33984614e-120 0.00000000e+000]
## Signif: [ True  True]
## 
## Params: [ 0.5788626   2.41212051 -0.08303408]
## p-vals: [1.35118939e-037 1.20383072e-231 1.05789523e-002]
## Signif: [ True  True  True]
```


2つのモデルにおけるXの偏回帰係数の推定値は共に有意となり, DAG, DGPと整合的な結果となりました.

次に, 2つ目のモデルにおけるXの偏回帰係数の推定値2.41212051が, 1つ目のモデルにおけるXの偏回帰係数の推定値2.53056834よりも小さな値となりました.

DAGによる真のパラメータは今回の場合2.5であることが明快であるため, ZをコントロールすることでXのバリエーションが制限されるような状況では, XがYに与える影響の推定精度が**悪化**する可能性があることをこの設例で示すことができました.

ただし, 数値の変動幅を考慮すると, 誤差による変動だと考えることも可能かと思います.

なお, 2つ目のモデルにおけるZの偏回帰係数の推定値が有意となりました.

DAGからすると, Zの偏回帰係数の推定値が有意であって欲しくないものです.

ですので, 標本サイズを1,000ではなく1,000,000に変した結果も確認してみましょう.

```python
# Set random seed for reproductivity 
np.random.seed(45)

# temporarily modify the number of samples: 
tmp_n_samples = 1000000

# assign values to the variables: 
x = np.random.uniform(0, 1, tmp_n_samples)
z = np.where(
  x + 0.1 * np.random.randn(tmp_n_samples) < 0.5, 
  1, 
  0
)
y = x * 2.5 +  np.random.uniform(0, 1, tmp_n_samples)

return_reg_summary()

```
```
## Params: [0.49992247 2.50047852]
## p-vals: [0. 0.]
## Signif: [ True  True]
## 
## Params: [ 5.01436658e-01  2.49869225e+00 -1.24122622e-03]
## p-vals: [0.         0.         0.23236849]
## Signif: [ True  True False]
```

今度は, (2つ目のモデルにおける)Zの偏回帰係数の推定値が有意となりませんでした.

なお, 設例のDGPでは$Y := 2.5x + \epsilon_y, \quad \epsilon_y \sim \mathcal{U}(0, 1)$のように, 目的変数Yの誤差項が一様分布に従うこととしました.

この誤差項を, 一様分布ではなく, 標準正規分布に従うようなDGPを仮定した場合には, N=1,000の標本サイズにおいて,

```python
# Set random seed for reproductivity 
np.random.seed(45)

# assign values to the variables: 
x = np.random.uniform(0, 1, N_SAMPLES)
z = np.where(
  x + 0.1 * np.random.randn(N_SAMPLES) < 0.5, 
  1, 
  0
)
y = x * 2.5 +  np.random.randn(N_SAMPLES)

return_reg_summary()

```
```
## Params: [0.03644642 2.35250388]
## p-vals: [5.56958548e-01 1.07556403e-85]
## Signif: [False  True]
## 
## Params: [ 0.16050041  2.20740808 -0.1017148 ]
## p-vals: [2.91281142e-01 5.72853550e-28 3.71569122e-01]
## Signif: [False  True False]
```


のような結果となり, (2つ目のモデルにおける)Zの偏回帰係数の推定値が有意とはいえなくなるようです.

### Model.15; Neutral control

下図のようなDAGを考えます.

```python
# Create a directed graph
g_15  = graphviz.Digraph(format='png')

## Add nodes
nodes_15 = ['Z', 'X', 'Y', 'U', 'W']
[g_15.node(n) for n in nodes_15] 

g_15.edges(['XY', 'XZ', 'ZW', 'UW', 'UY'])

# Render for print 
g_15.render('img/graph_15')
```

![graph_15.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/c67440ac-f9ac-07fb-0f01-eb5b17c03d25.png)


モデル14と同様に, このZはXの**Post treatment variable**と呼ばれます.

また, このケースにおけるZは**Neutral control**とされます. さらに, 標本サイズが有限なレジームにおいて, 上図のZをコントロールすることでXのバリエーションが制限されるような状況では, XがYに与える影響の推定精度が**悪化**する可能性があります.

設例として, 次のようなDGPを考えます.

```math
\begin{align} \\
&X \sim \mathcal{U}(0, 1) \\
&U \sim \mathcal{N}(0, 1) \\
&Z := 1\{ x + 0.1 \epsilon_z < 0.2\}, \quad \epsilon_z \sim \mathcal{N}(0, 1)\\
&Y := 2.5x - 2.5u\\ 
&W := u + z\\
\end{align}
```

ご覧のとおり, ZでコントロールするとXのバリエーションが小さくなります.

$y \sim x$, $y \sim x + z$のそれぞれで回帰した場合の所定のメトリクスを出力してみましょう.

```python
# Set random seed for reproductivity 
np.random.seed(45)

# assign values to the variables: 
x = np.random.uniform(0, 1, N_SAMPLES)
u = np.random.randn(N_SAMPLES)
z = np.where(
  x + 0.1 * np.random.randn(N_SAMPLES) < 0.2, 
  1, 
  0
)
y = x * 2.5 - u * 2.5
w = u + z
return_reg_summary()

```
```
## Params: [-0.22004632  2.83985173]
## p-vals: [1.56502049e-01 2.27213226e-24]
## Signif: [False  True]
## 
## Params: [-0.14263676  2.73243508 -0.11703056]
## p-vals: [5.35931034e-01 7.23474161e-14 6.49325703e-01]
## Signif: [False  True False]
```


2つのモデルにおけるXの偏回帰係数の推定値は共に有意となり, DAG, DGPと整合的な結果となりました.

次に, 2つ目のモデルにおけるXの偏回帰係数の推定値2.73243508が, 1つ目のモデルにおけるXの偏回帰係数の推定値2.83985173よりも小さな値となりました.

DAGによる真のパラメータは今回の場合2.5であることが明快であるため, ZをコントロールすることでXのバリエーションが制限されるような状況では, XがYに与える影響の推定精度が**悪化**する可能性があることはこの設例では示されませんでした.

合成データですので, 繰り返し回帰を行い推定値をヒストグラムにすることで, さらに掘り下げて確認することが可能と考えますが, ここでは割愛します.

念のため, 標本サイズを1,000ではなく1,000,000に変した結果も確認してみましょう.

```python
# Set random seed for reproductivity 
np.random.seed(45)

# temporarily modify the number of samples: 
tmp_n_samples = 1000000

# assign values to the variables: 
x = np.random.uniform(0, 1, tmp_n_samples)
u = np.random.randn(tmp_n_samples)
z = np.where(
  x + 0.1 * np.random.randn(tmp_n_samples) < 0.2, 
  1, 
  0
)
y = x * 2.5 - u * 2.5
w = u + z
return_reg_summary()

```
```
## Params: [-1.48565976e-03  2.50918182e+00]
## p-vals: [0.76649497 0.        ]
## Signif: [False  True]
## 
## Params: [-1.87736984e-03  2.50972410e+00  6.00379224e-04]
## p-vals: [0.79820641 0.         0.94190461]
## Signif: [False  True False]
```


上図のZが漸近的に**Neutral control**であることと整合的な結果となりました.

### Neutral controlのまとめ

このセクションで説明したDAG

![graph_08.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/0db0001e-47b1-c1d5-200f-40e21d169858.png)
![graph_09.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/4ee1a4b6-274e-571b-aedf-5f297a9ba6b0.png)
![graph_13.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/18e05236-8ab8-e476-e68a-2e6415186207.png)
![graph_14.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/67a652f1-9e0c-8ad4-4009-a96a4ec6f8ad.png)
![graph_15.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/80936808-71c7-891e-62a4-447d9a5d5673.png)


のいずれにおいても, Zが漸近的に**Neutral control**であることが確認できました.

一方で, 標本サイズを1,000とした場合, Zを説明変数に追加したときのXの偏回帰係数の推定値の変動幅・方向はモデルやDGPによってまちまちでした.

モデル8, モデル15のDGPでは,DGPが表す真の値に近づきましたが, モデル9, モデル14のDGPでは逆に, DGPが表す真の値から遠ざかりました.

さらに, モデル13のDGPではZを説明変数に追加してものXの偏回帰係数の推定値は大きく変動しませんでした

Zを説明変数に追加したときのXの偏回帰係数の推定値がDGPが表す真の値から遠ざかる理由として, ZがXのバリエーションを制限する場合の推定精度悪化について指摘しました.

しかし, ZがXのバリエーションを制限するようなケースでも, モデル15のように推定精度が結果として向上したと解釈できる事例も存在しました.

## Bad controls

### Model.7; Bad control(M-bias)

下図のようなDAGを考えます.

```python
# Create a directed graph
g_7  = graphviz.Digraph(format='png')

## Add nodes
nodes_7 = ['Z', 'X', 'Y', 'U1', 'U2']
[g_7.node(n) for n in nodes_7] 

g_7.edges(['XY', ('U1', 'X'), ('U2', 'Y'), ('U1', 'Z'), ('U2', 'Z')])

# Render for print 
g_7.render('img/graph_7')
```

![graph_7.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/b15d1086-4069-eb7b-e8c3-00cab1abd290.png)


このケースでZは**Bad control**とされ, このDAGにおけるZは**M-bias**とも呼ばれます.

設例として, 次のようなDGPを考えます.

```math
\begin{align} \\
&U_1 \sim \mathcal{U}(0, 1) \\
&U_2 \sim \mathcal{U}(0, 1) \\
&X := u_1 + \epsilon_x, \quad \epsilon_x \sim \mathcal{N}(0, 1) \\
&Z := 1\{u_1 + u_2 >0.8\} \\
&Y := 2u_2 + 1.5x \\
\end{align}
```

$y \sim x$, $y \sim x + z$のそれぞれで回帰した場合の所定のメトリクスを出力してみましょう.

```python
# Set random seed for reproductivity 
np.random.seed(45)

# assign values to the variables: 
u_1 = np.random.uniform(0, 1, N_SAMPLES)
u_2 = np.random.uniform(0, 1, N_SAMPLES)
x = u_1 + np.random.randn(N_SAMPLES)
z = np.where(u_1 + u_2 > 0.8, 1, 0)
y = 2 * u_2 + 1.5 * x

return_reg_summary()

```
```
## Params: [1.00322495 1.50299028]
## p-vals: [8.8844752e-277 0.0000000e+000]
## Signif: [ True  True]
## 
## Params: [0.52804617 1.45224423 0.73056826]
## p-vals: [1.94707838e-74 0.00000000e+00 1.43213979e-91]
## Signif: [ True  True  True]
```


2つのモデルにおけるXの偏回帰係数の推定値は共に有意となり, DAG, DGPと整合的な結果となりました.

次に, 2つ目のモデルにおけるXの偏回帰係数の推定値1.45224423が, 1つ目のモデルにおけるXの偏回帰係数の推定値1.50299028よりも小さな値となりました.

DAGによる真のパラメータは今回の場合1.5であることが明快であるため, Zを説明変数に追加したときのXの偏回帰係数の推定値の精度が**悪化**したと解釈することが可能です.

ただし, 数値の変動幅を考慮すると, 誤差による変動だと考えることもできます.

ですので, 標本サイズを1,000ではなく1,000,000にした時の結果も確認してみましょう.

```python
# Set random seed for reproductivity 
np.random.seed(45)

# temporarily modify the number of samples: 
tmp_n_samples = 1000000

# assign values to the variables: 
u_1 = np.random.uniform(0, 1, tmp_n_samples)
u_2 = np.random.uniform(0, 1, tmp_n_samples)
x = u_1 + np.random.randn(tmp_n_samples)
z = np.where(u_1 + u_2 > 0.8, 1, 0)
y = 2 * u_2 + 1.5 * x

return_reg_summary()
```
```
## Params: [0.99829881 1.50091951]
## p-vals: [0. 0.]
## Signif: [ True  True]
## 
## Params: [0.54458716 1.4524316  0.70253239]
## p-vals: [0. 0. 0.]
## Signif: [ True  True  True]
```


標本サイズを増加させても, Zを説明変数に追加したときのXの偏回帰係数の推定値の精度が**悪化**する事象は緩和されていないことが分かります.

これは, 漸近的には**Neutral**な変数とは異なる特徴といえるでしょう.

以上が, 文献で紹介されているDAGに関する考察でした.

ここで, 先ほどのDAGに修正を加え, XとYの因果パスを除去したDAG及びDGP

```python
# Create a directed graph
g_7_modified  = graphviz.Digraph(format='png')

## Add nodes
nodes_7_modified = ['Z', 'X', 'Y', 'U1', 'U2']
[g_7_modified.node(n) for n in nodes_7_modified] 

g_7_modified.edges([('U1', 'X'), ('U2', 'Y'), ('U1', 'Z'), ('U2', 'Z')])

# Render for print 
g_7_modified.render('img/graph_7_modified')
```

![graph_7_modified.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/8bad1af5-c755-6014-da2e-d40c4b806106.png)


```math
\begin{align} \\
&U_1 \sim \mathcal{U}(0, 1) \\
&U_2 \sim \mathcal{U}(0, 1) \\
&X := u_1 + \epsilon_x, \quad \epsilon_x \sim \mathcal{N}(0, 1) \\
&Z := 1\{u_1 + u_2 >0.8\} \\
&Y := 2u_2 \\
\end{align}
```


について, $y \sim x$, $y \sim x + z$のそれぞれで回帰した場合の所定のメトリクスを確認します.

```python
# Set random seed for reproductivity 
np.random.seed(45)

# assign values to the variables: 
u_1 = np.random.uniform(0, 1, N_SAMPLES)
u_2 = np.random.uniform(0, 1, N_SAMPLES)
x = u_1 + np.random.randn(N_SAMPLES)
z = np.where(u_1 + u_2 > 0.8, 1, 0)
y = 2 * u_2

return_reg_summary()

```
```
## Params: [1.00322495 0.00299028]
## p-vals: [8.8844752e-277 8.6245552e-001]
## Signif: [ True False]
## 
## Params: [ 0.52804617 -0.04775577  0.73056826]
## p-vals: [1.94707838e-74 8.13409786e-04 1.43213979e-91]
## Signif: [ True  True  True]
```


1つ目のモデルにおけるXの偏回帰係数の推定値は有意とはいえず, DAG, DGPと整合的な結果となりました.

一方で, 2つ目のモデルにおけるXの偏回帰係数の推定値は有意となり, DAG, DGPと平仄がとれない結果となりました.

**M-bias**を生起するZを説明変数に追加すると, Xの偏回帰係数の推定値が有意となることが分かりました.

標本サイズを1,000ではなく1,000,000にした時の結果も確認してみましょう.

```python
# Set random seed for reproductivity 
np.random.seed(45)

# temporarily modify the number of samples: 
tmp_n_samples = 1000000

# assign values to the variables: 
u_1 = np.random.uniform(0, 1, tmp_n_samples)
u_2 = np.random.uniform(0, 1, tmp_n_samples)
x = u_1 + np.random.randn(tmp_n_samples)
z = np.where(u_1 + u_2 > 0.8, 1, 0)
y = 2 * u_2

return_reg_summary()
```
```
## Params: [9.98298809e-01 9.19513511e-04]
## p-vals: [0.         0.09743119]
## Signif: [ True False]
## 
## Params: [ 0.54458716 -0.0475684   0.70253239]
## p-vals: [0. 0. 0.]
## Signif: [ True  True  True]
```


先ほどと同様に, 標本サイズを大きくしたからといってバイアスが削減される訳ではないようです.

**M-bias**を生起するZを説明変数に追加すると, 標本サイズを大きくしても, 依然としてXの偏回帰係数の推定値が有意となることが分かりました.

### Model.10; Bad control(bias amplification)

下図のようなDAGを考えます.

```python
# Create a directed graph
g_10  = graphviz.Digraph(format='png')

## Add nodes
nodes_10 = ['Z', 'X', 'Y', 'U']
[g_10.node(n) for n in nodes_10] 

g_10.edges(['ZX', 'UX', 'UY', 'XY'])

# Render for print 
g_10.render('img/graph_10')
```

![graph_10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/8901b6a7-1ca1-413a-1334-087cf17d9612.png)


このケースでZは**Bad control**とされ, このDAGにおけるZは**pre-treatment**に該当します.

さらに言えば, 交絡として除去すべきUが未観測であるような事例であると言えます.

設例として, 次のようなDGPを考えます.

```math
\begin{align} \\
&U \sim \mathcal{U}(0, 1) \\
&Z \sim \mathcal{N}(0, 1) \\
&X := 1\{u + z + \epsilon_x > 1.5\}, \quad \epsilon_x \sim \mathcal{N}(0, 1) \\
&Y := 5x -15u\\
\end{align}
```

$y \sim x$, $y \sim x + z$のそれぞれで回帰した場合の所定のメトリクスを出力してみましょう.

```python
# Set random seed for reproductivity 
np.random.seed(45)

# assign values to the variables: 
u = np.random.uniform(0, 1, N_SAMPLES)
z = np.random.randn(N_SAMPLES)
x = np.where(
  u + z + np.random.randn(N_SAMPLES) > 1.5, 
  1, 
  0)
y = 5 * x - 15 * u

return_reg_summary()

```
```
## Params: [-7.1089956   3.77264525]
## p-vals: [3.05230057e-247 3.21380023e-028]
## Signif: [ True  True]
## 
## Params: [-6.98433512  3.15143342  0.53814675]
## p-vals: [6.11849114e-236 1.71689503e-016 5.60556431e-004]
## Signif: [ True  True  True]
```

2つのモデルにおけるXの偏回帰係数の推定値は共に有意となり, DAG, DGPと整合的な結果となりました.

次に, 2つ目のモデルにおけるXの偏回帰係数の推定値3.15143342が, 1つ目のモデルにおけるXの偏回帰係数の推定値3.77264525よりも小さな値となりました.

DAGによる真のパラメータは今回の場合5であることが明快であるため, Zを説明変数に追加したときのXの偏回帰係数の推定値の精度が**悪化**したと解釈することが可能です.

ただし, 数値の変動幅を考慮すると, 誤差による変動だと考えることもできます.

ですので, 標本サイズを1,000ではなく1,000,000にした時の結果も確認してみましょう.

```python
# Set random seed for reproductivity 
np.random.seed(45)

# temporarily modify the number of samples: 
tmp_n_samples = 1000000

# assign values to the variables: 
u = np.random.uniform(0, 1, tmp_n_samples)
z = np.random.randn(tmp_n_samples)
x = np.where(
  u + z + np.random.randn(tmp_n_samples) > 1.5, 
  1, 
  0)
y = 5 * x - 15 * u

return_reg_summary()

```
```
## Params: [-7.15498711  3.54564107]
## p-vals: [0. 0.]
## Signif: [ True  True]
## 
## Params: [-7.03104009  3.0381356   0.43030811]
## p-vals: [0. 0. 0.]
## Signif: [ True  True  True]
```


標本サイズを増加させても, バイアスが緩和されていないことが分かります.

### Model.11; Bad control(overcontrol bias) 

下図のようなDAGを考えます.

```python
# Create a directed graph
g_11  = graphviz.Digraph(format='png')

## Add nodes
nodes_11 = ['Z', 'X', 'Y']
[g_11.node(n) for n in nodes_11] 

g_11.edges(['XZ', 'ZY'])

# Render for print 
g_11.render('img/graph_11')
```

![graph_11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/b719b3ea-7529-a9e0-b4bf-32a0a39cb206.png)


このケースでZは**Bad control**とされ, このDAGにおけるZは**Mediator**に該当します.

設例として, 次のようなDGPを考えます.

```math
\begin{align} \\
&X \sim \mathcal{U}(0, 1) \\
&Z := 1\{x+2\epsilon_x>0.5\} , \quad \epsilon_x \sim \mathcal{N}(0, 1) \\
&Y := z +\epsilon_y, \quad \epsilon_y \sim \mathcal{U}(0, 1)\\
\end{align}
```

$y \sim x$, $y \sim x + z$のそれぞれで回帰した場合の所定のメトリクスを出力してみましょう.

```python
# Set random seed for reproductivity 
np.random.seed(45)

# assign values to the variables: 
x = np.random.uniform(0, 1, N_SAMPLES)
z = np.where(
  x + 2 * np.random.randn(N_SAMPLES) > 0.5, 
  1, 
  0)
y = z + np.random.uniform(0, 1, N_SAMPLES)

return_reg_summary()

```
```
## Params: [0.89940823 0.20985761]
## p-vals: [1.06178607e-108 8.14551297e-004]
## Signif: [ True  True]
## 
## Params: [0.47595756 0.02987359 1.003875  ]
## p-vals: [4.74807167e-105 3.38502125e-001 1.12429053e-305]
## Signif: [ True False  True]
```


1つ目のモデルにおけるXの偏回帰係数の推定値は有意となり, DAG, DGPと整合的な結果となりました.

一方で, 2つ目のモデルにおけるXの偏回帰係数の推定値は有意ではなくなり, DAG, DGPと平仄がとれない結果となりました.

### Model.12; Bad control(overcontrol bias) 

モデル11では, Zが**mediator**となるDAGを想定しましたが, ここではZが**mediator**の直接の子となるようなDAG

```python
# Create a directed graph
g_12  = graphviz.Digraph(format='png')

## Add nodes
nodes_12 = ['Z', 'X', 'Y', 'M']
[g_12.node(n) for n in nodes_12] 

g_12.edges(['XM', 'MY', 'MZ'])

# Render for print 
g_12.render('img/graph_12')
```

![graph_12.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/64f92bd6-2693-1294-6d0d-d759fede3484.png)


を考えます.

このケースでZは**Bad control**とされます.

設例として, 次のようなDGPを考えます.

```math
\begin{align} \\
&X \sim \mathcal{U}(0, 1) \\
&M := x + \epsilon_m , \quad \epsilon_m \sim \mathcal{N}(0, 1) \\
&Z:= 1\{m > 0.5\}\\
&Y := 100m +\epsilon_y, \quad \epsilon_y \sim \mathcal{U}(0, 1)\\
\end{align}
```

$y \sim x$, $y \sim x + z$のそれぞれで回帰した場合の所定のメトリクスを出力してみましょう.

```python
# Set random seed for reproductivity 
np.random.seed(45)

# assign values to the variables: 
x = np.random.uniform(0, 1, N_SAMPLES)
m = x + np.random.randn(N_SAMPLES)
z = np.where(m > 0.5, 1, 0)
y = 100 * m + np.random.uniform(0, 1, N_SAMPLES)

return_reg_summary()

```
```
## Params: [ 9.27944503 86.43649901]
## p-vals: [1.35244784e-01 4.67846608e-15]
## Signif: [False  True]
## 
## Params: [-44.24459962  29.6050972  161.3400141 ]
## p-vals: [1.32257734e-026 1.48346888e-005 2.39639029e-214]
## Signif: [ True  True  True]
```


2つのモデルにおけるXの偏回帰係数の推定値は共に有意となり, DAG, DGPと整合的な結果となりました.

一方で, 2つ目のモデルにおけるXの偏回帰係数の推定値29.6050972が, 1つ目のモデルにおけるXの偏回帰係数の推定値386.43649901よりも小さな値となり, 大きく変動しました.

DAGに留意すると, 2つ目のモデルにおけるXの偏回帰係数の推定値は1つ目のモデルと比較して真のパラメータから遠ざかったと解釈することが可能です.

**mediator**の直接の子であるZを説明変数に追加すると, 関心下の説明変数の偏回帰係数の推定値が悪化しました.

標本サイズを1,000ではなく1,000,000にした時の結果も確認してみましょう.

```python
# Set random seed for reproductivity 
np.random.seed(45)

# temporarily modify the number of samples: 
tmp_n_samples = 1000000

# assign values to the variables: 
x = np.random.uniform(0, 1, tmp_n_samples)
m = x + np.random.randn(tmp_n_samples)
z = np.where(m > 0.5, 1, 0)
y = 100 * m + np.random.uniform(0, 1, tmp_n_samples)

return_reg_summary()

```
```
## Params: [ 0.55934886 99.63320567]
## p-vals: [0.00518799 0.        ]
## Signif: [ True  True]
## 
## Params: [-48.61674968  37.00067948 161.30502938]
## p-vals: [0. 0. 0.]
## Signif: [ True  True  True]
```


やはり, 2つ目のモデルにおいて, 標本サイズを増加させても, バイアスが緩和されていないことが分かります.　

### Model.16; Bad control(selection bias) 
下図のようなDAGを考えます.

```python
# Create a directed graph
g_16  = graphviz.Digraph(format='png')

## Add nodes
nodes_16 = ['Z', 'U', 'X', 'Y']
[g_16.node(n) for n in nodes_16] 

g_16.edges(['XY', 'XZ', 'UZ', 'UY'])

# Render for print 
g_16.render('img/graph_16')
```

![graph_16.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/1f7ea507-0c40-b80d-b8cc-b1bef4a53ed9.png)


このケースでZは**Bad control**とされ, この時に起こるバイアスは**selection bias**とも呼ばれます.

設例として, 次のようなDGPを考えます.

```math
\begin{align} \\
&U \sim \mathcal{U}(0, 1) \\
&X \sim \mathcal{N}(0, 1) \\
&Z:= 1.1u + x\\
&Y := u + x\\
\end{align}
```

$y \sim x$, $y \sim x + z$のそれぞれで回帰した場合の所定のメトリクスを出力してみましょう.

```python
# Set random seed for reproductivity 
np.random.seed(45)

# assign values to the variables: 
x = np.random.randn(N_SAMPLES)
u = np.random.uniform(0, 1, N_SAMPLES)
z = 1.1 * u + x
y = x + u

return_reg_summary()

```
```
## Params: [0.5041125  0.98273347]
## p-vals: [1.98120614e-305 0.00000000e+000]
## Signif: [ True  True]
## 
## Params: [4.33680869e-17 9.09090909e-02 9.09090909e-01]
## p-vals: [0.18650632 0.         0.        ]
## Signif: [False  True  True]
```


2つのモデルにおけるXの偏回帰係数の推定値は共に有意となり, DAG, DGPと整合的な結果となりました.

一方で, 2つ目のモデルにおけるXの偏回帰係数の推定値9.09090909e-02が, 1つ目のモデルにおけるXの偏回帰係数の推定値0.98273347よりも小さな値となり, 大きく変動しました.

DAGによる真のパラメータは今回の場合1であることが明快であるため, Zを説明変数に追加したときのXの偏回帰係数の推定値の精度が**悪化**したといえます.

標本サイズを1,000ではなく1,000,000にした時の結果も確認してみましょう.

```python
# Set random seed for reproductivity 
np.random.seed(45)

# temporarily modify the number of samples: 
tmp_n_samples = 1000000

# assign values to the variables: 
x = np.random.randn(tmp_n_samples)
u = np.random.uniform(0, 1, tmp_n_samples)
z = 1.1 * u + x
y = x + u

return_reg_summary()

```
```
## Params: [0.4996535  1.00000774]
## p-vals: [0. 0.]
## Signif: [ True  True]
## 
## Params: [-1.08672836e-14  9.09090909e-02  9.09090909e-01]
## p-vals: [0. 0. 0.]
## Signif: [ True  True  True]
```



やはり, 2つ目のモデルにおいて, 標本サイズを増加させても, バイアスが緩和されていないことが分かります.　

### Model.17; Bad control(selection bias)

下図のようなDAGを考えます.

```python
# Create a directed graph
g_17  = graphviz.Digraph(format='png')

## Add nodes
nodes_17 = ['Z', 'X', 'Y']
[g_17.node(n) for n in nodes_17] 

g_17.edges(['XY', 'XZ', 'YZ'])

# Render for print 
g_17.render('img/graph_17')
```

![](img/graph_17.png)

[モデル16](#model_16)と同様, このケースでZは**Bad control**とされ, この時に起こるバイアスは**selection bias**とも呼ばれます.

設例として, 次のようなDGPを考えます.

```math
\begin{align} \\
&X \sim \mathcal{N}(0, 1) \\
&Y := 3x + \epsilon_y, \quad \epsilon_y \sim \mathcal{U}(0, 1)\\
&Z:= x + y\\
\end{align}
```

$y \sim x$, $y \sim x + z$のそれぞれで回帰した場合の所定のメトリクスを出力してみましょう.

```python
# Set random seed for reproductivity 
np.random.seed(45)

# assign values to the variables: 
x = np.random.randn(N_SAMPLES)
y = 3 * x + np.random.uniform(0, 1, N_SAMPLES)
z = x + y

return_reg_summary()

```
```
## Params: [0.5041125  2.98273347]
## p-vals: [1.98120614e-305 0.00000000e+000]
## Signif: [ True  True]
## 
## Params: [ 4.82253126e-16 -1.00000000e+00  1.00000000e+00]
## p-vals: [7.47868963e-18 0.00000000e+00 0.00000000e+00]
## Signif: [ True  True  True]
```


2つのモデルにおけるXの偏回帰係数の推定値は共に有意となり, DAG, DGPと整合的な結果となりました.

一方で, 2つ目のモデルにおけるXの偏回帰係数の推定値-1.00000000e+00が, 1つ目のモデルにおけるXの偏回帰係数の推定値2.98273347よりも小さな値となり, 大きく変動しました.

DAGによる真のパラメータは今回の場合3であることが明快であるため, Zを説明変数に追加したときのXの偏回帰係数の推定値の精度が**悪化**したといえます.

標本サイズを1,000ではなく1,000,000にした時の結果も確認してみましょう.

```python
# Set random seed for reproductivity 
np.random.seed(45)

# temporarily modify the number of samples: 
tmp_n_samples = 1000000

# assign values to the variables: 
x = np.random.randn(tmp_n_samples)
y = 3 * x + np.random.uniform(0, 1, tmp_n_samples)
z = x + y

return_reg_summary()

```
```
## Params: [0.4996535  3.00000774]
## p-vals: [0. 0.]
## Signif: [ True  True]
## 
## Params: [-3.06276272e-15 -1.00000000e+00  1.00000000e+00]
## p-vals: [0. 0. 0.]
## Signif: [ True  True  True]
```


やはり, 2つ目のモデルにおいて, 標本サイズを増加させても, バイアスが緩和されていないことが分かります.　

### Bad controlのまとめ

このセクションで説明したDAG

![graph_7.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/a4a25275-36e9-b85d-9ef7-f467602efa39.png)
![graph_10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/71e38bec-9809-f077-5d24-806c3f1af0fc.png)
![graph_11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/97fa84d1-70f6-92d7-dd1b-8a7d2da16756.png)
![graph_12.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/92db547c-29b8-61fa-d61a-7995cc204e79.png)
![graph_16.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/2259c866-d967-e78b-3d20-c4c070415640.png)
![graph_17.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/6936020e-bd69-0791-3a2a-bc9f16541340.png)

のいずれにおいても, Zが漸近的に**Bad control**であることが確認できました.

さらにこれらは, 標本サイズを増加させてもバイアスが緩和されないことを確認しました.

とりわけモデル7は, X-Y間の因果関係が存在しないDAGでしたが, 例で用いた DGPでは, Zを説明変数に追加すると見せかけの相関が発生することが確認できました.

# まとめ

特にまとめることはないですが参考になれば幸いです.