「回帰分析から学ぶ計量経済学」を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 statsmodels.api as sm
import statsmodels.formula.api as smf

# トービット、ヘーキットモデル
from py4etrics.tobit import Tobit

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

# ワーニング非表示
import warnings
warnings.simplefilter('ignore')

3.10 ロジットを使った分析

■ ロジットと対数オッズ p.106
目的変数$${Y}$$が$${0}$$か$${1}$$で表されるカテゴリ変数（二値変数）のケースを取り扱います。
二値の例としては、Yes／No、あり／なし、該当する／しない、合格／不合格、受注／失注などがあります。

0 or 1 の二値の目的変数$${Y}$$を持つモデルでは、$${Y}$$の最小値をゼロ、最大値を１とする変換が必要となり、前の記事で取り扱ったロジスティック曲線を当てはめる ロジット を検討します。
ロジットの推定式は以下のようになります。

ロジスティック回帰の数式ですね！

テキストは上の数式を変形して、目的変数$${Y=1}$$の確率にかかる対数オッズの数式を紹介しています。

$${P(Y=1)}$$は$${Y=1}$$の確率です。
右辺が単回帰のような1次式になり、とてもシンプルになりました。

■ 景気拡大期の確率 p.107
景気の拡大期を１、後退期を０とする目的変数をロジット（というかロジスティック回帰）で分析します。
テキストの景気データを読み込みます。

### データの読み込み

# CSVファイルの読み込み
df5 = pd.read_csv('./data/03_05_economic_expansion.csv', index_col=0)
# 四半期日付のインデックスを設定
df5.index = pd.date_range(start='1994-06-30', periods=df5.shape[0], freq='QE')
df5.index.name = '四半期'
# データフレームの表示
print('df5.shape:', df5.shape)
display(df5.head())

【実行結果】
1994年6月期から2023年3月期までの 116 四半期の季節調整済み実質GDP前期比と景気の状態（内閣府が定める景気基準日付に基づく）のデータです。
景気の状態は、0：景気後退期、1：景気拡大期の二値変数です。

実質GDP前期比の時系列折れ線グラフを描画します。
景気拡大期を青色で示しています。

### 時系列プロットの描画

# 描画領域の設定
fig, ax = plt.subplots(figsize=(8, 3))
# GDP前期比の時系列プロットの描画
ax.plot(df5['GDP前期比'], color='royalblue')
# 景気拡大期（景気の状態＝１）を塗りつぶし
ax.fill_between(df5.index, -25, 10, where=df5['景気の状態'] == 1,
                color='lightblue', alpha=0.4, label='景気拡大期')
# 修飾
ax.set(ylim=(-25, 10), xlabel='時系列 [四半期]', ylabel='GDP前期比')
plt.grid(lw=0.5)
plt.legend()
plt.show()

【実行結果】

ロジスティック回帰を実行します。
テキストはツール gretl を用いていますが、この記事では statsmodels の logit を利用します。

### ロジスティック回帰の実行
result = smf.logit(formula='景気の状態 ~ GDP前期比', data=df5).fit()
display(result.summary())

【実行結果】

～ ちょいと寄り道 ～
Python の機械学習ライブラリ scikit-learn でロジスティック回帰を実践してみましょう！
係数が statsmodels の logit と一致するように引数を設定しています。

### （参考）scikit-learnのロジスティック回帰

# 追加インポート
from sklearn.linear_model import LogisticRegression

# モデルの学習 ※penalty:正則化項なし、tol:解を求める際の精度をlogitと近くする
clf = LogisticRegression(penalty=None, tol=1e-6)
clf.fit(X=df5[['GDP前期比']], y=df5['景気の状態'])

# 係数の表示
print('Intercept:', clf.intercept_[0])
print('GDP前期比:', clf.coef_[0][0])

【実行結果】
statsmodels logit の係数の値とほぼ同じになりました（当たり前？）

では statsmodels の分析結果に戻ります。
推定したモデルで 116 四半期分の予測値を算出し、テキストの「景気拡大期の確率」チャートを描画しましょう。

### 可視化

# 描画領域の設定
plt.figure(figsize=(8, 3))
# 景気の状態の実績値の散布図を描画
plt.scatter(df5['GDP前期比'], df5['景気の状態'], ec='white', s=40, alpha=0.4,
            label='実績値')
