その回帰分析、本当に因果関係ありますか? - 操作変数法の話 -

今回は統計学(計量経済学)の手法の1つ操作変数法について話したいと思います。その前に、具体例として、賃金関数と呼ばれる理論についてお話します。Jacob Mincer が70年代に発表した研究で、人の収入は、教育年数(≒学歴)である程度決まる、というものです。理屈としては、

• 教育を受け、より高度な知識を得ることで、より付加価値のある仕事ができる
• よって、対価として得られる給料も増える

というものです。
これだけでは、単なる空想にすぎません。これが現実世界の人の収入を決定しているのか、ということを確かめるため、統計的に検証する必要があります。そこで回帰分析です。
データは、国内のものがあればよかったのですが、個人レベルのデータを入手するのは難しいので、Blackburn と Neumark の 1992年の研究で使用されたデータ(の一部)を使います。こちらはインターネット上でダウンロードできます。
SAS で回帰分析を行うには、regプロシジャを使います。以前紹介した、 SAS University Edition を使ってやってみます。

/* .xls形式のデータセットを読み込む*/
proc import datafile='/folders/myfolders/grilic.xls'
   out     =work.grilic 
   dbms=xls replace;
 getnames=YES;
run;
/* 回帰分析*/
proc reg data=work.grilic;
 wage1:model lw = S;
quit;

これで結果が得られます。

LW は時給の対数で、Sは学校で教育を受けた年数です。これは、時給の対数としたほうが、モデルの当てはまりがよくなることが経験的に知られているからです。 また、従属変数だけを対数とすることで、係数が,「説明変数が1増えるごとに従属変数が何パーセント増加するか」を表すことになります。
さて、ここで次のような反論が現れるかもしれません。
「学校教育より現場で学んだ経験のほうが仕事をする上で重要だ」と。
そのとおりです。だから「教育年数イコール学歴」ではなく「≒学歴」と書きました。いわゆる OJT によって、人は仕事をしながらも学んでいくのです。そこで、Mincer は次のようなモデルを考えました。

では、このモデルで推定をやり直しましょう EXPRという変数に、その人がこれまで何年サラリーマンとして働いたかが記録されています。

proc reg data=work.grilic;
 wage2: model lw = S EXPR;
quit;

あらたに追加した経験年数の変数も有意でした。決定係数も大きくなっています。
ところで、給料を決定するのは果たして学校教育と職場経験の2つだけでしょうか? 例えば、特に勉強している様子のないのにやたらと成績のいい人を一度は見たことがあるでしょう。そして多分、そういう人は仕事でも要領よく立ち回れたりします。つまり、この回帰分析には、人の生まれつきの才能とか、能力を表す説明変数が欠けているのです。今回のデータセットには、「IQ」の変数も存在します。そこで、これも説明変数に加えましょう

proc reg data=grilic;
 wage3: model LW = S EXPR IQ;
quit;

結果は本当に正しいのか?

IQ を加えたことにより、係数がおよそ0.094となりました。もう1年学校で教育を受けたとしたら 9.4% だけ収入 (時給) が増えることになります。同様に、職業経験が1年増えると、収入が 4.6% 増加することになります。ただし、サンプルが10代~20代という若い世代に限定されているので注意しましょう(さらに言うと、このデータは70年代のアメリカのデータなので、この数字は日本でもピッタリ当てはまるとは限りません)。
さてこの結果は本当に正しいのでしょうか. ここまでで、説明変数を追加することで係数が少しづつ変わってきましたが、これで全てでしょうか?3つ目の説明変数として、「生まれつきの能力」が必要だと書きましたが、果たしてIQは能力の数値化として妥当でしょうか? IQ はもともと知的障碍のテストという医療目的で作られたもので、また、主に子どもを対象としたテストです。そのため、必ずしも能力の指標とはなりません。ただ、能力とある程度相関がありそうなので、 便宜的に説明変数として使われたに過ぎません(代理変数)。人間の能力を数値化するのは難しいのです。よって、IQ と実際の能力の間には、誤差があると予想できます。よって、この係数も正確な推定値ではなく、なんらかの偏りがあります。これを「観測誤差バイアス」と呼びます。
今回は IQ という変数がありましたが、そもそも能力の代理変数の情報を得られない場合があります。その場合、変数そのものが完全に欠落しているため、「欠落変数バイアス」と呼びます。経済分析では、加えて「同時決定バイアス」と呼ばれるバイアスを考慮する必要があり、これらいずれも「操作変数」を用いることで対処できますが、今回は賃金関数を例に、観測誤差バイアスの場合に限定して解説します。

