はじめに
本稿では"A Crash Course in Good and Bad Controls"というタイトルの文献(以下リンク)について, 簡単な具体例を用いてゆっくり確認していきます.
なお, 本稿で頻出するcontrol(statistical control)という用語は, 私が調べる限り, 関心下の説明変数が目的変数に与える効果を検証する際に, 関心のない変数の影響を固定するために敢えて説明変数に追加されるもののことを指すようです.
関連ライブラリのインポート
はじめに, 以下のように関連ライブラリをインポートしておきます.
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を描画します.
以下のようなコードを実行すると,
# 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
を作ることができました. ここで, Y, X, Zはそれぞれ被説明変数, 関心下の説明変数, コントロールしたい変数を表します.
つまり, Good controls, Bad controlsのいずれに該当するのか判断したいのは変数Zです.
実際にGood controlであることを確認したいので, このDAGに対して, 次のようなData generating process(DGP)を仮定しましょう.
\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に従って, 以下のコードのように変数をアサインします.
# 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%としたとき, 偏回帰係数が有意かどうか
を取得し, 表示させます.
それでは以下のとおりコードを実行して結果を出力させます.
# 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において,
\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
# 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')
を仮定します.
先ほどのDAGとは異なり, 観測できない(Unobservable)変数Uを追加しました.
更に, DGPとして,
\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で回帰するモデルを作成し, メトリクスを出力します.
# 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)についても同様のことが言えます.
# 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')
\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}
前述の手順を実施すると, 以下のようになり,
# 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
# 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')
を考えます.
これまでのDAGとは異なり, XはMを介してYと因果関係があります.
このDAGに基いて, XとYをいずれもバイナリ変数とした以下のDGP
\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$で回帰した場合のそれぞれについて所定のメトリクスを出力すると,
# 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が存在するケース
# 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')
について考えます.
更に, 先ほどのDGPを少し手直しした,
\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で回帰した場合のそれぞれについて所定のメトリクスを出力すると,
# 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を考えます.
# 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')
\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で回帰した場合のそれぞれについて所定のメトリクスを出力すると,
# 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
のいずれにおいても, ZがGood controlであることが確認できました.
これらのZはcommon cause又はforkとも呼ばれるようです。
Neutral controls
Model.8; Neutral control
下図のようなDAG, DGPを考えます.
# 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')
\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で回帰した場合のそれぞれについて所定のメトリクスを出力すると,
# 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に変更してみます.
# 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を考えます.
# 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')
\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で回帰した場合のそれぞれについて所定のメトリクスを出力すると,
# 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に変更すると, 結果は以下のようになります.
# 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を考えます.
# 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')
このケースにおけるZはモデル8, モデル9と同様に,漸近的にはNeutral controlとされます.
さらに, モデル8と同様, 標本サイズが有限なレジームにおいては, 上図のZを説明変数に追加することにより, XがYに与える影響の推定精度が向上する可能性があります.
設例として, 次のようなDGPを考えます.
\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$のそれぞれで回帰した場合の所定のメトリクスを出力してみましょう.
# 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に変更すると, 結果は以下のようになります.
# 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を考えます.
# 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')
このケースにおけるZはNeutral controlとされます.
なお, このZはXのPost treatment variableと呼ばれます.
さらに, 標本サイズが有限なレジームにおいて, 上図のZをコントロールすることでXのバリエーションが制限されるような状況では, XがYに与える影響の推定精度が悪化する可能性があります.
設例として, 次のようなDGPを考えます.
\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$のそれぞれで回帰した場合の所定のメトリクスを出力してみましょう.
# 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に変した結果も確認してみましょう.
# 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の標本サイズにおいて,
# 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を考えます.
# 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')
モデル14と同様に, このZはXのPost treatment variableと呼ばれます.
また, このケースにおけるZはNeutral controlとされます. さらに, 標本サイズが有限なレジームにおいて, 上図のZをコントロールすることでXのバリエーションが制限されるような状況では, XがYに与える影響の推定精度が悪化する可能性があります.
設例として, 次のようなDGPを考えます.
\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$のそれぞれで回帰した場合の所定のメトリクスを出力してみましょう.
# 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に変した結果も確認してみましょう.
# 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
のいずれにおいても, 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を考えます.
# 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')
このケースでZはBad controlとされ, このDAGにおけるZはM-biasとも呼ばれます.
設例として, 次のようなDGPを考えます.
\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$のそれぞれで回帰した場合の所定のメトリクスを出力してみましょう.
# 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にした時の結果も確認してみましょう.
# 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
# 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')
\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$のそれぞれで回帰した場合の所定のメトリクスを確認します.
# 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にした時の結果も確認してみましょう.
# 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を考えます.
# 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')
このケースでZはBad controlとされ, このDAGにおけるZはpre-treatmentに該当します.
さらに言えば, 交絡として除去すべきUが未観測であるような事例であると言えます.
設例として, 次のようなDGPを考えます.
\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$のそれぞれで回帰した場合の所定のメトリクスを出力してみましょう.
# 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にした時の結果も確認してみましょう.
# 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を考えます.
# 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')
このケースでZはBad controlとされ, このDAGにおけるZはMediatorに該当します.
設例として, 次のようなDGPを考えます.
\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$のそれぞれで回帰した場合の所定のメトリクスを出力してみましょう.
# 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
# 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')
を考えます.
このケースでZはBad controlとされます.
設例として, 次のようなDGPを考えます.
\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$のそれぞれで回帰した場合の所定のメトリクスを出力してみましょう.
# 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にした時の結果も確認してみましょう.
# 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を考えます.
# 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')
このケースでZはBad controlとされ, この時に起こるバイアスはselection biasとも呼ばれます.
設例として, 次のようなDGPを考えます.
\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$のそれぞれで回帰した場合の所定のメトリクスを出力してみましょう.
# 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にした時の結果も確認してみましょう.
# 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を考えます.
# 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')
モデル16と同様, このケースでZはBad controlとされ, この時に起こるバイアスはselection biasとも呼ばれます.
設例として, 次のようなDGPを考えます.
\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$のそれぞれで回帰した場合の所定のメトリクスを出力してみましょう.
# 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にした時の結果も確認してみましょう.
# 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
のいずれにおいても, Zが漸近的にBad controlであることが確認できました.
さらにこれらは, 標本サイズを増加させてもバイアスが緩和されないことを確認しました.
とりわけモデル7は, X-Y間の因果関係が存在しないDAGでしたが, 例で用いた DGPでは, Zを説明変数に追加すると見せかけの相関が発生することが確認できました.
まとめ
特にまとめることはないですが参考になれば幸いです.