# 景気の状態の予測値（確率値）の散布図を描画
plt.scatter(df5['GDP前期比'], result.predict(df5['GDP前期比']), color='tomato',
            ec='white', s=40, alpha=0.4, label='予測値')
# 景気の状態が0.5の水平線を描画
plt.axhline(0.5, color='black', ls='--', lw=0.5)
# 修飾
plt.xlabel('季節調整済み実質GDP前期比 [%]')
plt.ylabel('景気の状態')
plt.title('景気拡大期の確率\n景気後退期=0, 景気拡大期=1')
plt.grid(lw=0.5)
plt.legend()
plt.show()

【実行結果】
景気の状態$${=0.5}$$となるのは実質GDP前期比がおよそ$${-2\%}$$のときです。
実質GDP前期比が$${-2\%}$$を下回ると景気後退期になる可能性が高い、と読み取れるようです。

最後に、景気の状態の予測値（確率値）の時系列折れ線グラフを描画します。
景気拡大期を緑色で示しています。
確率値が 0.5 のラインに赤点線を描画します。

### 時系列プロットの描画

# 描画領域の設定
fig, ax = plt.subplots(figsize=(8, 3))
# GDP前期比の時系列プロットの描画
ax.plot(result.predict(df5['GDP前期比']), color='green')
# 景気拡大期（景気の状態＝１）を塗りつぶし
ax.fill_between(df5.index, -25, 10, where=df5['景気の状態'] == 1,
                color='lightgreen', alpha=0.2, label='景気拡大期')
# 景気の状態が0.5の水平線を描画
plt.axhline(0.5, color='tab:red', lw=0.5)
# 修飾
ax.set(ylim=(-0.1, 1.1), xlabel='時系列 [四半期]', ylabel='景気の状態（確率値）')
plt.grid(lw=0.5)
plt.legend()
plt.show()

【実行結果】

3.11 トービットとヘーキット

この節は難易度が高いです！

3.11.1 トービット

テキストは数学試験の点数を例にしてトービットモデルを検討します。
学生の点数は下限 0 点～上限 100 点の間に分布します。
打ち切り（censored）に該当するデータです。

① 簡単なテストのケース
試験が簡単すぎると 100 点の学生が増えます。
100点を超えて得点できないため、100点で打ち切られるのです。

② 難しいテストのケース
試験が難しすぎると 0 点の学生が増えます。
0点を下回って得点できないため、0点で打ち切られるのです。

■ 真の実力と試験の結果 p.108
テキストのデータを読み込みます。

### データの読み込み

# CSVファイルの読み込み
df6 = pd.read_csv('./data/03_06_tobit.csv')
# データフレームの表示
print('df6.shape:', df6.shape)
display(df6)

【実行結果】
19名の数学の真の学力値と、難易度別の数学試験の得点です。

真の学力と３つの難易度の散布図を描画します。

### 散布図による可視化
# 簡単すぎる試験：上限値100で打ち切り、難しすぎる試験：下限値0で打ち切り

# 描画領域の設定
fig, axes = plt.subplots(1, 3, figsize=(10, 3))
# 3つのグラフごとに散布図描画を繰り返し処理
for col, ax in zip(df6.columns[1:], axes):
    # 真の実力と実際の点数の散布図の描画
    ax.scatter(df6['真の学力'], df6[col])
    # 修飾
    ax.set(xlabel='数学の真の実力', ylabel='数学の点数', title=col,
           xlim=(-5, 105), ylim=(-5, 105))
    ax.grid(lw=0.5)
# 全体調整
plt.tight_layout()

【実行結果】
真ん中の簡単すぎる試験は 100 点の天井があるようです。
右側の難しすぎる試験は 0 点の床があるようです。

点数のヒストグラムを比べてみましょう。

### ヒストグラムによる可視化
# 簡単すぎる試験：上限値100で打ち切り、難しすぎる試験：下限値0で打ち切り

# 描画領域の設定
fig, axes = plt.subplots(1, 3, figsize=(10, 3))
# 3つのグラフごとに散布図描画を繰り返し処理
for col, ax in zip(df6.columns[1:], axes):
    # 点数のヒストグラムの描画
    ax.hist(df6[col], ec='white', alpha=0.7)
    # 修飾
    ax.set(xlabel='数学の点数', ylabel='人数', title=col, xlim=(-5, 105))
    ax.grid(lw=0.5)
