R,Python,Stataで重回帰分析：Wooldridge (2019) で追加されたExamplesを実装してみる

Rで計量経済学を勉強できる本の一つにHeiss(2016) Using R for Introductory Econometrics.があります*1。

この本はWooldridge Introductory Econometrics 6th edition.に完全準拠しています。この本は例題として、実際のdatasetを使った分析の結果を用いて説明をしているのですが、Heiss(2016)はこの実例の全てについてRによる再現を提供しています。

一方で、その後Wooldridgeの本は執筆時点で7th editionまで出ており、内容のアップデートと同時に例題も追加されています。
この記事はHeiss(2016)でカバーされていない例題を、Heiss(2016)と同様のR、勉強がてらのPython・Stataで再現したものです。この記事では追加された例題のうち、重回帰分析の例である例題3.7 / 4.11 / 7.13を扱います。
単回帰分析の記事はこちら。

7th Ed.で追加された例題
例題3.7 / 4.11 / 7.13
補足：R/Python/Stataのコード
- R
- Python
- Stata

7th Ed.で追加された例題

追加された例題は以下の4つです。（例題以外にも章構成が変わったり練習問題が多く追加されていたりします。別記事参照）

例題2.14 p.55

第2章:単回帰分析＞2-7 2値変数を用いた回帰分析＞2-7a 反事実アウトカム・因果性と政策分析

例題3.7 p.101

第3章:重回帰分析の推定＞3-7 重回帰分析が利用できるシナリオ＞3-7e 反事実アウトカム・因果性と政策分析

例題4.11 p.151

第4章:重回帰分析の統計的推測＞4-7 因果効果と政策分析への再訪

例題7.13 p.248

第7章:質的変数に対する回帰分析＞7-6 政策分析・プログラム評価に関する補足＞7-6a プログラム評価と回帰調整

（参照: Wooldridge(2019) Introductory Econometrics. 7th edition. 邦訳は筆者）

これらは全て、米国における職業訓練プログラムの評価をしているものです。職業訓練プログラムへの参加は労働者のスキルアップを促し、収入が増加すると予想されますが、4つの例題は2種類のデータセットを用いて、この関係が本当に生じているかを確かめています。
例題2.14を除く3題は利用しているデータセットが共通なのと、例題2.14だけ単回帰分析・ほかは重回帰分析なので、その2グループにまとめて再現しています。

例題3.7 / 4.11 / 7.13

この問題はJTRAIN98というデータセットを利用しています。このデータにおける職業訓練プログラムは1997年に行われたもので、処置の割り当ては以前の労働市場におけるアウトカムによって一部決定されますが、一部自発的な参加も含まれるという特徴があり、単回帰分析のJTRAIN2データセットの事例と違って完全なランダム割り当てができていない事例です。

データの読み込み

JTRAIN98データセットはWooldridge 7thのサポートページで公開されています。

多くのデータセットがcsv, R, STATA, EViewsの4フォーマットで公開されているのですが、なぜかJTRAIN98データセットはstataデータファイルしか公開されていません。最新版だからですかね…？
というわけで、サポートページからjtrain98.dtaファイルを入手しておきます。

Rでの実装

6thまでにテキスト内に登場したデータセットはRのパッケージとして提供されているのですが、JTRAIN98は7thで新しく加わったデータセットなのでパッケージに含まれていません。なので、上で入手したStataのデータファイルを読み込むことにします。
stata dtaファイルを読み込むにはパッケージが必要です。ここではhavenパッケージのread_dta()関数を使ってみます*2。

> library(haven)
> jtrain98 <- read_dta('jtrain98.dta')
> attach(jtrain98)
 以下のオブジェクトは jtrain98 (pos = 3) からマスクされています: 

     age, black, earn96, earn98, educ, hisp, married, train, unem96, unem98 

 以下のオブジェクトは jtrain98 (pos = 6) からマスクされています: 

     age, black, earn96, earn98, educ, hisp, married, train, unem96, unem98

