「回帰分析から学ぶ計量経済学」をPythonで写経 ～ 第７章「反実仮想の世界」①相関係数の問題点、ランダム化比較実験

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

書籍の著者　山澤成康 先生

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

第７章は反実仮想、因果推論に取り組みます。
今回のテーマは相関係数の問題点とランダム化比較実験（RCT）です。
では書籍を開いて回帰分析の旅に出発です🚀

はじめに

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

このシリーズは書籍「回帰分析から学ぶ計量経済学: 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.1 相関係数の問題点

テキストは相関係数で問題になりうる３つのことを指摘しています。

1. 相関係数は２つの変数の関係の強さを表すものだが、どちらが原因でどちらが結果か分からない。

2. 相関関係の裏側には第３の要因が関係している場合がある。

3. 逆の因果関係の場合がある。

この節では上記の 2. 「第３の要因が関係している場合」と 3. 「逆の因果関係」について事例を見ていきます。

■ 第３の要因が関係している場合 p.215
原因にも結果にも影響を与える第３の変数があることを交絡と呼び、この第３の変数を交絡変数、交絡因子と呼ぶそうです。
テキストは、体重と年収の散布図を題材にして、第３の変数「年齢」の存在を浮き彫りにします。

データを読み込みます。

### データの読み込み
# 2020年のデータ 体重: kg, 年収: 万円

df1 = pd.read_csv('./data/07_01_weight_wage.csv')
df1 = df1.drop(df1.index[0])   # 20～24歳を削除
df1['年齢'] = list(range(27, 58, 5))
print('df1.shape:', df1.shape)
display(df1)

【実行結果】

体重と年収の回帰直線付き散布図を描画しましょう。

### 可視化 体重と年収

# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 4))
# 回帰直線付き散布図の描画
sns.regplot(data=df1, x='体重 [kg]', y='年収 [万円]', ci=None,
            line_kws={'color': 'tomato'}, ax=ax)
# 修飾
ax.set(ylim=(0, 1200))
ax.grid(lw=0.5)
plt.show()

【実行結果】
体重が増えると年収が増加する関係があるように見えますが。。。
なんだか変です。疑似相関かもしれません。

３つの変数の相関係数を確認します。

### 相関係数
display(df1.corr(numeric_only=True).round(3))

【実行結果】
体重と年収の相関係数は$${0.907}$$。強い正の相関があるように見えます。
ただ、体重と年齢の相関係数が$${0.832}$$、年齢と年収の相関係数が$${0.984}$$です。
年齢と年収は因果関係がありそうな感じがします。

年齢と年収の回帰直線付き散布図を描画しましょう。

### 可視化 年齢と年収

# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 4))
# 回帰直線付き散布図の描画
sns.regplot(data=df1, x='年齢', y='年収 [万円]', ci=None,
            line_kws={'color': 'tomato'}, ax=ax)
# 修飾
ax.set(ylim=(0, 1200))
ax.set_xticks(ticks=df1['年齢'], labels=df1['年齢層'])
ax.grid(lw=0.5)
plt.show()

【実行結果】
年功序列が残っている場合、あるいは、年令に応じてスキルが高くなる場合、こちらの相関関係の方がしっくりきます。
年齢が原因、年収が結果、のように解釈できるかもしれませんね。

■ 逆の因果関係 p.216
テキストは「出火件数」と「消防職員数」の散布図を描いて、原因と結果を探ろうとしています。

データを読み込みます。

### データの読み込み
# 2020年のデータ

df2 = pd.read_csv('./data/07_02_fire.csv')
print('df2.shape:', df2.shape)
display(df2.head())

【実行結果】
都道府県別の消防職員数と出火件数です。

相関係数を確認します。

### 相関係数
display(df2.corr(numeric_only=True).round(3))

【実行結果】
消防職員数と出火件数との間には、とても強い正の相関関係があります。

回帰直線付き散布図を描いてみましょう。

### 可視化 消防職員と出火件数

# 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 3))
# 消防職員と出火件数の回帰直線付き散布図の描画
sns.regplot(data=df2, x='消防職員数【人】', y='出火件数【件】', ci=None,
            line_kws={'color': 'tab:red', 'lw': 1.5},
            scatter_kws={'alpha': 0.5}, ax=ax1)
# 出火件数と消防職員の回帰直線付き散布図の描画
sns.regplot(data=df2, x='出火件数【件】', y='消防職員数【人】', ci=None,
            line_kws={'color': 'tab:red', 'lw': 1.5},
            scatter_kws={'alpha': 0.5}, ax=ax2)
# 修飾
ax1.grid(lw=0.5)
ax2.grid(lw=0.5)
fig.suptitle('因果関係あり？なし？')
plt.tight_layout()
plt.show()

【実行結果】

消防職員が多いから出火件数が多いのでしょうか？
それとも、出火件数が多いから消防職員が多いのでしょうか？
後者のほうが理にかなっている感じがいたします。
散布図で相関が見られたとして、原因・結果の関係を探る際は、注意が必要そうです。

