上級マクロ経済学

第9回 Pythonによる回帰分析の実装

Author

荻巣嘉高

1 Pythonによる回帰分析

単回帰分析

y_i = \beta_0 + \beta_1 x_i + u_i

をモデルとして、回帰分析をする。

  • 推定値の\beta_1に関心がある場合が多い。

データの読み込み

  • データを用意してください。 wage_and_productivity.xlsx
import numpy as np
import pandas as pd
import scipy.stats as stat
import matplotlib.pyplot as plt
import statsmodels.api as sm # statsmodelsを使うためインポート

# データの読み込み
df = pd.read_excel("wage_and_productivity.xlsx", sheet_name=0)

df.head() # 最初の方のデータを表示
2020年 労働生産性（円/時間） 賃金所得（千円）
0 北海道 8997.702921 4232.5
1 青森 8946.026251 3667.9
2 岩手 8336.011457 3790.4
3 宮城 8875.070576 4459.4
4 秋田 8700.560140 3758.0

変数の割り当て

## 変数を設定

x = df["労働生産性（円/時間）"].values # 労働生産性の列の値をxとする
#x = df.iloc[:,1] # こっちでもOK

X = sm.add_constant(x) # 回帰式に定数項部分を足す

y = df["賃金所得（千円）"].values # 賃金所得をyにする


print(X[:5]) # Xの最初の方を表示する
[[1.00000000e+00 8.99770292e+03]
 [1.00000000e+00 8.94602625e+03]
 [1.00000000e+00 8.33601146e+03]
 [1.00000000e+00 8.87507058e+03]
 [1.00000000e+00 8.70056014e+03]]

statsmoelsでの回帰

statsmodelsでOLSする場合は、sm.OLS()とfit()がセット。

## statsmodelで回帰

## OLS()とfit()をセットで結果を出す
model = sm.OLS(y,X) # yをXにOLSで回帰
result = model.fit() # 先のmodelに.fit()をつけると、結果の計算がされる

#result = sm.OLS(y,X).fit() # こっちでもOK

result.summary() # 結果を表示する。
OLS Regression Results
Dep. Variable: y R-squared: 0.181
Model: OLS Adj. R-squared: 0.163
Method: Least Squares F-statistic: 9.978
Date: Sat, 13 Sep 2025 Prob (F-statistic): 0.00283
Time: 13:13:34 Log-Likelihood: -350.07
No. Observations: 47 AIC: 704.1
Df Residuals: 45 BIC: 707.8
Df Model: 1
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
const 2681.0906 550.477 4.870 0.000 1572.373 3789.808
x1 0.1830 0.058 3.159 0.003 0.066 0.300
Omnibus: 13.961 Durbin-Watson: 1.148
Prob(Omnibus): 0.001 Jarque-Bera (JB): 15.953
Skew: 1.089 Prob(JB): 0.000343
Kurtosis: 4.844 Cond. No. 8.45e+04


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 8.45e+04. This might indicate that there are
strong multicollinearity or other numerical problems.

結果の読み方

対応する値 回帰係数 標準誤差 t統計量 p値
表示 coef std err t P>|t|
const 2681.1 550.5 4.870 0.000
x1 0.183 0.058 3.159 0.003
  • ここのt統計量は、H_0:~ \beta = 0と仮定した時のt統計量が計算されています。
  • p値は、真の値として\beta_i=0を仮定した時に、推定された\hat{\beta}_iの大きさ（つまり|\hat{\beta}_i|）以上の値が現れる確率を表しています。
    • ただし、ここでの確率計算ではt分布を用いていることに注意が必要です。
    • Nが十分に大きければt分布は標準正規分布に近似されるので、あまり気にする必要はありません。

p値

  • \beta_i=0の帰無仮説を設定している場合、p値が0.05以下なら、5%の有意水準で帰無仮説が棄却される。

2025-09-13T13:13:34.505451 image/svg+xml Matplotlib v3.9.2, https://matplotlib.org/

個別に結果を表示

## statsmodelで回帰

## OLS()とfit()をセットで結果を出す
model = sm.OLS(y,X) # yをXにOLSで回帰
result = model.fit() # 先のmodelに.fit()をつけると、結果の計算がされる

#result = sm.OLS(y,X).fit() # こっちでもOK

beta = result.params # 回帰係数
SE = result.bse # 標準誤差
t = result.tvalues # t統計量

print("回帰係数は"+str(beta)+"です。")
print()

print("標準誤差は"+str(SE)+"です。")
print()

print("t統計量は"+str(t)+"です。")
回帰係数は[2.68109056e+03 1.82966021e-01]です。

標準誤差は[5.50477090e+02 5.79230982e-02]です。

t統計量は[4.87048528 3.15877477]です。

statsmoelsの分散