これでデータセットの固定まで済みました。

Pythonでの実装

pythonでStataデータファイルを使い時はpandasライブラリ内のread_stataを利用しました。
dtaファイルを読み込む関数は他にも、StatsModels内のiolib.foreignにあるgenfromdta()やStataReader()も使えそうな気配がします*3。

>>> import pandas as pd
>>> jtrain98 = pd.read_stata('jtrain98.dta')

Stataでの実装

Stataデータファイルなので1行で済みます。

use "jtrain98.dta",clear

試しに単回帰分析でTreatment Effectを測ってみる（例題3.7）

今回は職業訓練への参加が収入に与える影響を調べたいという文脈なので、少なくとも「収入」「参加の有無」の2変数が必要です。jtrain98のデータセットの中の変数trainは、1997年の職業訓練に参加していれば1、参加していなければ0をとる2値変数です。
Random Assignmentが成立していた例題2.14では、単回帰分析で処置の効果が測れました。ここでも同じように、1998年の年収（earn98）を非説明変数、trainを説明変数として単回帰分析をやってみると、職業訓練に参加した人の平均収入は参加していない人より$2,050も低いという結果になります。
例えば過去の収入が低い人ほど、収入を上げたいので職業訓練に参加する一方で、参加しない人たちは継続的に平均収入が高い、ということも十分あり得るので、係数マイナスもあり得る話です。同時に単回帰分析の係数を処置効果と解釈するのは厳しそうだということが分かります。

単回帰分析については例題2.14の実装と似たような形です。
参考: R,Python,Stataで単回帰分析：Wooldridge (2019) で追加されたExamplesを実装してみる - Sitting in an Armchair

Rでの実装

> SLR = lm(earn98 ~ train)
> summary(SLR)$coefficients
             Estimate Std. Error   t value      Pr(>|t|)
(Intercept) 10.609896  0.2794290 37.969913 6.732456e-204
train       -2.050053  0.4844142 -4.232026  2.504019e-05

Pythonでの実装

>>> import statsmodels.api as sma
>>> SLR = sma.OLS(jtrain98["earn98"], sma.add_constant(jtrain98["train"])).fit()
>>> print(SLR.summary().tables[1])
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         10.6099      0.279     37.970      0.000      10.062      11.158
train         -2.0501      0.484     -4.232      0.000      -3.001      -1.100
==============================================================================

Stataでの実装

. regress earn98 train, noheader
------------------------------------------------------------------------------
      earn98 |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       train |  -2.050053   .4844142    -4.23   0.000    -3.000507   -1.099599
       _cons |    10.6099    .279429    37.97   0.000     10.06164    11.15816
------------------------------------------------------------------------------

. display e(r2)
.01562954

重回帰分析

平均収入には処置以外にも、過去の平均年収earn96、教育年数educ、年齢age、結婚歴marriedなどが関与しそうですし、これらの属性は職業訓練に参加するかどうかの意思決定にも影響しそうです。なので、この4変数を説明変数に加えて重回帰分析をやってみます。

結果としては、職業訓練に参加した人は平均$2,410だけ平均年収が「高い」ことが判明し、単回帰分析と結果がひっくり返りました！関連しそうな変数を調整することで推定結果が大きく変わることがよく分かります。
加えた変数の係数については、符号は直感と一致しています：
- 1996年に年収が高いほど1998年にも年収が高い傾向
- 教育年数が長いほど収入が高い
- 結婚していると平均的に年収が高い

決定係数は0.016から0.405に大幅に向上し、まだ説明しきれない部分は大きいですが、4変数のコントロールはなかなか効いていると言えそうです。

Rでの実装

summary()の出力は$をつけて様々な要素にアクセスできます。例えば係数・t値などの表はcoefficients、決定係数R2はr.squaredで値を取り出します。

