「回帰分析から学ぶ計量経済学」を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.6 差の差の分析 p.224～

テキストは「実験が難しい場合、さまざまな統計から政策効果を把握する方法があります。差の差分析と呼ばれるものです」としています。
政策を実施した自治体と実施しなかった自治体のデータが入手できる場合に差の差分析が可能とし、２つの自治体グループに違いがあれば政策に効果があったと考えるようです。
手続き的に「２回、差を取る」ので「差の差」なのだそう。
・自治体ごとに政策実施前後の差を取る
・２つの自治体の差の差を取る

テキストの「まん延防止等重点措置」の実施による県内総生産の影響を追いかけます。

データを作成します。

### データの作成

# データの作成
df4 = pd.DataFrame(
    dict(県=['C県-措置あり', 'D県-措置なし'], 実施前=[500, 300], 実施後=[400, 250]))
df4 = df4.set_index(['県'])
# 結果の表示
display(df4)

【実行結果】
C県はまん延防止等重点措置を実施、D県は未実施です。
数値は県内総生産です。

差の差を算出します。

### 差の差の算出
df4['差'] = df4['実施後'] - df4['実施前']
df4 = pd.concat([df4, (df4.iloc[0] - df4.iloc[1]).rename('差').to_frame().T])
display(df4)

【実行結果】
差の差は$${-50}$$。
まん延防止等重点措置の実施による県内総生産への影響は$${-50}$$（50の減少）です。

7.7 平均値の差の検定 p.226～

テキストは２つの仮想データを用いて、少子化対策の効果を推定します。

■ 少子対策を実施した10の自治体の合計特殊出生率 p.226～
どの自治体も人口や経済規模が同じと仮定して、合計特殊出生率が単純に平均できるものとして、分析をします。
対応のある２標本のｔ検定を利用します。

データを読み込みます。

### データの読み込み

df5 = pd.read_csv('./data/07_03_child.csv', index_col=0)
df5['差'] = df5['政策後'] - df5['政策前']
print('df5.shape:', df5.shape)
display(df5)

【実行結果】

要約統計量を算出します。

### 要約統計量の算出
display(df5.describe().round(2))

【実行結果】
政策前の合計特殊出生率の平均値は$${1.36}$$、政策後の平均値は$${1.46}$$。
政策実施により合計特殊出生率が$${0.1}$$上昇しています。

各自治体の施策前後の合計特殊出生率の変化を可視化してみます。

### 可視化
sns.lineplot(data=df5.iloc[:, :-1].T)
plt.legend(bbox_to_anchor=(1, 1), title='自治体')
plt.ylabel('合計特殊出生率（仮想値）')
plt.grid(lw=0.5)
plt.show()

【実行結果】
合計特殊出生率が、施策後に「増えた」・「減った」・「横ばい」のケースが混在しています。

政策前後の合計特殊出生率の平均値の差が統計的に有意かどうかを対応のある２標本の平均のｔ検定（両側検定）で確かめます。
scipy.stats の ttest_1samp を利用します。

### 対応のある２標本の平均の差のｔ検定
result = stats.ttest_1samp(df5['差'], popmean=0, alternative='two-sided')
print(f't値: {result.statistic:.6f}')
print(f'p値: {result.pvalue:.6f}')
print('信頼区間 : ', result.confidence_interval(confidence_level=0.95))

【実行結果】
$${p}$$値は$${0.023441}$$であり、$${5\%}$$水準で平均値の差は有意である、との結果になりました。

ちなみに、pingouin 利用コードは以下のようになります。

### 対応のある２標本の母平均の差の検定 pingouin
# paired: 対応のある⇒True

pg.ttest(df5['政策後'], df5['政策前'], paired=True, correction=False,
         alternative='two-sided')

【実行結果】

政策を実施した自治体だけの分析では、政策以外の要因の影響を排除しきれず、因果関係の分析としては不十分、だそうです。
そこで、政策導入自治体と未導入自治体を比較する分析に進みます。

■ 政策を導入した自治体と導入しなかった自治体を比較 p.227～
10の少子対策導入自治体と10の少子対策未導入自治体の合計特殊出生率を分析します。
対応のない２標本の平均のｔ検定を利用します。

対策を導入しなかった自治体のデータを追加的に読み込みます。

### データの読み込み 対策導入しなかった自治体

df6 = pd.read_csv('./data/07_04_child2.csv', index_col=0)
df6['差'] = df6['期間２'] - df6['期間１']
print('df6.shape:', df6.shape)
display(df6)

【実行結果】

要約統計量を算出します。

### 要約統計量の算出
display(df6.describe().round(2))

【実行結果】
政策未導入の自治体の合計特殊出生率の平均値は、期間１から期間２にかけて$${0.02}$$上昇しています。

可視化してみます。

### 可視化
sns.lineplot(data=df6.iloc[:, :-1].T)
plt.legend(bbox_to_anchor=(1, 1), title='自治体')
plt.ylabel('合計特殊出生率（仮想値）')
plt.grid(lw=0.5)
plt.show()

【実行結果】
対策未導入の自治体は、ほぼ変化なし、のように見えます。

政策導入自治体の政策前後の「差」と政策未導入自治体の期間１・２の「差」について、分散が等しくないと仮定した２標本の平均の差の検定を実施します。
なお、テキストが用いるEXCELの$${p}$$値は「自由度を整数化して」求めており、一方で、Python等で算出するWelchの方法は「小数点の自由度」を用いているので、この記事の結果はテキストと異なることをあらかじめご了承くださいませ。

片側検定です。
scipy.stats の ttest_ind を利用します。

### 対応のない２標本の平均のｔ検定 等分散の仮定を置かない ★書籍の結果と異なる
# equal_var: 等分散を仮定しない⇒False,
# alternative: 差が大きい⇒greater
### 結論：有意水準5%で施策の効果あり

result =  stats.ttest_ind(df5['差'], df6['差'], equal_var=False,
                          alternative='greater')
print(f't値: {result.statistic:.6f}')
print(f'p値: {result.pvalue:.6f}')
print(f'信頼区間: {result.confidence_interval(confidence_level=0.95)}')

【実行結果】
$${p}$$値は$${0.025216}$$であり、$${5\%}$$水準で平均値の差は有意である、との結果になりました。
つまり、政策効果があった、ということです。

pingouin で計算してみます。

### 対応のある２標本の平均の差のｔ検定 pingouin ★書籍の結果と異なる
# correction: 等分散を仮定しない⇒True,
# alternative: 差が大きい⇒greater

pg.ttest(df5['差'], df6['差'], correction=True, alternative='greater')

【実行結果】
scipy.stats と同じ結果です。

今回の写経は以上です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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

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