# 全体調整
plt.tight_layout()

【実行結果】
真ん中の簡単すぎる試験は 100 点の人数が多いです。
右側の難しすぎる試験は 0 点の人数が多いです。

■ 試験の推計結果（通常の回帰分析） p.109
通常の回帰分析で３つの試験の回帰分析を実行します。
テキストによると、データに打ち切りが生じていることを「考慮せずに回帰分析を行うと、正しい係数が推定できません」とのことです。

まずは通常の試験です。

### 通常の回帰分析 通常の試験
result_normal = smf.ols(formula='通常の試験 ~ 真の学力', data=df6).fit()
display(result_normal.summary().tables[1])

【実行結果】

続いて簡単すぎる試験です。

### 通常の回帰分析 簡単すぎる試験
result_easy = smf.ols(formula='簡単すぎる試験 ~ 真の学力', data=df6).fit()
display(result_easy.summary().tables[1])

【実行結果】

最後に難しすぎる試験です。

### 通常の回帰分析 難しすぎる試験
result_difficult = smf.ols(formula='難しすぎる試験 ~ 真の学力', data=df6).fit()
display(result_difficult.summary().tables[1])

【実行結果】

通常の試験の「真の学力」の係数が 0.87 なのに対して、簡単すぎる試験は 0.66、難しすぎる試験は 0.68 という低い値になっています。
これが「正しい係数が推定できません」だそうです。

■ トービットモデル
テキストによると「打ち切りデータを推定する場合に使うのがトービットになります」。
そこでトービットモデルを構築して分析しましょう！
py4etrics ライブラリの Tobit を利用します。

◆ 簡単すぎる試験のケース ◆
簡単すぎる試験からトライします！
上限（right）を 100 で打ち切ります。
また Tobit の引数用に、データ点ごとに「上限に該当する／しない」を指定する変数 censor を作っています。

### トービットモデルに対応した回帰 簡単すぎる試験 上限100点で打ち切り

## 設定
# モデル式の設定
formula = '簡単すぎる試験 ~ 1 + 真の学力'
# 上限値の設定
right = 100
# 打ち切りの設定(Series型)：目的変数の各値が上限値に該当するとき=1, 上限でないとき=0
censor = df6['簡単すぎる試験'].apply(lambda x: 1 if x == right else 0)

## 回帰の実行
result_easy_tobit = Tobit.from_formula(
    formula=formula, cens=censor, right=right, data=df6
).fit()
display(result_easy_tobit.summary())

【実行結果】
真の実力の係数は 0.88 。
通常の試験の係数 0.87 と近くなりました。

散布図の上に通常の回帰分析とトービットモデルの回帰直線を可視化してみましょう。

### 簡単すぎる試験の線形回帰モデルとTobitモデルの予測値の可視化

# 観測値の散布図の描画
plt.scatter(data=df6, x='真の学力', y='簡単すぎる試験', color='tab:blue',
            s=70, alpha=0.5, label='観測値')
# 線形回帰の当てはめ値の描画
plt.plot(df6['真の学力'], result_easy.fittedvalues, color='tab:green',
         label='線形回帰モデル')
# Tobitモデルの当てはめ値の描画
plt.plot(df6['真の学力'], result_easy_tobit.fittedvalues, color='tab:orange',
         label='Tobitモデル(回帰直線)')
# Tobitモデルの当てはめ値の描画
plt.plot(df6['真の学力'], result_easy_tobit.fitted_endog, color='tab:red',
         ls='--', label='Tobitモデル(上限考慮)')
# 修飾
plt.xlabel('数学の真の学力')
plt.ylabel('数学の得点')
plt.title('簡単すぎる試験の線形回帰モデルとTobitモデルの予測値の比較')
plt.grid(lw=0.5)
plt.legend();

【実行結果】
緑色が通常の回帰分析の回帰直線、オレンジ色がトービットモデルの回帰直線、赤い点線がトービットモデルの上限を考慮した曲線です。
通常の回帰分析の方は 100 点の影響を受けて傾きが小さくなっています。
トービットモデルの方は 100 点以外の点をうまく捉えているようですね！
打ち切りの100点の影響を抑えています。