> MLR = lm(earn98~ train + earn96 + educ + age + married)
> summary(MLR)$coefficients # trainの係数は約2.41
              Estimate Std. Error   t value     Pr(>|t|)
(Intercept)  4.6670417 1.14528335  4.075010 4.924596e-05
train        2.4105467 0.43526249  5.538145 3.804549e-08
earn96       0.3725384 0.01862622 20.000750 2.260190e-76
educ         0.3628329 0.06404696  5.665107 1.865496e-08
age         -0.1810460 0.01887503 -9.591830 5.367737e-21
married      2.4817194 0.42626254  5.822044 7.578651e-09

> summary(SLR)$r.squared
[1] 0.01562954
> summary(MLR)$r.squared
[1] 0.404966

Pythonでの実装

pythonの回帰結果一部取り出しは回帰結果に直接.rsquaredや.paramsをつければ抽出が可能です。係数関連をまとめた表を一気に取り出したい場合は、summary()内のtablesを個別に抽出してprintすれば出力できます。

>>> MLR = sma.OLS(jtrain98["earn98"], sma.add_constant(jtrain98[["train","earn96","educ","age","married"]])).fit()
>>> print(MLR.summary().tables[1])
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          4.6670      1.145      4.075      0.000       2.420       6.914
train          2.4105      0.435      5.538      0.000       1.557       3.265
earn96         0.3725      0.019     20.001      0.000       0.336       0.409
educ           0.3628      0.064      5.665      0.000       0.237       0.488
age           -0.1810      0.019     -9.592      0.000      -0.218      -0.144
married        2.4817      0.426      5.822      0.000       1.645       3.318
==============================================================================

>>> print(SLR.rsquared)
0.01562954996359811
>>> print(MLR.rsquared)
0.4049660431630191

Stataでの実装

回帰結果については、係数関係だけ表示したいならregressコマンドにnoheaderをつければよいだけです*4。
また、回帰結果の個別数値に関してはe()に保存されています*5。再度regressを回すまでの一時的な保存です。

. regress earn98 train earn96 educ age married, noheader
------------------------------------------------------------------------------
      earn98 |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       train |   2.410547   .4352625     5.54   0.000     1.556528    3.264565
      earn96 |   .3725384   .0186262    20.00   0.000     .3359923    .4090845
        educ |   .3628329    .064047     5.67   0.000     .2371678    .4884979
         age |   -.181046    .018875    -9.59   0.000    -.2180803   -.1440118
     married |   2.481719   .4262625     5.82   0.000      1.64536    3.318079
       _cons |   4.667042   1.145283     4.08   0.000     2.419908    6.914176
------------------------------------------------------------------------------

. display e(r2)
.40496603

統計学的有意性を確認する（例題4.11）

上で確認した結果の表からt値等も確認できます。みたところ、単回帰分析の負の効果、重回帰分析の正の効果ともに有意な差があると言えそうです。t値は単回帰分析の係数は-4.27、重回帰分析の係数は5.47と、双方ともに有意なのです。この場合は直感と合っていそうな重回帰分析の結果を信じたいところです*6。

また、ここでのコントロール変数は処置効果（trainの係数）を不偏推定するために加えているので、コントロール変数の係数やその統計的有意性、それらに影響する多重共線性等は気にする必要がありません！とのことです。誤解しがちな部分だと思うので改めて確認できてよかった…。

treatment effect推定における「回帰調整」（例題7.13）