7.5 RCT（ランダム化比較実験）

テキストは、政策効果を行った場合と行わなかった場合を実験によって試すことができる場合、ランダム化比較実験（RCT）が採られるとし、「ゴミ収集が終了した看板を置いて不法投棄がどのくらい減るか」を調べた仮想データを分析しています。
２０カ所のごみ収集場について、介入群（看板設置：対策あり）と対照群（看板非設置：対策なし）をランダムに割り付けして不法投棄率を算出し、平均値の差のｔ検定で不法投棄減少の効果を確かめます。

データを読み込みます。

### データの読み込み

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

【実行結果】
対策あり１０か所、対策なし１０か所をランダムに割り付けしています。
値は不法投棄率です。

要約統計量を確認します。

### 要約統計量の表示
display(df3.describe().T.round(3))

【実行結果】
対策ありの不法投棄率の平均値は$${19.8\%}$$、対策なしは$${23.2\%}$$。
不法投棄率の平均値の差は$${3.4}$$ポイントです。

不法投棄率の平均値を可視化しておきましょう。

### 可視化
plt.figure(figsize=(6, 4))
sns.barplot(data=df3, orient='y', alpha=0.7)
plt.xlabel('不法投棄率')
plt.grid(lw=0.7, axis='x');

【実行結果】

平均値の差のｔ検定を実行します。
scipy.stats の ttest_ind を利用して、片側検定（下側）と両側検定を実行します。

こちらは片側検定（下側）です。

### 対応のない２標本の母平均の差の検定 scipy.stats
# equal_var: 等分散を仮定しない⇒False,
# alternative: 片側検定（下側）⇒less

result =  stats.ttest_ind(df3['対策あり'], df3['対策なし'], equal_var=False,
                          alternative='less')
print(f't値: {result.statistic:.6f}')
print(f'p値: {result.pvalue:.6f}')
print(f'信頼区間: {result.confidence_interval(confidence_level=0.95)}')

【実行結果】
$${p}$$値は$${0.000012}$$であり、$${5\%}$$水準で平均値の差は有意です。

こちらは両側検定です。

### 対応のない２標本の母平均の差の検定 scipy.stats
# equal_var: 等分散を仮定しない⇒False,
# alternative: 両側検定⇒two-sided

result =  stats.ttest_ind(df3['対策あり'], df3['対策なし'], equal_var=False,
                          alternative='two-sided')
print(f't値: {result.statistic:.6f}')
print(f'p値: {result.pvalue:.6f}')
print(f'信頼区間: {result.confidence_interval(confidence_level=0.95)}')

【実行結果】
$${p}$$値は$${0.000023}$$であり、$${5\%}$$水準で平均値の差は有意です。

平均値の差のｔ検定によると、看板の設置によって不法投棄率が低下する効果があった、ということでしょう。

平均値の差のｔ検定を pingouin でやってみます。
もちろん上記の結果と同じになります。
こちらは片側検定（下側）です。

### 対応のない２標本の母平均の差の検定 pingouin
# correction: 等分散を仮定しない⇒True
# alternative: 片側・小さい⇒less

pg.ttest(df3['対策あり'], df3['対策なし'], correction=True, alternative='less')

【実行結果】

こちらは両側検定です。

### 対応のない２標本の母平均の差の検定 pingouin
# correction: 等分散を仮定しない⇒True
# alternative: 両側⇒two-sided

pg.ttest(df3['対策あり'], df3['対策なし'], correction=True, alternative='two-sided')

【実行結果】

RCTにより因果分析の仮定を満たしているものとして、回帰分析で因果効果を推定してみます。
書籍「Pythonで学ぶ効果検証入門」を参考にしました。

まずはデータの前処理です。データを縦持ちに変換します。

### 回帰分析による因果効果算出（因果分析の仮定を満たしていると仮定）
# 参考：Pythonで学ぶ効果検証入門

## データの前処理
df3reg = df3.melt(value_vars=['対策あり', '対策なし'], var_name='対策',
                  value_name='不法投棄率')
df3reg['対策'] = df3reg['対策'].map({'対策あり': 1, '対策なし': 0})
display(df3reg)

【実行結果】
対策＝0 は対策なし、対策＝1 は対策あり（看板設置）です。

回帰分析を実行します。
目的変数は不法投棄率、説明変数は対策（有無）です。

## 回帰分析の実行
result = smf.ols(formula='不法投棄率 ~ 対策', data=df3reg).fit()
display(result.summary())

【実行結果】
対策（有無）の係数$${-3.4}$$が効果の推定値です（$${5\%}$$水準で有意です）。
対策する（対策＝１）ことで、不法投棄率が$${3.4}$$ポイント低下する、という推定です。
不法投棄率の平均値の差$${3.4}$$ポイントと同じ結果になりました。

今回の写経は以上です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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

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