デフォルトの回帰分析では、

  • 誤差項が均一分散
    • u_iの分散がX_iの大きさによって変動しない。

であることを仮定して標準誤差を計算しています。

一方で、近年の分析では

  • 誤差項が不均一分散である場合にも対応できる標準誤差を計算することがほとんどです。
  • statsmodelでは、fit()のオプションにcov_typeがあり、これを"HC0"と指定すれば不均一分散にも対応できるような標準誤差が計算できます。
    • これを不均一分散に頑健な標準誤差（Whiteの標準誤差）と言います。
  • 基本はこちらを用いましょう。

不均一分散に頑健な標準誤差

## statsmodelで回帰

model = sm.OLS(y,X)
result0 = model.fit() # 均一分散を仮定
result1 = model.fit(cov_type="HC0") # 不均一分散に頑健


SE0 = result0.bse # 均一分散仮定のSE
SE1 = result1.bse # 不均一分散に頑健なSE

t0 = result0.tvalues
t1 = result1.tvalues


print("---均一分散を仮定すると---")
print("標準誤差は"+str(SE0))
print("t統計量は"+str(t0))
print()

print("---不均一分散に頑健に計算すると---")
print("標準誤差は"+str(SE1))
print("t統計量は"+str(t1))
---均一分散を仮定すると---
標準誤差は[5.50477090e+02 5.79230982e-02]
t統計量は[4.87048528 3.15877477]

---不均一分散に頑健に計算すると---
標準誤差は[5.20905463e+02 5.52152006e-02]
t統計量は[5.14698108 3.31368934]

対数をとった回帰

こういった推計をする際に、説明変数と被説明変数についてそれぞれ対数を取ってから回帰を行うことがあります。

\log Y = \beta_0 + \beta_1 \log X + u

  • このとき、\beta_1の値は、「Xが1%上昇した時に、Yが何%上昇するか」を表すと解釈されます。
    • なぜでしょうか？
  • 単位の異なる変数を用いて回帰分析を行うことも多い。
    • %で解釈ができるようになる対数を取った回帰は頻出します。

対数を取った回帰の実装

## 変数を設定

x = np.log(df["労働生産性（円/時間）"].values) # 労働生産性の列の値をxとする
#x = np.log(df.iloc[:,1].values) # こっちでもOK

X = sm.add_constant(x) # 回帰式に定数項部分を足す

y = np.log(df["賃金所得（千円）"].values) # 賃金所得をyにする


print(X[:5]) # Xの最初の方を表示する

## OLS()とfit()をセットで結果を出す
model = sm.OLS(y,X) # yをXにOLSで回帰
result = model.fit(cov_type="HC0") # 先のmodelに.fit()をつけると、結果の計算がされる
# cov_type="HC0"を指定

#result = sm.OLS(y,X).fit(cov_type="HC0") # こっちでもOK

result.summary() # 結果を表示する。
[[1.         9.10472459]
 [1.         9.09896472]
 [1.         9.02834014]
 [1.         9.09100157]
 [1.         9.07114269]]
OLS Regression Results
Dep. Variable: y R-squared: 0.210
Model: OLS Adj. R-squared: 0.193
Method: Least Squares F-statistic: 13.41
Date: Sat, 13 Sep 2025 Prob (F-statistic): 0.000656
Time: 13:13:34 Log-Likelihood: 46.323
No. Observations: 47 AIC: -88.65
Df Residuals: 45 BIC: -84.95
Df Model: 1
Covariance Type: HC0
coef std err z P>|z| [0.025 0.975]
const 4.5367 1.051 4.315 0.000 2.476 6.597
x1 0.4208 0.115 3.662 0.000 0.196 0.646
Omnibus: 7.225 Durbin-Watson: 1.200
Prob(Omnibus): 0.027 Jarque-Bera (JB): 6.106
Skew: 0.774 Prob(JB): 0.0472
Kurtosis: 3.849 Cond. No. 765.


Notes:
[1] Standard Errors are heteroscedasticity robust (HC0)

2 重回帰分析

重回帰分析のモデル

回帰分析の説明変数が複数ある場合、この回帰分析を重回帰分析という。

  • 一般に、n個の説明変数X_1, X_2, \cdots, X_nが存在している場合、モデルは次のように表される。 Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_n X_n + u

この\beta_iについても、前回までの内容と同様に係数の推計と仮説検定が可能。

pythonでの実装

  • wage_and_others.xlsxをDL、pythonで読み込みをしてください。
import numpy as np
import pandas as pd
import statsmodels.api as sm