平均処置効果（average treatment effect: ATE）の推定の際に、重回帰分析で他の条件・変数をコントロールすることを、プログラム評価の文脈では「回帰調整」（Regression Adjustment: RA）というらしいです。上の例題二つで確認した通り、回帰調整した方が妥当なATE推定になりそうです。
一方で、コントロール変数を含めた重回帰分析でのtrainの係数が平均処置効果ATEを表している、というときの前提として、職業訓練の効果は他の条件によらず一定という暗黙の仮定があります。ただ、例えば収入が低かった人の方が収入を上げる効果が大きいとか、教育を受けてきた人の方が職業訓練を効果的に吸収できるとか、条件によってtreatment effectが違ってくることも十分考えられます。
その制約を緩める方法として、各変数と処置ダミーの交差項を入れるという方法があります。詳しくは他記事で。交差項を入れることで制約を外しているので、交差項なしの重回帰分析をRestricted Regression Adjustment (RRA)、交差項を入れたものをUnrestricted Regression Adjustment (URA)と呼ぶそうです。
この例題では交差項の導入、及び「処置効果が属性によって異なる」という仮説に対する検定を行なっています。交差項を入れると、ATEはRRAの$2,410に対してURAでは$3,106となり、t値もRRAとURAでほぼ同レベルに統計的有意性があるという結果になりました。一方で、職業訓練の効果は個人の属性によって差がない（＝交差項の係数は全てゼロ）という仮説のF統計量はp=0.113と有意ではありませんでした。

Rでの実装

F検定にはcarパッケージのlinearHypothesis関数を使いました。制約の表記法にはいくつかありますが、制約を表す文字列をベクトルにして放り込んでやれば検定してくれるので、かなり楽です。

> ## 変数作成
> adj.earn96 = (earn96-mean(earn96))*train
> adj.educ = (educ-mean(educ))*train
> adj.age = (age-mean(age))*train
> adj.married = (married-mean(married))*train

> ## 重回帰分析
> ADJedMLR = lm(earn98~ train + earn96 + educ + age + married
+               + adj.earn96 + adj.educ + adj.age + adj.married)
> summary(ADJedMLR)$coefficients
               Estimate Std. Error    t value     Pr(>|t|)
(Intercept)  5.07757354 1.38948699  3.6542793 2.698517e-04
train        3.10620223 0.53211027  5.8375160 6.932947e-09
earn96       0.35336167 0.02004416 17.6291579 1.404337e-61
educ         0.37783461 0.07793423  4.8481215 1.421998e-06
age         -0.19606086 0.02278920 -8.6032359 2.581569e-17
married      2.75907670 0.54508712  5.0617170 4.853862e-07
adj.earn96   0.13275474 0.05415212  2.4515151 1.437754e-02
adj.educ    -0.03487017 0.13657627 -0.2553164 7.985256e-01
adj.age      0.05806490 0.04079468  1.4233449 1.549147e-01
adj.married -0.99320920 0.88273635 -1.1251482 2.607673e-01

> ## F検定（H0:全ての交差項の係数=0）
> library(car)
> linearHypothesis(ADJedMLR, c("adj.earn96 = 0","adj.educ = 0","adj.age = 0","adj.married = 0"))
Linear hypothesis test

Hypothesis:
adj.earn96 = 0
adj.educ = 0
adj.age = 0
adj.married = 0

Model 1: restricted model
Model 2: earn98 ~ train + earn96 + educ + age + married + adj.earn96 + 
    adj.educ + adj.age + adj.married

  Res.Df   RSS Df Sum of Sq      F Pr(>F)
1   1124 40143                           
2   1120 39876  4    266.79 1.8733 0.1128

Pythonでの実装

F検定はStatsModels内RegressionResults.compare_f_test関数を使います*7
。制約付きモデルと制約なしモデルの両方を推定して比較する形式で実行する形なので、制約を指定するRの関数とはスタイルが違います。この関数は2つのモデルがNestedである必要があります。(F値、p値、自由度)の3要素のタプルを返します。2番目の要素が0.1128なので非有意ということがわかります。

>>> jtrain98['adj.earn96'] = ( jtrain98['earn96'] - jtrain98['earn96'].mean() ) * jtrain98['train']
>>> jtrain98['adj.educ'] = ( jtrain98['educ'] - jtrain98['educ'].mean() ) * jtrain98['train']
>>> jtrain98['adj.age'] = ( jtrain98['age'] - jtrain98['age'].mean() ) * jtrain98['train']
>>> jtrain98['adj.married'] = ( jtrain98['married'] - jtrain98['married'].mean() ) * jtrain98['train']