各変数は、矢印の方向に影響を及ぼしていると見てください。そのうち、能力の代理変数 IQは、能力を完全には説明できておらず、誤差部分があります。この誤差部分は測定できないので、誤差項の一部として扱われてしまいます。すると、誤差項とIQの相関が発生することになります。このような変数を「内生変数」と呼びます。 内生変数がモデルに含まれている時、最小二乗法による推定は意味をなしません。そこで、内生変数に影響し、かつ誤差項には影響しない、「操作変数」Zを新たに用意する必要があります。いわば、このモデルには、本来の回帰式に加え、内生変数に対して回帰する、2つの回帰式が存在するのです。
「欠落変数バイアス」「同時決定バイアス」の流れもも基本的にこのような図で表現できます。
使用できる変数の中で、このZになりうるのは、「母親の教育年数」です。親の教育年数が長い=高学歴ならば、遺伝的に子どもも能力がある、と仮定すれば、IQと関係があるはずです。さらに、母親の教育は、おそらく子どもの賃金に影響しないでしょう。よって、操作変数として使えそうです。

2段階最小二乗法

2段階最小二乗法は、名前の通り、2段階の最小二乗法を行います。1段階目が、内生変数に対する回帰です。ここでの説明変数は、操作変数と、賃金関数の他の説明変数です。 この推定から、IQの予測値を計算し、当初の賃金関数のIQに置き換え,回帰分析を行います。これが2段階目の最小二乗法です。よって、「2段階最小二乗法」と呼ばれます。
SAS では、SYSLINプロシジャでこの手法が使えますが、残念ながら SAS/ETS に属するプロシジャのため、 SAS University Edition では使用できません (入力補完の候補には出てきますが、実行するとエラー扱いになります)。 ですが、2段階最小二乗法の原理を知っていれば、regプロシジャだけで結果を出せます。

/* 1段階目*/
proc reg data=work.grilic;
 first_stage: model IQ = S EXPR MED;
 output out=work.grilic2 predicted=IQhat residual=resid;
run;
/* 2段階目 */
proc reg data=work.grilic2;
 second_stage: model LW = S EXPR IQhat / hcc acovmethod=1;
run;