df = pd.read_excel("wage_and_others.xlsx", index_col=0)
df.head()
賃金所得（千円） 労働生産性（円/時間） 人口密度 生産年齢人口比率
2020年
北海道 4232.5 8997.702921 66.6 57.20614
青森 3667.9 8946.026251 128.3 55.72851
岩手 3790.4 8336.011457 79.2 55.41224
宮城 4459.4 8875.070576 316.1 60.18364
秋田 3758.0 8700.560140 82.4 52.83574
  • 縦方向にサンプルを並べて、横方向に説明変数X1, X2, \cdots, X_nを並べる。

回帰モデルの設定

\begin{aligned} \log (賃金) = \beta_0 &+ \beta_1 \log(労働生産性) \\ &+ \beta_2 \log(人口密度) \\ &+ \beta_3 \log(生産年齢人口比率) + u \end{aligned}

で、帰無仮説H_0と対立仮説H_1をそれぞれ以下のように設定する。 H_0 :~ \beta_1 = 0 H_1 :~ \beta_1 \neq 0

データフレームでの説明変数設定

  • Pythonのstatsmodelsの説明変数では、ndarrayではなくpandasのDataFrameを指定することもできる。
  • 結果が読みやすくなるので、可能ならこちらを用いるのもおすすめ。

重回帰の実装

y = np.log(df.iloc[:,0].values) # 第0列目（賃金）の対数をyとする

x = np.log(df.iloc[:,1:]) # 1列目以降（労働生産性、人口密度、生産年齢人口比率）の対数をxとする。
X = sm.add_constant(x) # xに定数項を加える。

print("Xのチェック、最初の3行") 
print(X.head(3)) # データをチェックする。

result = sm.OLS(y,X).fit(cov_type="HC0") # OLSを実行

print(result.summary()) # 結果の表示
Xのチェック、最初の3行
       const  労働生産性（円/時間）      人口密度  生産年齢人口比率
2020年                                        
北海道      1.0     9.104725  4.198705  4.046661
青森       1.0     9.098965  4.854371  4.020492
岩手       1.0     9.028340  4.371976  4.014801
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.791
Model:                            OLS   Adj. R-squared:                  0.776
Method:                 Least Squares   F-statistic:                     121.2
Date:                Sat, 13 Sep 2025   Prob (F-statistic):           5.34e-21
Time:                        13:13:34   Log-Likelihood:                 77.558
No. Observations:                  47   AIC:                            -147.1
Df Residuals:                      43   BIC:                            -139.7
Df Model:                           3                                         
Covariance Type:                  HC0                                         
===============================================================================
                  coef    std err          z      P>|z|      [0.025      0.975]
-------------------------------------------------------------------------------
const           1.4936      0.766      1.949      0.051      -0.009       2.996
労働生産性（円/時間）     0.3491      0.078      4.494      0.000       0.197       0.501
人口密度            0.0406      0.011      3.696      0.000       0.019       0.062
生産年齢人口比率        0.8563      0.305      2.808      0.005       0.259       1.454
==============================================================================
Omnibus:                       23.205   Durbin-Watson:                   1.322
Prob(Omnibus):                  0.000   Jarque-Bera (JB):               44.407
Skew:                          -1.393   Prob(JB):                     2.28e-10
Kurtosis:                       6.862   Cond. No.                     1.94e+03
==============================================================================

Notes:
[1] Standard Errors are heteroscedasticity robust (HC0)
[2] The condition number is large, 1.94e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

結果の読み方

変数 係数 標準誤差 z値（t値） p値
表示 coef std err z P>|z|
const 1.4936 0.766 1.949 0.051
労働生産性（円/時間） 0.349 0.078 4.494 0.000
人口密度 0.0406 0.011 3.696 0.000
生産年齢人口比率 0.8563 0.305 2.808 0.005
  • 列の順番に、x1、x2、x3に対応している。
    • DataFrameで説明変数を入れれば、変数にラベルを入れてくれる。
  • 他は単回帰分析とほぼ同じ。
    • t値（標準正規分布に従う確率変数として、z値として書かれている）も、係数の真の値が0である、という帰無仮説の元で計算されている。

重回帰モデルの注意点

  • モデルがあるから、係数が意味を持つ。
    • 仮説を設定して、仮説を検定する、という流れを大切に。
  • 興味のある変数を絞る。
    • 他の「興味はないけど入れる変数」はコントロール変数と呼ばれるもの。
    • コントロール変数の係数は解釈しない。
    • より詳しいお話は計量経済学のテキストを読んで勉強する必要があります。

まとめ

  • 回帰分析はたいてい仮説検定がセット。
    • 検定のフローを頭に叩き込むべし。
    • 検定の基本は、真の\beta_1=0を仮定して、背理法で\beta_1 \neq 0を示す、ということ。
  • statsmodelsを用いた回帰分析の基本を行なった。
    • 最初は違和感があるかもしれないが、慣れれ特に難しいことはない。
    • コーディングでは適切に調べながら自分の行いたい操作を実装すればよい。
  • Pythonは重回帰分析についても実装が簡単。