>>> ADJedMLR = sma.OLS(jtrain98["earn98"], sma.add_constant(jtrain98[["train","earn96","educ","age","married","adj.earn96","adj.educ","adj.age","adj.married"]])).fit()
>>> print(ADJedMLR.summary().tables[1])
===============================================================================
                  coef    std err          t      P>|t|      [0.025      0.975]
-------------------------------------------------------------------------------
const           5.0776      1.389      3.654      0.000       2.351       7.804
train           3.1062      0.532      5.838      0.000       2.062       4.150
earn96          0.3534      0.020     17.629      0.000       0.314       0.393
educ            0.3778      0.078      4.848      0.000       0.225       0.531
age            -0.1961      0.023     -8.603      0.000      -0.241      -0.151
married         2.7591      0.545      5.062      0.000       1.690       3.829
adj.earn96      0.1328      0.054      2.452      0.014       0.027       0.239
adj.educ       -0.0349      0.137     -0.255      0.799      -0.303       0.233
adj.age         0.0581      0.041      1.423      0.155      -0.022       0.138
adj.married    -0.9932      0.883     -1.125      0.261      -2.725       0.739
===============================================================================

>>> ADJedMLR.compare_f_test(MLR)
(1.873336195416279, 0.11281450612498053, 4.0)

Stataでの実装

仮説検定にはtestコマンドを使うのですが、regressコマンドの後にtestコマンドで帰無仮説を列挙すればいいだけで、とてもシンプルです。書き方のバリエーションは様々あります*8。

. egen avg_earn96 = mean(earn96)
. gen adj_earn96 = (earn96 - avg_earn96) * train
. egen avg_educ = mean(educ)
. gen adj_educ = (educ - avg_educ) * train
. egen avg_age = mean(age)
. gen adj_age = (age - avg_age) * train
. egen avg_married = mean(married)
. gen adj_married = (married - avg_married) * train

. regress earn98 train earn96 educ age married adj_earn96 adj_educ adj_age adj_married

      Source |       SS           df       MS      Number of obs   =     1,130
-------------+----------------------------------   F(9, 1120)      =     86.09
       Model |   27586.969         9  3065.21878   Prob > F        =    0.0000
    Residual |  39875.9225     1,120  35.6035022   R-squared       =    0.4089
-------------+----------------------------------   Adj R-squared   =    0.4042
       Total |  67462.8915     1,129   59.754554   Root MSE        =    5.9669

------------------------------------------------------------------------------
      earn98 |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       train |   3.106202   .5321103     5.84   0.000     2.062157    4.150248
      earn96 |   .3533617   .0200442    17.63   0.000     .3140333      .39269
        educ |   .3778346   .0779342     4.85   0.000     .2249211    .5307481
         age |  -.1960609   .0227892    -8.60   0.000    -.2407752   -.1513465
     married |   2.759077   .5450871     5.06   0.000      1.68957    3.828584
  adj_earn96 |   .1327547   .0541521     2.45   0.014     .0265037    .2390058
    adj_educ |  -.0348702   .1365763    -0.26   0.799    -.3028443     .233104
     adj_age |   .0580649   .0407947     1.42   0.155    -.0219777    .1381075
 adj_married |  -.9932092   .8827363    -1.13   0.261    -2.725212     .738794
       _cons |   5.077574   1.389487     3.65   0.000     2.351283    7.803864
------------------------------------------------------------------------------

. test adj_earn96 adj_educ adj_age adj_married

 ( 1)  adj_earn96 = 0
 ( 2)  adj_educ = 0
 ( 3)  adj_age = 0
 ( 4)  adj_married = 0

       F(  4,  1120) =    1.87
            Prob > F =    0.1128

補足：R/Python/Stataのコード

R

