「回帰分析から学ぶ計量経済学」を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.10 震災の復興需要（パネルデータ）

宮城県の復興需要の大きさを推定します。
次の３ステップで実践します。

• 震災前の県内総生産の動きが宮城県と似ていて、かつ、震災の影響が小さな都道府県をチェック（答えは「高知県」）

• 震災前の宮城県と高知県の回帰分析を実行

• 震災の影響が無い宮城県の県内総生産（ベースライン）の推定、および、ベースラインと実際の宮城県の県内総生産の差を復興需要として把握

■ データの確認・動きが似ている都道府県を探す p.234
データを読み込みます。

### データの読み込み
# 宮城県と高知県は2006～2009年度までの相関係数が高い

df9 = pd.read_csv('./data/07_07_restoration.csv', index_col=0)
print('df9.shape:', df9.shape)
display(df9)

【実行結果】

2006～2009年度の両県の相関係数を確認します。

### 2006～2009年度までのデータの相関係数
display(df9.loc[2006:2009].corr())

【実行結果】
確かに宮城県と高知県の県内総生産の間に強い正の相関関係があります。

■ 回帰分析の実行 p.234
テキストの回帰の式に基づき、回帰分析を実行します。

### 回帰分析の実行

# 2006～2009年までのデータで回帰分析を実行
result = smf.ols(formula='宮城県 ~ 高知県', data=df9.loc[2006:2009]).fit()
# 回帰分析の結果を表示
display(result.summary())

【実行結果】
推定した回帰モデルです。

■ ベースラインの推定、復興需要の推定p.235
このモデルを用いて、2006～2018年度の宮城県の「復興需要の影響が無い」県内総生産を推定します。
テキストはこの県内総生産の推定値を「ベースライン」と呼んでいます。
あわせて、宮城県の実際の県内総生産とベースラインの差＝復興需要の大きさを計算します。

### 回帰分析の係数で宮城県を推定

# 推定値を予測
df9['宮城県_推定値'] = result.predict(df9['高知県'])
# 実績値と推定値の差を計算
df9['復興需要_推定値'] = df9['宮城県'] - df9['宮城県_推定値']
# 結果の表示
display(df9.round(0))

【実行結果】
宮城県の推定値がベースライン、復興需要の推定値が宮城県の実際の県内総生産とベースラインの差です。

これらの数値を可視化しましょう。

### 可視化

# 描画領域の設定
plt.figure(figsize=(8, 3))
# ３変数の折れ線グラフの描画
sns.lineplot(data=df9[['宮城県', '宮城県_推定値', '復興需要_推定値']] / 10000)
# 修飾
plt.axhline(0, color='black', lw=0.7)
plt.ylim(-100, 1000)
plt.ylabel('県内総生産 [兆円]')
plt.grid(lw=0.5)
plt.show()

【実行結果】
青い線が宮城県の実際の県内総生産、オレンジの点線が復興需要が無いと仮定した時の宮城県の県内総生産の推定値（ベースライン）、緑の点線が復興需要の推定値です。
2012年度から復興需要による県内総生産の増加（影響）が増加基調で見られます。

7.10 回帰不連続デザイン

■ 回帰不連続デザイン p.236
テキストによると「回帰不連続デザインは、政策が実行された日の前後で、変数に非連続的な変化が起こることを利用して、政策の効果を測るもの」です。

テキストの例題は次のシナリオで「A高校の学力増進効果」を推定しています。

• A高校の入試の最低合格点を300点とする。

• B君は300点でA高校に合格した。

• C君は299点でA高校に合格できず、他の高校に行った。

• B君とC君の大学入試時の成績を見ることでA高校の学力増進効果があったかどうかを分析する。

300点前後の回帰分析の定数項の変化＝水準の変化が A高校の学力増進効果です。

仮想データを生成します。
このコードの乱数生成結果がテキストと異なるため、テキストの推定結果と異なります。

### データの作成

## 設定と準備
N = 30                                # 標本サイズ
score = 300                           # A高校の合格最低点
rng = np.random.default_rng(seed=0)   # 乱数生成器

## データの作成
# 高校入試の得点
x = np.arange(286, 286 + N)
# A高校の合否
y = np.where(x >= 300, 1, 0)
# 大学入試時の成績
base = np.hstack([np.arange(2, 16), np.arange(25, 41)])
z = base + 40 + rng.standard_normal(size=N)

