「回帰分析から学ぶ計量経済学」をPythonで写経 ～ 第７章「反実仮想の世界」③回帰分析、操作変数法

第７章「反実仮想の世界」

書籍の著者　山澤成康 先生

この記事は、書籍「回帰分析から学ぶ計量経済学」第７章「反実仮想の世界」の Python写経活動 を取り扱います。

第７章は反実仮想、因果推論に取り組みます。
今回のテーマは 回帰分析 と 操作変数法 です。
では書籍を開いて回帰分析の旅に出発です🚀

はじめに

書籍「回帰分析から学ぶ計量経済学」のご紹介

このシリーズは書籍「回帰分析から学ぶ計量経済学: Excelで読み解く経済のしくみ」（オーム社、「テキスト」と呼びます）の Python 写経です。

テキストは、2023年11月に発売され、副題「Excelで読み解く経済のしくみ」のとおり、主に Excel を用いて、計量経済学を平易に学べる素晴らしい書籍です。
テキストの「はじめに」に著者の先生が執筆の動機を書かれています。

私にとって計量経済学は高嶺の花ですが、このテキストでさまざまな回帰分析のアプローチを知ることができました。
また、書籍の Excel 処理を Python に置き換える「寄り道写経」の実践を通じて、回帰分析のお気持ちに少し近づけた感じがいたします。

回帰分析に慣れ親しむのに丁度良いレベル感と内容ですので、これはぜひともブログにしたい！と思って現在に至ります。
計量経済学の色を薄め、データ分析の色を濃いめに書いてまいります！

引用表記

この記事は、出典に記載の書籍に掲載された文章と配布データを引用し、適宜、掲載文章・配布データを改変して書いています。
【出典】
「回帰分析から学ぶ計量経済学: Excelで読み解く経済のしくみ」
第１版第１刷、著者 山澤成康、オーム社

記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！

第７章 反実仮想の世界

この記事は第７章の以下の節を取り扱います。

記事に用いるデータは、オーム社の書籍紹介サイトからダウンロードできる Excel ファイル内のデータをもとにしてCSVファイルを作成し、data フォルダに格納しています。

第７章で用いるライブラリをインポートします。

### インポート

# 数値計算
import numpy as np
import pandas as pd

# 統計処理
import scipy.stats as stats
import pingouin as pg
import statsmodels.formula.api as smf

# 効果検証
from linearmodels import IV2SLS                        # 操作変数法
import rdrobust                                        # 回帰不連続デザイン
import rddensity                                       #  同上
from causallib.estimation import PropensityMatching    # 傾向スコアマッチング
from causallib.estimation import IPW                   # 傾向スコア・IPW法
from causallib.evaluation import evaluate              # 傾向スコアの評価

# 機械学習
from sklearn.linear_model import LogisticRegression    # ロジスティック回帰

# グラフ描画
from graphviz import Digraph

# 描画
import matplotlib.pyplot as plt
import seaborn as sns
import japanize_matplotlib

7.7 回帰分析

テキストは「回帰分析によっても、条件がそろえば因果関係の分析ができます」とし、職業訓練（原因変数$${T}$$）が年収（結果変数$${Y}$$）に影響を与えるかどうかを仮想データで分析しています。

■ データの確認 p.229
データを読み込みます。

### データの読み込み
# 年収⇒万円、訓練⇒1:職業訓練を受けた、2:受けなかった

df7 = pd.read_csv('./data/07_05_jobtrain.csv')
print('df7.shape:', df7.shape)
display(df7)

【実行結果】
10名分のデータであり、年収（万円）が結果変数$${Y}$$、訓練が 0：未実施、1：実施を示す原因変数$${T}$$です。

なお、年齢、就学年数は共変量として扱われます。
共変量はテキストによると「原因変数$${T}$$以外で結果変数$${Y}$$の両方に影響する変数。原因変数と結果変数以外の変数全てという意味でも使われる」です。

ではテキストにならって、年収を目的変数、訓練（有無）を説明変数にして回帰分析を実行します。

### 回帰分析の実行 年収＝定数項＋β×職業訓練ダミー
# 回帰結果は芳しくない：決定係数が小さい、F値のp値が大きい、係数のp値が大きい

result = smf.ols(formula='年収 ~ 訓練', data=df7).fit()
display(result.summary())

【実行結果】
訓練（有無）の係数の推定値は$${-23.2}$$です。
訓練を実施すると年収が$${23.2}$$万円下がるって、何か変ですね。。。
ただ、係数の$${p}$$値は$${0.724}$$であり、有意ではなさそうです。

職業訓練の有無別の各変数の平均値を見ます。