1度目のregプロシジャで内の outputプロシジャで、predicted=オプションを使用して S の予測値を grilic2 というデータセットに出力しました。これにはgrilicの元の変数も含まれています。residual=は同様に残差を出力するオプションです。何に使うかはこの後説明します。 この grilic2を利用して第2段階の推定を行います。
まず1段階目の回帰分析では、使用した操作変数が、適切かどうかを確認します。操作変数は、内生変数の説明変数とならなければならないので、相関のある変数でなければなりません。相関係数の小さい操作変数を「弱相関操作変数」といい、使用することでかえって推定の精度を悪化させます。操作変数と内生変数の相関が弱すぎる (あるいはまったく関係ない変数を操作変数として使ってしまった) ことが原因なので、1段階目の回帰分析の結果から、適切な操作変数かどうかを判断できます。1段階目の回帰分析で使用する変数は決まっているので、個別のt値ではなくモデル全体の有意性を見るF値を確認します (有意でない変数が1つ2つあってもF値が大きいなら問題なし)。今回はF値は十分に大きく、対応するp値が非常に小さいので、操作変数として適切と判断できます。
2段階目では、目当ての回帰式である賃金関数を回帰分析します。IQ は、予測値 IQhatに置き換えています。また、hcc、acovmethod=1という標準誤差の計算方法を変更するオプションが必要です。なぜ必要なのかというと話が長くなるので、大雑把に「推定された値を使用してさらに推定しているので、誤差が大きくなるから」と考えてください。実際、結果を見ると、S と IQ の標準誤差が上昇しています。通常の回帰分析(最小二乗法)と2段階最小二乗法の結果を比較して、操作変数を適用した変数の係数だけが変化していること、また決定係数がゼロやマイナスなど不自然な値になっていない、ということが、分析がうまくいっている目安になるのですが、今回は S の値も大きく変化しました。これは 「S も内生変数の可能性がある」事を意味します。「才能は仕事だけでなく学校の成績にも影響する」と考えれば、才能のある人間ほど進学する、という相関が生まれます。よって、学校教育Sによる収入アップ効果が、生まれつきの才能と混同されている、というシナリオが考えられます。
よって、S からもバイアスを排除する必要があります。2段階最小二乗法は、内生変数が複数ある場合にも対応できます。
その前に1点、重要なのが内生性の検定です。内生性の検定として、Hausman-Wu 検定が有名で、model プロシジャで実行できますが、これも University Edition では使用できないので、別の方法を用います。これは Wooldridge (2010) の131頁に詳しく書かれている手法です。1段階目の回帰分析の残差を、2段階目の回帰分析に説明変数として加えます。残差の係数が有意なら、IQは内生変数の可能性があります。厳密には「外生変数ではない」という帰無仮説を棄却することになるので、これは「内生性の検定」ではありません。これは Hausman-Wu 検定も同じです。そのため、「この変数は内生変数なので操作変数が必須である」という判断には使えません。このことから「内生性の検定は存在しない」という格言のようなものもあります。

/* 内生性テスト */
proc reg data=work.grilic2;
 endo_test: model LW = S EXPR IQ resid;
run;

複数の内生変数であっても、基本的には各内生変数に対して1段階目の回帰分析を行えばいいです。

/* 年ダミー作成 */
data work.grilic;
 set work.grilic;
 
 y_67 = (year=67);
 y_68 = (year=68);
 y_69 = (year=69);
 y_70 = (year=70);
 y_71 = (year=71);
 y_72 = (year=72);
 y_73 = (year=73);
run;
/* 1段階目 */
proc reg data=work.grilic;
 first_stage: model S IQ = EXPR TENURE RNS SMSA y_: MED KWW AGE;
 output out=work.grilic2 predicted=Shat IQhat residual=S_r IQ_r;
quit;
/* 2段階目*/
proc reg data=work.grilic2;
 TSLS: model LW = Shat EXPR TENURE IQhat RNS SMSA y_: / hcc acovmethod=1;
quit;

ここまでは、説明のために必要最低限の説明変数だけを使いましたが,最後なので、ついでに他の説明変数も追加して推定します。追加するのは,

• 現職の勤続年数 TENURE
• アメリカ南部在住を表すダミー RNS
• 都市圏に在住しているかのダミー SMSA
• 年度ダミー y_67, ... y_73

です。最初のデータステップは、年度を表す変数 YEAR を元に年度ダミー変数を作成しています。クロスセクションデータの場合、よくこのような時間のダミー変数を使います。金融市場のような長期間の時系列データの場合はより丁寧な分析が必要ですが、今回は短い期間の話なので、ダミー変数で代用できます。
さらに、新たに2つの操作変数を追加します。

• 職業知識テストのスコア KWW
• 年齢 AGE

です。
model ステートメントでは, 左辺を複数の変数を書くと, それぞれに対して回帰分析を実行してくれます。これで, 内生変数 S、IQの両方を回帰します。
さて結果ですが、1段階目では見慣れないメッセージが表示されました。