◆ 難しすぎる試験のケース ◆
続いて難しすぎる試験のトービットモデルです。
下限（left）を 0 で打ち切ります。
また Tobit の引数用に、データ点ごとに「下限に該当する／しない」を指定する変数 censor を作っています。

### トービットモデルに対応した回帰 難しすぎる試験 下限0点で打ち切り

## 設定
# モデル式の設定
formula = '難しすぎる試験 ~ 1 + 真の学力'
# 下限値の設定
left = 0
# 打ち切りの設定：目的変数の各値が下限値に該当するとき=-1, 下限でないとき=0
censor = df6['難しすぎる試験'].apply(lambda x: -1 if x == left else 0)

## 回帰の実行
result_difficult_tobit = Tobit.from_formula(
    formula=formula, cens=censor, left=left, data=df6
).fit()
display(result_difficult_tobit.summary())

【実行結果】
真の実力の係数は 0.78 。
通常の試験の係数 0.87 と近くなりました。

散布図の上に通常の回帰分析とトービットモデルの回帰直線を可視化してみましょう。

### 難しすぎる試験の線形回帰モデルとTobitモデルの予測値の可視化

# 観測値の散布図の描画
plt.scatter(data=df6, x='真の学力', y='難しすぎる試験', color='tab:blue',
            s=70, alpha=0.5, label='観測値')
# 線形回帰の当てはめ値の描画
plt.plot(df6['真の学力'], result_difficult.fittedvalues, color='tab:green',
         label='線形回帰モデル')
# Tobitモデルの当てはめ値の描画
plt.plot(df6['真の学力'], result_difficult_tobit.fittedvalues, color='tab:orange',
         label='Tobitモデル(回帰直線)')
# Tobitモデルの当てはめ値の描画
plt.plot(df6['真の学力'], result_difficult_tobit.fitted_endog, color='tab:red',
         ls='--', label='Tobitモデル(下限考慮)')
# 修飾
plt.xlabel('数学の真の学力')
plt.ylabel('数学の得点')
plt.title('難しすぎる試験の線形回帰モデルとTobitモデルの予測値の比較')
plt.grid(lw=0.5)
plt.legend();

【実行結果】
緑色が通常の回帰分析の回帰直線、オレンジ色がトービットモデルの回帰直線、赤い点線がトービットモデルの上限を考慮した曲線です。
通常の回帰分析の方は 0 点の影響を受けて傾きが小さくなっています。
トービットモデルの方が 0 点以外の点をうまく捉えているようですね！
打ち切りの0点の影響を抑えています。

3.11.2 ヘーキット

テキストによると、へ－キットはサンプルにセレクション・バイアスがある場合に使う、ようです。

セレクション・バイアスの理解を深めるために、次の２つのブログ情報を引用させていただき、ヘーキットモデルを検討します。
ありがとうございます！

１つ目
[Rによるデータ分析入門] 離散選択モデルの様々（３）：ヘーキットモデル

２つ目
Pythonを使った計量経済分析－制限従属変数モデル

■ セレクション・バイアスのあるデータ p.110
１つ目のブログの掲載情報をお借りしてセレクション・バイアスのあるデータの回帰分析を確認します。

データの内容は、健康状態と賃金です。
健康状態が良好であり働ける人は企業が提示した賃金で働くとしています。
実際に観測できるのはこの「実際に働いている人」の賃金データです。
一方で健康状態が良くなくて企業が提示した賃金では働けない人がいるとします。
この「働けない人」の賃金データは観測されません。

健康状態と賃金を回帰分析で推定するときに、観測できる賃金データだけでは「健康良好な就労者」だけに偏ることになります。
この偏りがセレクション・バイアスに相当するようです。

以下のコードは、健康良好な就労者だけのセレクション・バイアスありデータによる回帰分析と、健康状態が良くない未就労者の企業提示賃金を含むセレクション・バイアスなしデータによる回帰分析を可視化して比較するものです。

### バイアスあり／なしの回帰分析

## データの作成
x_data = np.array([25, 60, 101, 140, 160, 171, 180, 190, 218, 248, 285, 306, 320])
y_data = np.array([800, 890, 950, 1030, 1150, 1240, 1330, 1503, 1600, 1800, 2030,
                   2150, 2340])