## データフレーム化
df10 = pd.DataFrame(dict(高校入試の得点=x, A高校の合否=y, 大学入試時の成績=z))
print('df10.shape: ', df10.shape)
display((df10[(df10['高校入試の得点'] >= 297) & (df10['高校入試の得点'] <= 303)]))

【実行結果】
高校入試の得点が 299 点と 300 点との間で非連続な変化が起きています。

仮想データを可視化します。

### 可視化

# 回帰直線付き散布図の描画
sns.lmplot(data=df10, x='高校入試の得点', y='大学入試時の成績', 
           hue='A高校の合否', palette=['royalblue', 'tab:green'],
           height=3, aspect=1.5, line_kws=dict(color='tomato'))
# A高校合格最低点の垂直線の描画
plt.axvline(300, color='tab:red', ls='--')
# 修飾
plt.grid(lw=0.5)

【実行結果】
300 点をしきい値（カットオフ値）とし、299点以下と300点以上の２本の回帰直線を見ると、大学入試時の成績の「水準」が違う感じがいたします。
この水準の変化＝定数項（切片）の変化が A高校の学力増進効果です。

書籍「Pythonで学ぶ効果検証入門」と「因果推論」を参考にして、A高校の学力推進効果を推定します。

まずは、カットオフ値 300 点付近の分布が連続的になっているかを確認します。
ヒストグラムを描画して連続性を確認します。

### カットオフ値付近の強制変数の分布の確認

# ビンの幅の設定
bins=np.arange(285, 317)
# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 4))
# 高校入試の得点のヒストグラムの描画
sns.histplot(data=df10, x='高校入試の得点', bins=bins, alpha=0.3,
             hue='A高校の合否', palette=['royalblue', 'green'], ax=ax)
ax.axvline(300, color='tab:red', lw=3, ls='--')
# 凡例処理
sns.move_legend(ax, loc='upper left', bbox_to_anchor=(1, 1),
                labels=['否', '合格'])
# 修飾
ax.set(title='高校入試の得点のヒストグラム')
plt.show()

【実行結果】
めっちゃ一様分布してます！
赤い点線の 300 点付近の分布は同じであり「連続性あり」と考えられます。

さらにカットオフ値付近の連続性の検定（McCraryの検定）で確認します。
rddensity ライブラリの rddensity で連続性の検定を実行します。
帰無仮説は「連続である」です。

### カットオフ値付近の連続性の検定 McCraryの検定
# 参考: 書籍「因果推論」 p.103
# 帰無仮説：カットオフ値近傍において「確率密度関数が連続である」
# MethodのRobustのp値1.0・・・帰無仮説を棄却できず、カットオフ値付近は「連続である」

rddensity.rddensity(X=df10['高校入試の得点'], c=300)

【実行結果】

注目点は「Method:」「Robust」の「P > |T|」（$${p値}$$）であり、値は$${1.0}$$です。
$${5\%}$$水準で帰無仮説を棄却できず、データは「連続である」を受容します（カットオフ値付近で連続性がある、ということ）。

カットオフ値付近の連続性を確認できました。
続いて、A高校の学力推進効果を推定します。
rdrobust ライブラリの rdrobust を利用します。

### Sharp RDDによる効果の推定
# 参考: Pythonで学ぶ効果検証入門 p.239
# Robustの係数、標準誤差、t値・p値を見る

result = rdrobust.rdrobust(
    y=df10['大学入試時の成績'], x=df10['高校入試の得点'], c=300, all=True)
display(result)

【実行結果】

注目点は「Method」「Robust」の「Coef.」です。
A高校の学力推進効果の推定値は $${8.084}$$ 点です。
ただ、推定値の$${p}$$値は$${0.09}$$であり、$${5\%}$$水準で有意とは言えないという結果になりました。。。（$${10\%}$$水準なら有意？）

最後に rdrobust ライブラリの rdplot を利用して、回帰不連続デザインの RD プロットを描画します。

### RDプロットの描画
# 参考: 書籍「因果推論」 p.107
rdrobust.rdplot(y=df10['大学入試時の成績'], x=df10['高校入試の得点'], c=300,
                title='Causal Effects of A High School',
                x_label='High School Entrance Exam Scores',
                y_label='Grades at University Entrance Exam')

【実行結果】
カットオフ値 300 点付近の不連続になっている部分が、局所的な平均因果効果を示しています。

■ 回帰屈折デザイン p.237
ある時点（カットオフ値）前後で、水準が変わるのではなく（つまり回帰不連続デザインではない）、傾きが変わるものを回帰屈折デザイン、または、回帰ねじれデザインと呼ぶそうです。

