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で再現したものです。この記事では追加された例題のうち、単回帰分析の例である例題2.14を扱います。
重回帰分析の記事はこちら。
scientia-socialis-f-discolor.hatenablog.com
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グループにまとめて再現しています。
例題2.14
この問題はJTRAIN2というデータセットを利用しています。1970年代に行われた職業訓練プログラムに関するデータセットなのですが、この職業訓練プログラムは処置をランダムに割り当てているという特徴があり、要は珍しくRCTができた事例です。
（元論文：R.J. Lalonde (1986), “Evaluating the Econometric Evaluations of Training Programs with Experimental Data,” American Economic Review 76, 604-620. ）
データの読み込み
Rでの実装
6thまでにテキスト内に登場したデータセットはRのパッケージとして提供されています。JTRAIN2はパッケージに収録されているので、そこから読み込むことにします*2。
> # パッケージのインストール（初回のみ） > # install.packages("wooldridge") > library(wooldridge) # パッケージの読み込み Warning message: パッケージ ‘wooldridge’ はバージョン 3.6.1 の R の下で造られました > data('jtrain2')
これでjtrain2という名前でデータセットが読み込めました。
Stataでの実装
Wooldridgeの6thまでに登場したデータセットは、Stataデータファイルがweb上でダウンロードできます。Stataデータセットを開くuse関数に直接URLを打ち込んで実行すればOKです。
. use "http://fmwww.bc.edu/ec-p/data/wooldridge/jtrain2.dta"
Pythonでの実装
Rと違ってデータセットがパッケージで提供されているわけではなく、Stataと違ってオンラインに転がっているデータセットを引っ張ってこられもしなかった*3ので、素直にデータセットをWooldridgeのサポートページからダウンロードして使います。
データセットはStata, EViews, R, Excel 97-2003の4種類が提供されていますが、ここは素直にExcelを使ってみます。
xlsファイルの読み込みにはpandasモジュールのpandas.read_excel()関数を使うのですが、この関数は内部でxlrdというライブラリを使っているので、事前にインストールが必要です。
表頭1行目からデータが格納されているので、pandas.read_excel関数内ではheader = Noneを指定して表頭からデータとして読み込みます。Excelデータには列名が含まれていないので外付けで指定します（面倒ですがｗ）。
>>> import pandas as pd >>> jtrain2 = pd.read_excel('jtrain2.xls', header = None) >>> jtrain2.columns = ['train', 'age', 'educ', 'black', 'hisp', 'married', 'nodegree', 'mosinex', 're74', 're75', 're78', 'unem74', 'unem75', 'unem78', 'lre74', 'lre75', 'lre78', 'agesq', 'mostrn']
参考：pandasでExcelファイル（xlsx, xls）の読み込み（read_excel） | note.nkmk.me
変数の内容の確認
今回は職業訓練への参加が収入に与える影響を調べたいという文脈なので、「収入」「参加の有無」の2変数が必要です。jtrain2のデータセットの中の変数trainは、職業訓練に参加していれば1、参加していなければ0をとる2値変数です。職業訓練に参加した人の人数は、「trainの中の1の数を数える」「trainの総和を計算する」などの方法でチェックできます。ここではどの言語でも総和計算でなく0、1それぞれをカウントしています。
Rでの実装
変数の中身のカウントには、総数ならnrow()で列数をカウント、数値ごとにカウントならtable()を使います。
> nrow(jtrain2) [1] 445 > ## サンプルサイズは n = 445 > table(jtrain2$train) 0 1 260 185 > ## 185人が職業訓練に参加していて、260人は参加していない
Pythonでの実装
カウントにはcollectionsライブラリのCounterメソッドを使います。
>>> len(jtrain2) ## サンプルサイズはn=445 445 >>> import collections >>> collections.Counter(jtrain2["train"]) ## 185人が職業訓練プログラムに参加、260人が不参加 Counter({0: 260, 1: 185})
Stataでの実装
obs数はデータをインポートした時点でシステム変数_Nに入っています。数値を表示したいときはdisplayしてあげる必要があります。
カウントにはRと同じくtableコマンドを使います。
. display _N /** サンプルサイズは n = 445 **/ 445 . table train /** 185人が職業訓練プログラムに参加、260人が不参加 **/ ---------------------- =1 if | assigned | to job | training | Freq. ----------+----------- 0 | 260 1 | 185 ----------------------
単回帰分析
ここでは1978年の実質収入（re78）に対する1970年代の職業訓練の影響を調べたいので、re78を目的変数、trainを説明変数にして回帰分析を行います。
回帰分析の結果から、職業訓練プログラムに参加した人は参加していない人に比べて平均で$1,790収入が多いことがわかります*4。実際、プログラムに参加していない人の平均収入と参加した人の平均収入を比較すると、差は同じ値になります。
これは参加していない人と比べて39.3%も大きいので、費用便益比較はいったん置いておくとしても、無視できない効果があるといえるでしょう。
Rでの実装
回帰分析はlm関数を使います。
attach()関数は、detachするまでデータセットを固定でき、データフレーム名を毎回入力する必要がなくなります。継続して同じデータセットを扱う場合に便利です。挙動としてはStataに近くなるでしょうか。
> attach(jtrain2) > result = lm(re78 ~ train) # 単回帰分析 > summary(result)$coefficients # 回帰分析の結果のうち係数関係を表示 Estimate Std. Error t value Pr(>|t|) (Intercept) 4.554802 0.4080460 11.162474 1.154113e-25 train 1.794343 0.6328536 2.835321 4.787524e-03
プログラム参加（train = 1）、不参加（train = 0）ごとに平均収入を計算するなど、「○○ごとに△△の□□を計算する」といった処理にはapply関数の仲間が便利です。tapply(re78,train,mean)で、「trainごとにre78のmeanを計算する」という処理になります。
無処置群の平均収入は回帰結果の定数項、平均収入の群間の差は上の回帰結果内のtrainの係数とそれぞれ一致していることがわかります。
参考：applyファミリー | R で同じ処理を”並列的”に実行する関数
> tapply(re78,train,mean) 0 1 4.554802 6.349145 > tapply(re78,train,mean)[2]-tapply(re78,train,mean)[1] 1 1.794343
増分を不参加者の平均収入で割れば、割合で表した増加分が出ます。
> coef <- summary(result)$coefficients[,1] > coef[2]/coef[1]*100 train 39.39453
attach()関数を使ったので、最後に固定を解除して終わります。
> # 最後にデータセットの固定を解除して終了 > detach(jtrain2)
Pythonでの実装
- 回帰分析にはStatsModelsモジュールを利用しています*5。
- statsmodels.apiのOLS()関数でモデルを指定し、.fit()を後につけることで推定を行うという構造の模様。OLS()関数の中身は被説明変数,説明変数という感じの表記になるのですが、定数項を含めて推定したいときは説明変数にadd_constant()関数をひっかけてやる必要があります。
- summary()関数を使えば回帰結果に関する様々な指標が出てきます。ここでは係数結果だけ参照したかったので、tablesのうち係数関係の表のみ抜粋して示しています。
>>> import statsmodels.api as sma >>> result = sma.OLS(jtrain2["re78"], sma.add_constant(jtrain2["train"])).fit() >>> print(result.summary().tables[1]) ## 回帰分析結果のうち係数関係のみ表示 ============================================================================== coef std err t P>|t| [0.025 0.975] ------------------------------------------------------------------------------ const 4.5548 0.408 11.162 0.000 3.753 5.357 train 1.7943 0.633 2.835 0.005 0.551 3.038 ==============================================================================
群ごとの平均収入計算にはgroupby()メソッドを使います。この計算結果はpandasのdataframeで出力されるので、ilocメソッドでdataframeの一部を抜粋したものをfloat型に変換してあげれば、その数値自体を計算に使えます。
>>> means = jtrain2[["re78","train"]].groupby(["train"]).mean() >>> float(means.iloc[1]) - float(means.iloc[0]) 1.7943430604261987
回帰結果に.paramsをくっつけることで係数のみ参照できます。係数はリストで返ります。係数の比から、比率で表した平均収入増分効果が計算できます。
>>> coef = result.params >>> coef[1] / coef[0] * 100 39.39453232062628
Stataでの実装
- 回帰分析はregressコマンドを使います。regと省略しても認識してくれます。オプションにnoheaderをつけると、係数関連のResultだけ表示してくれます。
- 増分（trainの係数）を不参加者の平均収入（定数項）で割れば増加割合が出ます。stataではシステム変数_b、_constを使えば係数・定数項の数値にアクセスできて、これを使えば数値を直接ベタ打ちする必要はありません。
- 「○○ごとに△△の□□を計算する」といった処理にはtabstat コマンドを使います。オプションにnototalをつけると全体の平均値を省略してくれます。
. regress re78 train, noheader /** 回帰結果のうち係数関係のみ表示 **/ ------------------------------------------------------------------------------ re78 | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- train | 1.794343 .6328536 2.84 0.005 .5505748 3.038111 _cons | 4.554802 .408046 11.16 0.000 3.752856 5.356749 ------------------------------------------------------------------------------ . tabstat re78, by(train) nototal Summary for variables: re78 by categories of: train (=1 if assigned to job training) train | mean ---------+---------- 0 | 4.554802 1 | 6.349145 -------------------- . display _b[train] / _b[_cons] * 100 /** training program による収入増加割合は約39.4%！ **/ 39.394533
まとめ・感想
以上が3言語の単回帰実装です。
この例題はself-selectionのない政策介入における回帰分析による政策効果推定の例でした。重回帰分析のほうでは完全なrandom assignmentになっていない事例に重回帰分析を適用していく例になっています。
補足：各言語のコード
R
# データセットの読み込み library(wooldridge) data('jtrain2') # 変数の内容確認 nrow(jtrain2) ## サンプルサイズはn=445 table(jtrain2$train) ## 185人が職業訓練プログラムに参加、260人が不参加 # 単回帰分析 attach(jtrain2) ## データセットの固定 result = lm(re78 ~ train) summary(result)$coefficients ## 回帰分析結果のうち係数関係のみ表示 tapply(re78,train,mean) tapply(re78,train,mean)[2]-tapply(re78,train,mean)[1] ## 各群の平均収入の差は回帰係数と一致 coef <- summary(result)$coefficients[,1] coef[2]/coef[1]*100 ## training program による収入増加割合は約39.4%！ # 最後にデータの固定を解除して終了 detach(jtrain2)
Python
# ライブラリの読み込み import statsmodels.api as sma import pandas as pd import collections # データセットの読み込み jtrain2 = pd.read_excel('jtrain2.xls', header = None) jtrain2.columns = ['train', 'age', 'educ', 'black', 'hisp', 'married', 'nodegree', 'mosinex', 're74', 're75', 're78', 'unem74', 'unem75', 'unem78', 'lre74', 'lre75', 'lre78', 'agesq', 'mostrn'] # 変数の内容確認 len(jtrain2) ## サンプルサイズはn=445 collections.Counter(jtrain2["train"]) ## 185人が職業訓練プログラムに参加、260人が不参加 # 単回帰分析 result = sma.OLS(jtrain2["re78"], sma.add_constant(jtrain2["train"])).fit() print(result.summary().tables[1]) ## 回帰分析結果のうち係数関係のみ表示 means = jtrain2[["re78","train"]].groupby(["train"]).mean() float(means.iloc[1]) - float(means.iloc[0]) ## 各群の平均収入の差は回帰係数と一致 coef = result.params coef[1] / coef[0] * 100 ## training program による収入増加割合は約39.4%！
Stata
/* データセットの読み込み */ use "http://fmwww.bc.edu/ec-p/data/wooldridge/jtrain2.dta" /* 変数の内容確認 */ display _N /** サンプルサイズは n = 445 **/ table train /** 185人が職業訓練プログラムに参加、260人が不参加 **/ /* 単回帰分析 */ regress re78 train, noheader /** 回帰結果のうち係数関係のみ表示 **/ tabstat re78, by(train) nototal save display _b[train] / _b[_cons] * 100 /** training program による収入増加割合は約39.4%！ **/
*1:2020年5月末現在、第2版が公開されオンラインで読めるようになっています。2nd editionは2019年のWooldridge 7th editionに準拠して書いているとありますが、全問題には対応してないみたいです。
*2:他の方法として、6thまでに登場するデータセットはオンラインで公開されているので、リンク先から直接インポートするという手もあります。その場合のコードは割愛しますが、結局パッケージは必要ですし、インターネット環境が必要で処理が重いです。
*3:Web上のコンテンツを入手できるrequestモジュールのget関数（参照）や、dtaファイルを読み込むpandasモジュールのread_dtaもあるのですが、どうにも読み込めなかったので諦めました。
*4:これがテキストで述べられている2値変数の係数の解釈です。
*5:機械学習向けモジュールのscikit-learnを使ってもできるようですが、StatsModelsのほうが統計表出力のスタイルなどの取り回し方が馴染む気がしました。統計学・計量経済学畑の人はこちらという意識なのかもしれません