wage_df = pd.DataFrame({'健康状態': x_data, '賃金': y_data})

## 回帰分析
# 健康状態が良好のデータのみで回帰分析
result_observed = smf.ols(formula='賃金~健康状態', data=wage_df.iloc[3:, :]).fit()
# 健康状態が不良を含む全データで回帰分析
result_all = smf.ols(formula='賃金~健康状態', data=wage_df).fit()

## 描画
# 設定
x_line = np.linspace(x_data.min(), x_data.max(), 1000)  # 回帰直線用のx軸の値
# 散布図の描画：健康状態が良好で就労しているデータ
plt.scatter(x_data[3:], y_data[3:], color='tab:blue', s=70, alpha=0.7,
            label='健康良好・就労の実際賃金')
# 散布図の描画：健康状態が不良で就労していないデータ
plt.scatter(x_data[:3], y_data[:3], color='tab:red', s=70, alpha=0.7,
            label='健康不良・未就労の提示賃金')
# 回帰直線の描画：健康状態が良好のデータのみで回帰分析
plt.plot(x_line, result_observed.predict(exog=dict(健康状態=x_line)),
         color='tab:blue', label='健康良好・就労のみの回帰直線')
# 回帰直線の描画：健康状態が不良のデータを含む全データで回帰分析
plt.plot(x_line, result_all.predict(exog=dict(健康状態=x_line)),
         color='dimgray',ls='--', label='健康不良・未就労を含めた回帰直線')
# 修飾
plt.xlabel('健康状態', fontsize=14)
plt.ylabel('賃金', fontsize=14)
plt.xlim(0, 350)
plt.ylim(0, 2500)
plt.grid(lw=0.5)
plt.legend();

【実行結果】
セレクション・バイアスありの回帰直線（青線）とセレクション・バイアスなしの回帰直線（グレー点線）は傾きが異なります。
セレクション・バイアスがある場合、通常の回帰分析では適切な推定ができないことが分かりました。

■ 逆ミルズ比 p.111
ヘーキットモデルは「逆ミルズ比」という指標をもちいて、セレクション・バイアスに対応するようです。
テキストの数式をお借りします。

【第１段階】働くかどうかをプロビットモデルで推計して逆ミルズ比を算出

おそらく、$${Pr(M_i=1)}$$は働く確率、$${F(\cdot)}$$は標準正規分布の累積分布関数、だと思われます。

【第２段階】説明変数に逆ミルズ比を含めた線形回帰モデルで推定

テキストの図表「逆ミルズ比の図解」を描いてみます。

### 逆ミルズ比の計算

## データの作成
# 確率変数の値
x = np.linspace(-3, 3, 1001)
# 標準正規分布の確率密度関数
prob_density = stats.norm.pdf(x=x, loc=0, scale=1)
# 標準正規分布の累積分布関数
cum_dist = stats.norm.cdf(x=x, loc=0, scale=1)

## 逆ミルズ比の計算
inv_mills_ratio = prob_density / cum_dist

## 可視化
# 描画領域の設定 twinxで左右２軸化
fig, ax = plt.subplots(figsize=(6, 3))
twin_ax = ax.twinx()
# 累積分布関数の描画
ax.plot(x, cum_dist, color='tab:blue', label='働く可能性')
# 逆ミルズ比の曲線の描画
twin_ax.plot(x, inv_mills_ratio, color='tab:red', label='逆ミルズ比')
# 凡例
handles1, labels1 = ax.get_legend_handles_labels()
handles2, labels2 = twin_ax.get_legend_handles_labels()
ax.legend(handles=handles1 + handles2, labels=labels1 + labels2,
          bbox_to_anchor=(0.5, 1), ncols=2, loc='lower center')
# 修飾
ax.set(ylabel='働く可能性')
twin_ax.set(ylabel='逆ミルズ比')
ax.grid(lw=0.5)

【実行結果】
テキストのとおり、逆ミルズ比は「働く可能性が高いときにはゼロに近く」なり、「$${Y_i}$$には$${X_i}$$のみが反映」されます（逆ミルズ比の項は反映されにくいとのこと）。
一方で、「働く可能性が低い場合には逆ミルズ比が大きくなり、$${c \lambda_i}$$分だけ$${Y_i}$$が補正されます」になります。