### 職業訓練を受けた／受けなかった人の平均
display(df7.groupby(by=['訓練'])[['年収', '年齢', '就学年数']].mean().round(1))

【実行結果】
訓練を受けた（=1）のほうが年収が低くなっています。
共変量である年齢と就学年数が、訓練を受けたほうが小さな値であることがどうやら、因果推論の鍵を握っているようです。

データを可視化してみましょう。

### 散布図の描画
sns.scatterplot(data=df7, x='年齢', y='年収', alpha=0.7, 
                hue='訓練', palette=['tab:blue', 'tab:red'],
                size='就学年数', sizes=(40, 100));

【実行結果】
横軸に年齢、縦軸に年収をとり、散布図の点は青色が訓練未実施、赤色が訓練実施です。
また、点の大きさは就学年数の大きさに比例しています。

訓練を実施したのは $${35}$$歳以下、かつ、就学年数$${16}$$年以下に「偏って」います。
$${35}$$歳以下、かつ、就学年数$${16}$$年以下に限定してみると、訓練を受けたほう（赤い点）が年収が高くなっています。
また、訓練を受けたうち、就学年数が高いほうが年収が高くなっているようにも見えます。

共変量である年齢と就学年数を重回帰分析の変数に加えることで、偏り＝セレクション・バイアスの影響を「調整」（コントロール）することができるそうです。

■ 共変量を加えた回帰分析による訓練の効果の推定 p.230
では早速、訓練（有無）、年齢、就学年数を説明変数にした回帰分析を実行しましょう。

### 回帰分析の実行 年収＝定数項＋β1×職業訓練ダミー + β2×年齢 + β3×就学年数
# 共変量を加える・・・共変量の影響を取り除く

result = smf.ols(formula='年収 ~ 訓練 + 年齢 + 就学年数', data=df7).fit()
display(result.summary())

【実行結果】
訓練（有無）の係数が大きく変わりました！
年齢・就学年数が同じと仮定すると、訓練を実施することで年収が$${100.7}$$万円増加すると推定されました。

ちなみに、自由度調整済み決定係数は$${1}$$（！）と当てはまりがよく、訓練（有無）・年齢・就学年数の係数推定値の$${p}$$値はほぼ０、$${5\%}$$水準で有意となっています。

7.9 操作変数法の具体例

操作変数法の説明はやや難しいです。。。

■ 操作変数法とは p.231
テキストは年収に影響する変数として「教育：就学年数」と「持って生まれたもの：個人の能力」を挙げています。
個人の能力は、測定するのが難しく、かつ、年収と就学年数の両方に影響を与える交絡因子となっています。

図示します。
graphviz ライブラリを利用します。

### DAGの可視化 by graphviz

## 設定
g = Digraph(engine='neato')

## node：頂点の作成、posで位置固定, shapeで形状を指定
g.node('X就学年数', pos='2, 1!', shape='oval')   # 処置変数:就学年数
g.node('Y年収', pos='0, 1!', shape='oval')       # 結果変数:年収
g.node('U個人の能力', pos='2, 0!', shape='oval') # 共変量:個人の能力

## edge：辺の作成
g.edge('U個人の能力', 'Y年収')
g.edge('U個人の能力', 'X就学年数')
g.edge('X就学年数', 'Y年収')

## グラフの描画
g

【実行結果】

個人の能力が就学年数に与える影響を断ち切り、就学年数が年収に与える効果を正確に把握しようという方法が「操作変数法」です。
「操作変数」を導入し、第１段階の回帰分析では操作変数で説明変数を推定、第２段階の回帰分析では推定した説明変数で目的変数を推定するものです。

テキストは操作変数に「大学までの距離」を採用しています。
大学までの距離が近ければ就学年数が長くなるという可能性と、大学までの距離は個人の能力と関係ないという点がポイントです。

図示します。

### DAGの可視化 by graphviz

## 設定
g = Digraph(engine='neato')

## node：頂点の作成、posで位置固定, shapeで形状を指定
g.node('Z大学までの距離', pos='5, 1!', shape='oval')  # 操作変数:大学までの距離
g.node('X就学年数', pos='2, 1!', shape='oval')       # 処置変数:就学年数
g.node('Y年収', pos='0, 1!', shape='oval')           # 結果変数:年収
g.node('U個人の能力', pos='2, 0!', shape='oval')      # 共変量:個人の能力

## edge：辺の作成
g.edge('Z大学までの距離', 'X就学年数')
g.edge('U個人の能力', 'Y年収')
g.edge('X就学年数', 'Y年収')

## グラフの描画
g

【実行結果】