library(haven)
library(car)
jtrain98 <- read_dta('jtrain98.dta')

attach(jtrain98)

# 単回帰分析（例題3.7/4.11）
SLR = lm(earn98 ~ train)
summary(SLR)$coefficients # trainの係数は約-2.05

# 重回帰分析（例題3.7/4.11）
MLR = lm(earn98~ train + earn96 + educ + age + married)
summary(MLR)$coefficients # trainの係数は約2.41

summary(SLR)$r.squared # 単回帰のR^2は約0.016
summary(MLR)$r.squared # 重回帰のR^2は約0.405

# 重回帰分析と回帰調整（例題7.13）
## 変数作成
adj.earn96 = (earn96-mean(earn96))*train
adj.educ = (educ-mean(educ))*train
adj.age = (age-mean(age))*train
adj.married = (married-mean(married))*train

## 重回帰分析
ADJedMLR = lm(earn98~ train + earn96 + educ + age + married
              + adj.earn96 + adj.educ + adj.age + adj.married)
summary(ADJedMLR)$coefficients # trainの係数は約3.11

## F検定（H0:全ての交差項の係数=0）
linearHypothesis(ADJedMLR, 
                 c("adj.earn96 = 0","adj.educ = 0",
                   "adj.age = 0","adj.married = 0"))

detach(jtrain98)

Python

# ライブラリの読み込み
import statsmodels.api as sma
import pandas as pd

# データセットの読み込み
jtrain98 = pd.read_stata('jtrain98.dta')

# 単回帰分析 （例題3.7/4.11）
SLR = sma.OLS(jtrain98["earn98"], sma.add_constant(jtrain98["train"])).fit()
print(SLR.summary().tables[1]) # trainの係数は約-2.05

# 重回帰分析 （例題3.7/4.11）
MLR = sma.OLS(jtrain98["earn98"], sma.add_constant(jtrain98[["train","earn96","educ","age","married"]])).fit()
print(MLR.summary().tables[1]) # trainの係数は約2.41

print(SLR.rsquared) # 単回帰のR^2は約0.016
print(MLR.rsquared) # 重回帰のR^2は約0.405

# 回帰調整付き重回帰分析
## 変数作成
jtrain98['adj.earn96'] = ( jtrain98['earn96'] - jtrain98['earn96'].mean() ) * jtrain98['train']
jtrain98['adj.educ'] = ( jtrain98['educ'] - jtrain98['educ'].mean() ) * jtrain98['train']
jtrain98['adj.age'] = ( jtrain98['age'] - jtrain98['age'].mean() ) * jtrain98['train']
jtrain98['adj.married'] = ( jtrain98['married'] - jtrain98['married'].mean() ) * jtrain98['train']
## 回帰分析
ADJedMLR = sma.OLS(jtrain98["earn98"], sma.add_constant(jtrain98[["train","earn96","educ","age","married","adj.earn96","adj.educ","adj.age","adj.married"]])).fit()
print(ADJedMLR.summary().tables[1]) # trainの係数は約3.11

ADJedMLR.compare_f_test(MLR) # p-value > 0.1

Stata

/* データセットの読み込み */
use "jtrain98.dta",clear

/* 単回帰分析 */
regress earn98 train, noheader
display e(r2)

/* 重回帰分析 */
regress earn98 train earn96 educ age married, noheader
display e(r2)

/* 回帰調整付き重回帰分析 */
/** 変数作成 **/
egen avg_earn96 = mean(earn96)
gen adj_earn96 = (earn96 - avg_earn96) * train
egen avg_educ = mean(educ)
gen adj_educ = (educ - avg_educ) * train
egen avg_age = mean(age)
gen adj_age = (age - avg_age) * train
egen avg_married = mean(married)
gen adj_married = (married - avg_married) * train

/** 回帰分析 **/
regress earn98 train earn96 educ age married adj_earn96 adj_educ adj_age adj_married

test adj_earn96 adj_educ adj_age adj_married