テキストは新型コロナウイルス感染拡大初期の緊急事態宣言（東京都）が新規感染者数にどの程度影響したかを分析しています。
・緊急事態宣言発令：2020年4月7日
・緊急事態宣言解除：2020年5月25日

データを読み込みます。

### データの読み込み

df11 = pd.read_csv('./data/07_08_infected.csv', index_col=0, parse_dates=[0])
print('df11.shape:', df11.shape)
display(df11)

【実行結果】

データの前処理を行います。
・トレンド変数の追加（2020年2月15日を１とし、日々１カウントアップ）
・新規感染者数の対数変換
・緊急事態宣言ダミー変数の追加（0：宣言未出、1：宣言発出）

### データの前処理

# トレンド変数の追加 1～連番
df11['トレンド'] = range(1, len(df11)+1)
# 新規感染者の対数の追加 新規感染者=0の場合は0を設定
df11['新規感染者_対数'] = (
    df11['新規感染者'].apply(lambda x: np.log(x) if x>0 else 0))
# 緊急事態宣言発令あり・なしの列の追加
df11['緊急事態宣言'] = np.where(df11.index >= '2020-04-07', 1, 0)
# 結果の表示
print('df11.shape:', df11.shape)
display(df11)

【実行結果】

緊急事態宣言の発出前（～2020年4月6日）と発出後（2020年4月7日～5月25日）の２つの期間別に回帰分析を実行します。
まずは発出前。

### 2020/3/24～2020/4/6の期間の回帰分析の実行

# 2020/3/24～2020/4/6の期間のデータを抽出
df11_before = (
    df11[(df11.index >= '2020-03-24') & (df11.index <= '2020-04-6')].copy())
# 回帰分析の実行
result = smf.ols(formula='新規感染者_対数 ~ トレンド', data=df11_before).fit()
# 回帰係数の取得
intercept_b, slope_b = result.params
# 回帰の結果の表示
display(result.summary())

【実行結果】
発出前のトレンド変数の傾きは$${0.1041}$$です（$${5\%水準で有意}$$）。

続いて発出後。

### 2020/4/7～2020/5/25の期間の回帰分析の実行

# 2020/4/7～2020/5/25の期間のデータを抽出
df11_after = (
    df11[(df11.index >= '2020-04-07') & (df11.index <= '2020-05-25')].copy())
# 回帰分析の実行
result = smf.ols(formula='新規感染者_対数 ~ トレンド', data=df11_after).fit()
# 回帰係数の取得
intercept_a, slope_a = result.params
# 回帰の結果の表示
display(result.summary())

【実行結果】
発出後のトレンド変数の傾きは$${-0.0808}$$です（$${5\%水準で有意}$$）。

緊急事態宣言の発出によって、新規感染者数のトレンドは増加から減少に転じているようです。

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

### 可視化

# 緊急事態宣言発令前、発令中の予測値の算出
pred_before = (
    intercept_b + slope_b * df11[df11['緊急事態宣言'] == 0]['トレンド'])
pred_after = (
    intercept_a + slope_a * df11[df11['緊急事態宣言'] == 1]['トレンド'])

# 描画領域の設定
fig, ax = plt.subplots(figsize=(8, 4))
# 新規感染者の散布図の描画
sns.scatterplot(data=df11, x=df11.index, y='新規感染者_対数', s=50, 
                hue='緊急事態宣言', palette=['tab:blue', 'coral'], alpha=0.8, 
                legend=False, ax=ax)
# 緊急事態宣言前の回帰直線の描画
sns.lineplot(x=df11[df11['緊急事態宣言'] == 0].index, y=pred_before,
             color='tab:red')
# 緊急事態宣言中の回帰直線の描画
sns.lineplot(x=df11[df11['緊急事態宣言'] == 1].index, y=pred_after,
             color='tab:red')
# 緊急事態宣言発令日の垂直点線の描画
ax.axvline(18358.5, color='tab:red', ls='--')
# 修飾
ax.tick_params(axis='x', labelrotation=30)
ax.set_ylim(-1, 7)
ax.grid(lw=0.5)
plt.show()

【実行結果】
赤い垂直点線が 緊急事態宣言発出日の2020年4月7日です。
緊急事態宣言の発出前は増加トレンド、発出後は減少トレンドになっています。
緊急事態宣言の発出は、新規感染者数（対数ベース）増加の抑制に効果があったように思えます。

今回の写経は以上です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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の教科書です。
よかったらぜひ、お試しくださいませ。

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