■ 仮想データの確認
テキストの仮想データは以下の式で作られています。

データを読み込みます。

### データの読み込み
# 年収⇒万円、本人の能力100～ランダム、大学までの距離1～5までランダム
# 年収＝就学年数×30＋本人の能力×3＋標準正規分布乱数で作成
# 就学年数＝本人の能力＋大学までの距離×2で算出

df8 = pd.read_csv('./data/07_06_wage.csv')
print('df8.shape:', df8.shape)
display(df8)

【実行結果】
10人分のデータです。

データの相関係数を見てみましょう。

### 相関係数の表示
display(df8.corr().round(3))

【実行結果】
年収と就学年数・本人の能力に強い正の相関関係が見られます。
就学年数と本人の能力・大学までの距離にも強めの正の相関関係が見られます。

### 相関係数のヒートマップの描画
sns.heatmap(data=df8.corr(), annot=True, cmap='Greens',
            annot_kws=dict(fontsize=14))

【実行結果】

■ 本人の能力が測定できると仮定したときの回帰分析 p.232
本来は測定できない本人の能力が測定できたと仮定したときの回帰分析です。

### 回帰分析の実行 年収、就学年数、本人の能力
# データ生成の仮定に合った結果を得られるが、本人の能力データは現実的に入手不能

result = smf.ols(formula='年収 ~ 就学年数 + 本人の能力', data=df8).fit()
display(result.summary())

【実行結果】
係数に着目します。
就学年数は真値$${30}$$に近い$${30.5}$$、本人の能力は真値$${3}$$に近い$${2.97}$$です。

■ 測定不能な本人の能力を説明変数から除外する回帰分析 p.233
本人の能力を説明変数から除外して、年収を目的変数、就学年数だけを説明変数にする回帰分析では、脱落変数バイアスが生じるそうです。
脱落変数は共変量のうち、変数として採用されていないか、観測できないものです。

### 回帰分析の実行 年収、就学年数
# 本人の能力は現実的に入手不能のため説明変数から外す
# 就学年数がデータ生成仮定から乖離してしまう

result = smf.ols(formula='年収 ~ 就学年数', data=df8).fit()
display(result.summary())

【実行結果】
就学年数の係数の推定値は$${45.6}$$となり、真値から大きく乖離してしまいました。

■ 操作変数法による推定 p.233
操作変数法で就学年数の影響を推定します。
２段階の回帰分析を実行します。

まずは第１段階目の回帰分析。
就学年数を目的変数、大学までの距離を説明変数にします。

### 操作変数法による２段階の回帰分析

## １段階目：就学年数を大学までの距離で回帰して推定値を得る
result1 = smf.ols(formula='就学年数 ~ 大学までの距離', data=df8).fit()
# 就学年数の推定値を取得
studty_est = result1.fittedvalues
# 回帰分析の結果を表示
display(result1.summary())

【実行結果】

続いて、上のモデルを用いて就学年数の予測値を推定します。

## ２段階目：年数を就学年数の推定値で回帰して推定値を得る

# データの前処理
df8reg = df8.copy()
df8reg['就学年数_推定値'] = studty_est
display(df8reg)

【実行結果】
最右列が就学年数の推定値です。

最後に、年収を目的変数、就学年数の推定値を説明変数にして、第２段階目の回帰分析を実行します。

# 回帰分析の実行
result2 = smf.ols(formula='年収 ~ 就学年数_推定値', data=df8reg).fit()
# 回帰分析の結果を表示
display(result2.summary())

【実行結果】
就学年数の係数の推定値は$${30.5}$$。真値$${30}$$に近似しました！
就学年数（１年）が年収に与える影響は$${30.5}$$です。

■ linearmodels ライブラリの操作変数法
linearmodels ライブラリの IV2SLS を用いることで、２段階の回帰分析を「１度に」実行できます。

### linearmodelsのIV2SLSを用いた操作変数法の実行
# 参考：Pythonで学ぶ効果検証入門 152ページ
# 就学年数の係数30.5万円(p値=0.003)はデータ生成の仮定に合った結果である

iv = IV2SLS.from_formula(
    formula='年収 ~ 1 + [就学年数 ~ 大学までの距離]', data=df8).fit()
display(iv.summary)

【実行結果】
就学年数の係数の推定値は$${30.5}$$です！うまく推定できています！
就学年数（１年）が年収に与える影響は$${30.5}$$です。

自由度調整済み決定係数は$${0.78}$$であり、まずまずの当てはまり具合です。
また、就学年数の係数の推定値の$${p}$$値は$${0.0037}$$であり、$${5\%}$$水準で有意です。

今回の写経は以上です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。