■ ヘーキットモデルの例題（mrozデータセット）
py4etrics ライブラリの Heckit を用いてヘーキットモデルを実践します。
Heckit の使い方は２つ目のブログに準拠いたしました。
mrozデータセットは１つ目・２つ目の両方のブログで分析しています。
この記事では１つ目のブログに準拠して分析します。

mrozデータセットは女性の労働参加に関するデータセットだそうです。
今回のモデルで使用する変数を１つ目のブログの説明をお借りして表にします。

Pythonでヘーキットモデルを構築しましょう！

追加ライブラリをインポートします。

### インポート
import wooldridge                     # mrozデータセットを含むライブラリ
from py4etrics.heckit import Heckit   # ヘーキット

続いてmrozデータセットを読み込みます。

### データセットの読み込み
mroz = wooldridge.data('mroz')
print('mroz.shape:', mroz.shape)
display(mroz)

【実行結果】
サンプルサイズ 753 のデータセットです。

２つ目のブログにならって、労働参加有無と対数賃金の関係を確認します。

### 労働参加有無inlfと対数賃金lwageの関係の確認
# inlfは1975年の既婚女性の労働参加有無であり、参加した場合は1，参加しなかった場合は0
# 参加有無による対数賃金lwageの値の違いを確認する。
# ・inlf=1（参加した）の場合，lwageはNaNではない。
# ・inlf=0（参加しなかった）の場合，lwageはNaNである。

print('inlf==1の個数: ', mroz.query('inlf == 1').shape[0])
print('inlf==1の場合のlwageのNaNの個数: ',
      mroz.query('inlf == 1')['lwage'].isna().sum())
print('inlf==0の個数: ', mroz.query('inlf == 0').shape[0])
print('inlf==0の場合のlwageのNaNの個数: ',
      mroz.query('inlf == 0')['lwage'].isna().sum())

【実行結果】
労働参加している場合、対数賃金には欠損値が無く、値が入っています。
労働参加していない場合、対数賃金は欠損値となっています。
Heckit は、目的変数：対数賃金が欠損値かどうかで、労働参加している／していないを把握するようです。

分析に使用する変数を選択します。

### 分析に使用する変数の作成

# 第２段階の目的変数の作成：対数賃金
endog = mroz.loc[:,'lwage']

# 第２段階の説明変数の作成：教育年数、経験年数、経験年数の二乗
exog = mroz.loc[:,['educ', 'exper', 'expersq']]
exog['Intercept'] = 1.0  # 定数項

# 第１段階の説明変数の作成：
# 教育年数、経験年数、経験年数の二乗、妻の所得、6歳未満子ども有無、6歳以上子ども有無
exog_select = mroz.loc[:,['educ', 'exper', 'expersq','nwifeinc', 'kidslt6',
                          'kidsge6']]
exog_select['Intercept'] = 1.0  # 定数項

２つ目のブログにならって、上のコードで作成した変数の行数が同一になっていることを確認します。

### endog、exog、exog_selectの行数が同じになっていることを確認
print('endog.shape      : ', endog.shape)        # 第２段階の目的変数
print('exog.shape       : ', exog.shape)         # 第２段階の説明変数
print('exog_select.shape: ', exog_select.shape)  # 第１段階の説明変数

【実行結果】
データセットの行数と同じ 753 となっています。

いよいよ分析に入ります。
注目する変数は 教育年数 educ です。
まず、通常の回帰分析を実行します。

### 通常の回帰分析（最小二乗法）による推定：heckitモデルの第２段階の変数を使用
result_ols = sm.OLS(endog, exog, missing='drop').fit()
print(result_ols.summary())

【実行結果】
教育年数 educ の係数は$${0.1075}$$です（有意です）。

ヘーキットモデルで推定します。

### ヘーキットモデルによる推定の実行
# 推定を行う際，fit()にオプションを追加し不均一分散頑健標準誤差を指定する。
# cov_type_1：第１段階推定でのオプション
# cov_type_2：第２段階推定でのオプション

result_heckit = Heckit(endog, exog, exog_select).fit(cov_type_2='HC1')
print(result_heckit.summary())

【実行結果】
係数表記ゾーンの内、上側が第２段階、下側が第１段階です。

今回の写経は以上です。
難解でしたね・・・

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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

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