これは、多重共線性と呼ばれる現象で、説明変数どうしの相関係数が大きいと、計算ができなくなる場合があります(特異行列の問題です)。今回は年度ダミー変数どうしで多重共線性が発生したようです。そのままでは計算できないので、SASが自動でダミー変数のうち一つを除外したようです。今回は問題にすることではないので先に進みましょう。
1段階目の回帰のF値を見ても弱相関操作変数の問題はなさそうです。では、2段階目の結果へ行きましょう。

新たな操作変数を追加したことにより、決定係数が上昇しています。また、S の係数も、0.16 へと上昇しています。つまり、1年の学習が16%の収入上昇につながることになります。最初の結果が10%程度という結果を出したことと比較すると、かなり開きが有ります。
今回、内生変数2つに対し、操作変数を3つ使用しました(厳密には、他の説明変数も操作変数なので、「操作変数が新たに3つ追加された」というのが正しいです)。内生変数より多くの操作変数を用いた場合、過剰識別の問題が発生します。操作変数と内生変数の数が一致する場合にはない問題です。操作変数の数が内生変数より多くなっても, 操作変数としての条件を満たしているかという問題で, これを検定するのが、 Sargan の J 検定と呼ばれるものです。これも例によって自力で計算しましょう。

/* Sargan J stat */
/* 変数の数と残差を得るために再実行 */
proc reg data=work.grilic outest=work.sargan_inst(keep= _P_) RSQUARE noprint;
 first_stage: model S IQ = EXPR TENURE RNS SMSA y_: MED KWW AGE;
 output out=work.grilic2 predicted=Shat IQhat residual=S_r IQ_r;
quit;
proc reg data=work.grilic2 outest=work.sargan_endo(keep=_P_) RSQUARE noprint;
 TSLS: model LW = Shat EXPR TENURE IQhat RNS SMSA y_:;
 output out=work.sargan residual=resid;
quit;
/* R2 を計算*/
proc reg data=work.sargan outest=work.sargan_R(keep=_EDF_ _P_ _RSQ_) rsquare noprint;
 Sargan: model resid = EXPR TENURE RNS SMSA y_: MED KWW AGE;
quit;
/* 検定量とp値を表示 */
data _null_;
 merge sargan_R
    sargan_inst(rename=(_P_ = _P_inst_))
    sargan_endo(rename=(_P_ = _P_endo_));
 N=_EDF_ + _P_;
 Sargan=N*_rsq_;
 dof=_P_inst_ - _P_endo_;
 p_value=cdf('CHISQUARE',Sargan, dof);
 R2 = _RSQ_;
 put N= R2= Sargan= p_value= dof=;
 stop;
run;

Sargan の検定は尤度比検定の一種で、2段階目の回帰分析の残差に対して操作変数で回帰したモデルの決定係数×サンプルサイズが、帰無仮説のもとで、自由度が (追加した操作変数の数) - (内生変数の数) のカイ二乗分布に従うことが知られています。ただし、帰無仮説は「過剰識別制約を満たす」なので、操作変数がうまく機能している時、このテストは棄却できないことになります。
回帰したモデルのパラメータの推定値をデータセットに出力するには、regプロシジャに outest= オプションを追加します。さらに、RSQURE オプションを追加することで決定係数などの指標も出力できます。しかし、outest= では推定に用いたオブザベーションの数を出力できません。データセットのオブザベーション数を求めるには、contentsプロシジャが使えますが、欠損値がある場合少々面倒です(このデータセットでは欠損値はありません)。そこで、ちょっとしたテクニックなのですが、上記のコードでは、定数項を含む説明変数の数 (_P_) + 誤差項の自由度 (_EDF_) がオブザベーション数に一致することを利用して計算してます (後から見つけたのですが、公式でもこのテクニックは紹介されていました:REGプロシジャで解析に使用されたオブザベーション数をSASデータセットへ出力する方法)。
各プロシジャの結果をデータセットに保存し、最後のデータステップではそれらを読み込んでサンプルサイズ(N)、決定係数 (R2)、Sargan検定量 (Sargan)、p値 (p_value)、カイ二乗分布の自由度 (dof)を計算しています。
今回は、単なる線形回帰でしたが、実際の研究ではモデルの当てはまりを良くするためにいろいろな形状の式が考えられています(この辺りの話は Bazen の "Econometric Method for Labour Economics" にまとめて書かれています。関心を持たれたら是非読んでみてください。)。

