「回帰分析から学ぶ計量経済学」をPythonで写経 ～ 第１章「データを関連づける」②回帰分析をはじめる

第１章「データを関連づける」

書籍の著者　山澤成康 先生

この記事は、書籍「回帰分析から学ぶ計量経済学」第１章「データを関連づける」の Python写経活動 を取り扱います。

今回は 回帰分析と第１章の課題 に取り組みます！
では書籍を開いて回帰分析の旅に出発です🚀

はじめに

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

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

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

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

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

引用表記

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

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

第１章 データを関連づける

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

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

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

### インポート

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

# 機械学習(線形回帰・OLS)
from sklearn.linear_model import LinearRegression

# 描画
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'

1.5 回帰分析

■ 回帰分析の仕組み

テキストの数式をお借りして単回帰モデルを眺めてみます。

① 単回帰モデル：$${\alpha, \beta}$$は係数、$${u_i}$$は誤差項 または かく乱項
※テキストは$${Y}$$を結果または被説明変数、$${X}$$を原因または説明変数と呼んでいます。

② 理論値の計算：$${\hat{\alpha}, \hat{\beta}}$$はデータから推定した係数

③ 残差：実績値と理論値の差

テキスト p.38 の回帰分析の仕組みの散布図を描画します。

データを読み込みます。

### データの読み込み
df7 = pd.read_csv('./data/01_07_reg1.csv', index_col=0)
print('df7.shape:', df7.shape)
display(df7)

【実行結果】
サンプルサイズ 5 （5行）のデータです。

回帰分析を実行します。
scikit-learn の LinearRegression クラスを利用します。

### 回帰分析の実行

# 説明変数の準備 ※縦ベクトル化
X = df7['X'].values.reshape(-1, 1)
# 回帰分析の実行
reg = LinearRegression().fit(X, df7['Y'])
# 予測値の取得
y_pred = reg.predict(X)
# 回帰係数の表示
print(f'切片α: {reg.intercept_:.6f}')
print(f'傾きβ: {reg.coef_[0]:.6f}')

【実行結果】

散布図と回帰直線を描画します。

### 散布図と回帰直線の描画

# 描画領域の設定
plt.figure(figsize=(6, 4))
# 散布図の描画
plt.scatter(df7['X'], df7['Y'])
# 回帰直線の描画
plt.plot(df7['X'], y_pred, color='tomato')
# 残差の垂直点線の描画
for x, y1, y2 in zip(df7['X'], df7['Y'], y_pred):
    plt.plot([x, x], [y1, y2], color='black', ls=':')
# 修飾
plt.xlabel('X')
plt.ylabel('Y')
plt.grid(lw=0.5)
plt.show()

【実行結果】
青いデータ点は変数 X と変数 Y の実績値です。
赤い直線が回帰直線、黒い点線の長さは残差を示しています。

■ アルバイト収入と支出の単回帰分析 p.40
単回帰分析をやってみます。
アルバイト収入と支出額データを読み込みます。

### データの読み込み
df8 = pd.read_csv('./data/01_08_reg2.csv', index_col=0)
print('df8.shape:', df8.shape)
display(df8)

【実行結果】
5人分のデータです。
アルバイト収入を説明変数 X に、支出額を目的変数 Y にします。

散布図を描画します。

### 散布図の描画

# 描画領域の設定
plt.figure(figsize=(6, 4))
# 散布図の描画
plt.scatter(df8['アルバイト収入'], df8['支出額'])
# 修飾
plt.xlabel('アルバイト収入')
plt.ylabel('支出額')
plt.grid(lw=0.5)
plt.show()

【実行結果】
アルバイト収入と支出額との間には直線的な関係がありそうです。
相関が強そうな感じです。

回帰分析を実行して、支出額の理論値と残差を算出します。

### 回帰分析

## 回帰分析の準備
# 説明変数の準備 ※縦ベクトル化
X = df8['アルバイト収入'].values.reshape(-1, 1)

## 回帰分析
# 回帰分析の実行
reg = LinearRegression().fit(X, df8['支出額'])
# 予測値の取得
y_pred = reg.predict(X)