補足: 賃金関数の研究について

この節では、操作変数法のテクニカルな話から離れて、賃金関数の研究について補足説明をします。
今回は操作変数法の説明をテーマにしたため、賃金関数の説明が粗っぽくなっています。ですので、これを読んだからと言って「じゃあさっそく会社をやめて大学に入り直せば給料が上がるじゃないか」などと考えてはいけません(笑)。なぜなら、これはあくまで「平均的な」教育の効果を表しているだけだからです。
「高学歴ほど給料が上がる」は確かに事実ですが、実はこれとは別に「高学歴ほど収入の振れ幅が大きい」という研究も存在します。つまり平均して高学歴ほど給料が高くなりますが、一方で給料が低くなる「リスク」も大きくなります。また、世の中には様々な業種があり、全く異なる業種に簡単に移動するのは難しいでしょう。
加えて、今回のデータは男性だけのデータです。女性の場合、結婚を機に仕事をやめる、という人が男性より多いため、そもそも賃金を記録できません。よって、このままでは就業していない女性を完全に無視した分析となってしまいます。さらにアメリカの場合、白人と有色人種、特にアフリカ系の格差が昔から社会問題となっていることもあり、性別や人種の差というのかなり以前から研究者の問題意識として存在したように思われます。
現在は、このような仕事の「質の違い」を考慮したり、「平均」以外の切り口から見たりする研究も進んでいます。最近話題になっているトマ・ピケティの著作「21世紀の資本」も、そのようなここ最近(といっても10年スパン程度の話ですが)の「格差の研究」をまとめたものとなっているようです (私は残念ながらまだ読んでいません)。

データの出処

実は、これは計量経済学の教科書である、 Fumio Hayashi の "Econometrics" でも取り上げられている問題で、林先生のWebサイトからダウンロードすることが出来ます。それ以外の海外の大学のサイトでも教材としてダウンロードできるようにしているところがあります。

• https://sites.google.com/site/fumiohayashi/home/links-related-to-hayashi-econometrics/data-for-empirical-exercises 林先生の個人サイト。 grilic.xls が今回使用したデータです。
• https://ideas.repec.org/p/boc/bocins/griliches76.html こちらからダウンロードできるのは STATA というソフトウェアで使用される .dta 形式なので、冒頭の import プロシジャのオプションを dbms=dta に設定し、getnames=を削除する必要があります:

  
  proc import datafile='/folders/myfolders/griliches76.dta'
     out     =work.grilic 
     dbms=dta replace;
  run;
  

  中身はエクセルファイルと同じですが、こちらをインポートした場合、変数にラベルが付きます。

参考文献

• SAS公式 syslin プロシジャリファレンス: http://support.sas.com/documentation/cdl/en/etsug/60372/HTML/default/viewer.htm#etsug_syslin_sect007.htm
• SAS公式 output ステートメント: http://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_reg_sect015.htm
• SAS公式 reg プロシジャリファレンス
• Hayashi, Fumio [2000] "Econometrics," Princeton University Press
• Blackburn, McKinley and David Neumark [1992] "Unobserved Ability、Efficiency Wages、and Interindustry Wage Differentials," Quarterly Journal of Economics、Vol. 107, pp. 1421-1436
• Bazen, Stephen [2011] "Econometric Methods for Labour Economics," Oxford University Press
• Wooldridge, Jeffrey M. [2010] "Econometric Analysis of Cross Section and Panel Data," 2nd Ed. MIT Press