## データのまとめ
# データフレームに予測値と残差列を追加
df8_2 = df8.copy()
df8_2['支出額の理論値'] = y_pred
df8_2['残差'] = df8_2['支出額の理論値'] - df8_2['支出額']
# 結果表示
display(df8_2.round(1))

【実行結果】
テキスト p.40 の表と同じ結果になりました。

決定係数と自由度調整済み決定係数を算出します。
これらの指標はデータ（実績値）に対して推定した回帰式の当てはまり度合いを示します。
決定係数は scikit-learn の LinearRegression クラスの score メソッドで取得できます。
自由度調整済み決定係数はメソッド等で取得できないので、以下の公式で算出します。

### 決定係数の算出

# 準備 標本サイズnと説明変数の数kの算出
n, k = X.shape
# 決定係数の算出
R2 = reg.score(X, df8['支出額'])
# 自由度調整済み決定係数の算出
R2_adjusted = 1 - (1 - R2) * (n - 1) / (n - k - 1)  # 書籍の数式と異なる
# 結果表示
print('決定係数R²               = ', reg.score(X, df8['支出額']))
print('自由度調整済み決定係数R²* = ', R2_adjusted)

【実行結果】
ものすごく当てはまり度合いが高いです！

■ さまざまな決定係数とそのばらつき p.43
決定係数が$${0, 0.5, 0.9, 0.99}$$のデータを読み込んで、散布図と回帰直線を描画し、決定係数を算出します。

### データの読み込み
df9 = pd.read_csv('./data/01_09_R2.csv')
print('df9.shape:', df9.shape)
display(df9)

【実行結果】
サンプルサイズ 5 のデータです。

散布図と回帰直線を描画します。

### 散布図の描画 by matplotlib

## 設定と準備
# 説明変数 ※縦ベクトル化
X = df9['x'].values.reshape(-1, 1)
# 回帰直線用のXの値
X_pred = np.array([[X.min()], [X.max()]])

## 描画処理
# 描画領域の設定
fig, axes = plt.subplots(2, 2, figsize=(8, 6), sharex=True, sharey=True)

# ４つのyごとに回帰分析実行と散布図描画を繰り返し処理
for col, ax in zip(df9.columns[1:], axes.flat):
    
    ## 回帰分析
    # 回帰分析の実行
    reg = LinearRegression().fit(X, df9[col])
    # yの予測値の算出 ※回帰直線に使用
    y_pred = reg.predict(X_pred)
    # 決定係数の算出
    R2 = reg.score(X, df9[col])
    
    ## 描画
    # 散布図の描画
    ax.scatter(df9['x'], df9[col], s=60, ec='white', alpha=0.7)
    # 回帰直線の描画
    ax.plot(X_pred, y_pred, color='tab:red', ls='--')
    # 修飾
    ax.set(title=f'$R^2$ = {R2:.2f}', xlim=(0, 15), ylim=(0, 9),
           xticks=range(0, 17, 5))
    ax.grid(lw=0.5)

# 全体修飾
fig.supxlabel('x')
fig.supylabel('y')
plt.tight_layout()
plt.show()

【実行結果】
決定係数の値が大きくなるほど、データ点と回帰直線の当てはまり具合いが良くなっていることが分かります。

テキストによると「予測値と実績値の相関係数$${R}$$」がExcelの「重相関R」（重相関係数）であり、重「相関係数を二乗したものが決定係数$${R^2}$$」とのことです。

課題1 ヒストグラム

気象庁のWebサイトから1875年から2023年までの東京都の月次平均気温データを取得して、７月の平均気温のヒストグラムを描画します。
予め気象庁のWebサイトからダウンロードしたデータを読み込みます。

### データの読み込み
# https://www.data.jma.go.jp/stats/etrn/view/monthly_s3.php?prec_no=44&block_no=47662

df10 = pd.read_csv('./data/01_ex1.csv', index_col=0)
print('df10.shape:', df10.shape)
display(df10)

【実行結果】

７月の平均気温データの基本統計量を pandas の describe() メソッドで算出します。

### ７月の平均気温データの基本統計量
df10['7月'].describe().to_frame().round(1)

【実行結果】
平均値は 25.1 ℃、最大値は 28.7 ℃です。

７月の平均気温のヒストグラムを描画します。

### ７月の平均気温のヒストグラムを描画

# 描画領域の設定
plt.figure(figsize=(6, 4))
# ヒストグラムの描画
plt.hist(df10['7月'], bins=range(21, 31), ec='white', alpha=0.7)
# 修飾
plt.xlabel('気温 [℃]')
plt.ylabel('頻度 [年]')
plt.grid(lw=0.5)
plt.show()

【実行結果】

課題2 平均値と中央値

テキストの「ある国」の国民の年収データについて、平均値、中央値などの基本統計量を計算します。
年収データを読み込みます。

### データの読み込み
df11 = pd.read_csv('./data/01_ex2.csv')
print('df11.shape:', df11.shape)
display(df11.head())

【実行結果】
サンプルサイズ 21 のデータです。

テキストが示す基本統計量を算出する関数を定義します。

### 基本統計量算出関数の定義

def calc_stats(x):  # x: pandas.Series
    
    # 設定と準備
    x = x.dropna()  # １行目のNaNを削除
    result = []     # 結果を格納するリストの初期化
    N = len(x)      # 標本サイズ

    # 統計量の算出
    result.append(x.mean())                    # 平均値
    result.append(x.std(ddof=1) / np.sqrt(N))  # 標準誤差
    result.append(x.median())                  # 中央値
    result.append(x.mode()[0])                 # 最頻値
    result.append(x.std(ddof=1))               # 標準偏差
    result.append(x.var(ddof=1))               # 分散
    result.append(x.kurt())                    # 尖度
    result.append(x.skew())                    # 歪度
    result.append(x.max() - x.min())           # 範囲
    result.append(x.min())                     # 最小
    result.append(x.max())                     # 最大
    result.append(x.sum())                     # 合計
    result.append(N)                           # データの個数

    # ラベルの設定
    labels = ['平均値', '標準誤差', '中央値', '最頻値', '標準偏差', '分散',
              '尖度', '歪度', '範囲', '最小', '最大', '合計', 'データの個数']

    return result, labels

年収データの基本統計量を算出します。

### 基本統計量の算出 自作関数を使う

# 基本統計量の算出
results, labels = calc_stats(df11['年収（万円）'])
# データフレーム化
stats_df10 = pd.DataFrame(results, index=labels, columns=['統計量'])
# 結果表示
display(stats_df10.round(2))

【実行結果】
歪度が正の値であり、左側に偏っていることが推察され、平均値と中央値は乖離しています。

ヒストグラムで分布を確認しましょう。
seaborn の histplot でヒストグラムを描画します。

### ヒストグラムの描画
sns.histplot(data=df11, kde=True, ec='white')
plt.axvline(df11['年収（万円）'].mean(), color='tab:green', lw=3, label='平均値')
plt.axvline(df11['年収（万円）'].median(), color='tab:red', lw=3, label='中央値')
plt.grid(lw=0.5)
plt.legend();

【実行結果】
分布は左側に偏っていて、平均値よりも中央値の方が小さくなっている様子が分かりました。

課題3 相関係数

実質GDP、実質民間最終消費、実質民間設備投資、実質財貨サービスの輸出の相関係数を算出します。
データを読み込みます。

### データの読み込み
df12 = pd.read_csv('./data/01_ex3.csv', index_col=0)
print('df12.shape:', df12.shape)
display(df12.head())

【実行結果】
2018年度から2022年度までの年額（単位：兆円）のデータです。

相関係数を算出します。
pandas の corr() メソッドで相関行列を求めます。

### 相関行列の算出

# 相関行列のデータフレームを作成
corr_df12 = df12.corr()

# 相関係数の絶対値が0.9以上のセルに色を付けて表示
# エラーになる場合、.map を .applymap に変更してください
color = 'background-color: lightpink'
(corr_df12
 .style.map(lambda x: color if (abs(x)>=0.9) & (abs(x)!=1)  else '')
 .format(precision=3))

【実行結果】
色付けしたセルに該当する２変数は強い相関関係があります。

相関行列をヒートマップで可視化してみます。
seaborn の heatmap でヒートマップを描画します。

### ヒートマップで可視化
sns.heatmap(corr_df12, vmax=1, cmap='Greens', 
            annot=True, fmt='.3f', annot_kws={"size": 14});

【実行結果】
ヒートマップ形式で相関の強さを確認できます。

今回の写経は以上です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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

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