SASによる解析：操作変数法

今回は以前に書いた操作変数法に関する概略の続きとして、実際に操作変数法を用いて因果効果の推定が行われた研究を、後追いしつつ、紹介していこうと思います。前記事はこちら。

norihirosuzuki.hatenablog.com

使用するデータはangrist and krueger(1991)での教育年数と賃金に関するデータです。この研究は、操作変数法を利用した研究としては経済分野では非常に有名で、様々な文献、書籍等でその名前があげられています。

操作変数法の概要（再掲）

本題に入る前に、まずは操作変数法（IV estimation）という因果推論の方法論についての復習をしようかと思います。より詳細に、厳密には前記事や専門書などを参考にしていただければいいかなとは思いますが、ざっくりと理解する分にはこの記事だけでも十分かなと思います。なお説明の際にはDAGsも使っているので念のため、DAGの紹介記事だけおいておきます。とりあえずこの記事の中では矢印があれば、矢印の元が原因、先が結果になっている関係だと、ひとまず理解しておいてもらえればいいかなと思います。

norihirosuzuki.hatenablog.com

また、操作変数については知ってるから実際に研究した例を見たいという方は、この項目を飛ばして分析例紹介の項目までスキップしてください。

操作変数とは

操作変数法は、ある興味のある介入、曝露（A）がアウトカム（Y）に対して与える因果効果を算出するにあたって、操作変数（Instrumental Variable）と呼ばれる変数を利用します。ある変数Zが操作変数となるためには以下の3つの条件が必要です。（より正確には後述の2項目のうちいずれか一方も追加で必要）

1. Aと関連性を持つ（relevance condition）
2. Aを通してのみYに効果を与える（除外制約, exclusion restriction）
3. Yと共通の原因（Common cause）を持たない

この3つの条件についてそれぞれ簡単に説明を行います。

Aと関連性を持つ（relevance condition）

一つ目の条件は、操作変数Zと、曝露（介入）を示す変数Aが関連性を持つというものです。これはZがAに対して因果効果を与えている必要はなく、もっと緩い仮定であるZがAと関連していること、つまりCov(Z,A)≠0ということ、もしくはＺとＡが独立でないことを意味しています。この条件については、操作変数が1つである場合にはt検定で真値が0であるか、複数である場合にはF検定によって全変動が0であるかを帰無仮説と置くことで、実際のデータから検証が可能です。

ただ相関していれば（帰無仮説が棄却されていれば）なんでもいいわけではなく、その相関の”強さ”にも注意する必要があります。おおよそですが、t統計量が3.2、F統計量であれば10を超えるというのが目安とされています。

Aを通してのみYに効果を与える（除外制約, exclusion restriction）

二つ目は除外制約、exclusion restrictionと呼ばれるものです。only through conditionとか唯一経路条件とかとも呼ばれるかなと。これは操作変数ZはAを通した経路以外でYに対して効果を与えるということを指していて、DAGsでは下図のようにＺからＹに対して直接の矢印が存在しないことを意味します。

これを数式で表現すると、すべての個人i、操作変数ｚ, z'、介入aに対して

が成り立つことを意味します。これは個人レベルでの除外制約を表現したものですが、この仮定がもう少し緩まったpopulationレベルでの話（期待値をとったもの）の制約を論じているものもあります。

ただこれは一つ目の条件とは異なり、検証が不可能な仮定です。というよりも一つ目の条件だけが検証できるといったほうが伝わるかなと思います。ちなみにこの項目については、二重盲検が行われる研究だと成り立つことが予想されています。

Yと共通の原因（Common cause）を持たない

3つ目の条件は、操作変数ZとアウトカムYが共通の原因を持たないこと、つまり操作変数Zが、介入AとアウトカムYの交絡因子Uと関連性を持たない（独立である）というものです。

実はこの条件が、回帰モデルや傾向スコアを用いた諸手法などと大きく差別化が図れる点です。操作変数以外の手法だと、因果効果の推定のためにはこのAとYの共変量Uはすべて測定される必要があります。もし仮に観測されない共変量が存在してしまう場合には、推定量自体に一定のバイアスが含まれてしまいます。しかし操作変数法では検証できないにしろ、共変量Uと操作変数Zが独立であるのならば、Uはすべて測定する必要なく因果効果を推定することができます。こういった点から操作変数法は共変量が測定できない状況などで使われているわけです。

追加の項目

上記の3項目が世間一般でよく操作変数の仮定として挙げられますが、より厳密には以下の2項目のうちのどちらか一方が追加で必要な仮定です。

均一性（homogeneity）

均一性（homogeneity）は一般的には、治療群と対象群の平均因果効果（ATE）がZに関わらず等しいというものです。この条件を追加することによってIV estimandが平均因果効果に等しいことが証明だれます。もし証明が気になる方がいれば、参考の項目にあるWhat IfのTechnical Point16.2に記載があるのでそちらも是非ご参照ください。

ただ、この均一性という条件は現実的には保証がされないケースがしばしばあるため、その利用に関しては批判もされています。そのため解決策として考えられているのが

①ベースラインでの共変量をモデルに組み込む

②均一性（homogeneity）ではなく、単調性（monotonicity）を使う

というものです。

単調性（monotonicity）

単調性を考えるにあたって、まず最初に割り当て（Z）と、治療（A）の潜在アウトカムを考えます。これはAとYの潜在アウトカムを考えたときと同様の概念です。（Z,Aはバイナリー）

すると、この潜在アウトカムの取りうるパターンによって研究集団は下図のように四種類に分類することが可能です。

ここで単調性は、上図の右下にあるDefierがいないという内容を指した仮定になります。この単調性は、均一性に比べてより現実的なものですが、その利用には注意が必要です。先ほどの均一性の仮定を使った場合には、IV estimandは平均因果効果（ATE）を推定していることになりますが、単調性を使った場合には、Complier集団での因果効果を推定していることになります。

じゃあこのComplier集団での因果効果、つまり局所的平均因果効果（LATE）がATEであるのかというと基本的には異なります。これら二つが一致するのはあくまでも集団全体がComplierであるという場合に限り、その時のみ、ATE=ATE=ITT（Intention to Treat）となるわけです。

推定量に関して

操作変数法ではZとAが二値である場合、IV estimandを次のような式で算出し、上記の諸仮定が満たされる場合には平均因果効果としてみなすことが可能です。

なおここではZとAは二値変数であるとして考えていますが、連続変数である場合には、期待値ではなく共分散で考えればOKです。

この式のイメージとしてはこんな感じ。

1. 操作変数ZはAを介してのみYに対して影響を与える。（①＋②）
2. ZはAに対して影響を与える。（①）
3. 興味ある対象は②であるから、（①＋②）/ ①とすることで②を算出

また、２段階最小二乗法を用いてIV estimandを算出することも可能です。

ここにまとめようかと思ったのですが、わかりやすくまとめていらっしゃる方がいたので、そちらの資料を紹介します。

https://www.nishiyama.kier.kyoto-u.ac.jp/2017/jugyochukei3.pdf

分析例紹介

今回の紹介する研究は、angrist and kruegerらによって1991年に行われた
「DOES COMPULSORY SCHOOL ATTENDANCE AFFECT SCHOOLING AND EARNINGS?」（義務教育への出席率は学校教育と収入に影響を与えるか？）というものです。

論文の概要であったり、考察に関しては今回の記事のテーマである「SASでの操作変数法の実装」という観点から離れてしまうのでそこは触れていません。もし興味がある方がいましたら、参考の欄から適宜元論文などに飛んでもらえればいいかなと思います。

使用データ説明

使用データは上記の通り、angrist and krueger（1991）らの研究データです。元データについては、マサチューセッツ工科大学のサイトで公開されていますので、こちらからダウンロードが可能です。今回の解析で使用しているのは必要最低限のデータが含まれているasciiqob.zipに含まれているものです。

データの内訳としては、サンプルサイズが329,509、項目は以下の通りです。

• lwklywge（log weekly wage

週当たりの賃金の対数を取った値

• educ（education）

最高学年修了

• yob（year of birth）

生年月日

• qob（quater of birth）

一年間を４つに分けたとき（四半期）にいつ生まれたか

• pob（place of birth）

出生地（1980年の国勢調査の州コード）

解析内容

操作変数法の利用にあたって操作変数とするのは、QOB（どの時期に生まれたか）です。この変数を操作変数とする妥当性は次の観点から得られています。

1. アメリカでは６歳になる年に義務教育が始まり、１６歳になるまで退学が不可能である。このことから誕生した月日によって各個人が受ける教育の期間が異なり、QOBとeducaionの間には関連がある。
2. 各個人がいつ生まれるのかと、個人の能力には直接的な関係はない。

以下はAngrist and Krueger（2001）より引用した、出生時期と教育年数の関係を示した図であり、正の相関があることがうかがえます。

とういうわけで、出生時期を操作変数とした場合の解析を以後見ていきます。解析に当たってはSASのsyslinプロシジャを用います。このプロシジャは線形回帰モデルを同時に実行することが可能であり、後ほど出てきますが操作変数法にあたり２段階回帰を行うのですが、その作業を一つのプロシジャで完結させることが可能です。構文、オプションについてはこちらのヘルプを参考にどうぞ。

documentation.sas.com

angrist and krueger（1991）では操作変数をQOB, YOBを利用しつつ、年齢やいくつかのダミー変数を加えたモデルを考えていますが、基本的な考え方としてはこの記事の序盤でまとめていたものと変わりません。

まずは単純に教育を説明変数、目的変数を賃金とした場合の重回帰モデルを考えます。なお yr21-yr29は下記コード中のマクロで作成した操作変数です。おそらくですがこれらの変数を説明変数に入れた理由としては、それぞれの操作変数の教育を通しての間接効果を除いた効果（直接効果）がどの程度かを測るためなのかなと個人的には理解しています。（つまり０になればOK）

この回帰モデル（OLS, Ordinary Least Squares regression ）からは、教育の偏回帰係数は0.071081と推定がされました。つまり教育年数が１年増加すれば、おおよそ賃金の対数値が0.07増加するといった関係性がうかがえるということですが、ここには交絡による影響が含まれていると推測がされます。

次に操作変数法として、２段階最小二乗法（TSLS; two-stage least squares）によるIV estimatorの推定を行います。ステップとしては次の流れになります。

1. 操作変数Zから介入A（今回は教育）への回帰モデルを考えOLS推定を行う。この回帰によってAの推定値を得る
2. 1の操作で得たＡの推定値を用いて、説明変数をアウトカムＹ（対数賃金）とした回帰モデルを考えOLS推定を行う

この２段階回帰モデルはSASのsyslinプロシジャを用いて簡単に行うことが可能です。今回のセッティングで得られた結果はこちら。

ここでeducの係数推定量は0.089115となり、先ほどの単純な重回帰モデルよりややその値が増加したことが分かります。

なおこの研究の問題点としては、今回用いた操作変数が弱い操作変数であることがあげられていますが、結果の解釈も含め、今回はその部分に関しては取り扱いません。

コード

参考の欄にあるMITのサイトよりTextファイルをダウンロード後、テキストファイルの保存先のパス名をコピーして、下記コード中の赤文字の部分に入れてください。

MITのサイトで ak91.sasという読み取り、解析用のファイルが公開されていますが冒頭の部分でエラーが起きたので、一部変更したものが下記コードになります。またメモリに関するエラーが出る場合がありますが、動いてる分には特に問題ありません。

data zero;
infile "テキストファイルのパス名";
input lwklywge educ yob qob pob;
run;

proc means;
title '1980 qob extract';
data one;
 set zero;

%MACRO QTRSYR;
             QTR220 QTR320 QTR420        QTR221 QTR321 QTR421
             QTR222 QTR322 QTR422
             QTR223 QTR323 QTR423        QTR224 QTR324 QTR424
             QTR225 QTR325 QTR425
             QTR226 QTR326 QTR426        QTR227 QTR327 QTR427
             QTR228 QTR328 QTR428        QTR229 QTR329 QTR429
%mend;

LENGTH QTR120 3
       %QTRSYR 3
       YR20-YR29 3
       QTR1-QTR4 3;

ARRAY AQTR QTRSYR;
  DO OVER AQTR; AQTR = 0 ; END;

YR20 = *1;
YR21 = *2;
YR22 = *3;
YR23 = *4;
YR24 = *5;
YR25 = *6;
YR26 = *7;
YR27 = *8;
YR28 = *9;
YR29 = *10;

QTR1 = (QOB=1);
QTR2 = (QOB=2);
QTR3 = (QOB=3);
QTR4 = (QOB=4);

QTR120 = QTR1 * YR20; QTR220 = QTR2 * YR20; QTR320 = QTR3 * YR20;
QTR420 = QTR4 * YR20;
QTR121 = QTR1 * YR21; QTR221 = QTR2 * YR21; QTR321 = QTR3 * YR21;
QTR421 = QTR4 * YR21;
QTR122 = QTR1 * YR22; QTR222 = QTR2 * YR22; QTR322 = QTR3 * YR22;
QTR422 = QTR4 * YR22;
QTR123 = QTR1 * YR23; QTR223 = QTR2 * YR23; QTR323 = QTR3 * YR23;
QTR423 = QTR4 * YR23;
QTR124 = QTR1 * YR24; QTR224 = QTR2 * YR24; QTR324 = QTR3 * YR24;
QTR424 = QTR4 * YR24;
QTR125 = QTR1 * YR25; QTR225 = QTR2 * YR25; QTR325 = QTR3 * YR25;
QTR425 = QTR4 * YR25;
QTR126 = QTR1 * YR26; QTR226 = QTR2 * YR26; QTR326 = QTR3 * YR26;
QTR426 = QTR4 * YR26;
QTR127 = QTR1 * YR27; QTR227 = QTR2 * YR27; QTR327 = QTR3 * YR27;
QTR427 = QTR4 * YR27;
QTR128 = QTR1 * YR28; QTR228 = QTR2 * YR28; QTR328 = QTR3 * YR28;
QTR428 = QTR4 * YR28;
QTR129 = QTR1 * YR29; QTR229 = QTR2 * YR29; QTR329 = QTR3 * YR29;
QTR429 = QTR4 * YR29;

keep yr21-yr29 %qtrsyr educ lwklywge qtr1-qtr3;

*******************************************************;

proc means;
title 'working data set';proc syslin data=one;
title 'OLS';

model lwklywge = yr21-yr29 educ;

/* part I: 30 instrument case */

proc syslin data=one 2sls;
     title '30-instrument case: regular 2sls';
     instruments yr21-yr29 %qtrsyr;
     endogenous lwklywge educ;
     model lwklywge = yr21-yr29 educ/overid;
run;

参考

• CAUSAL INFERENCE What If. Miguel A. Hernan, James M. Robins

www.hsph.harvard.edu

• Massachusetts Institute of Technology Economics

https://economics.mit.edu/faculty/angrist/data1/data/angkru1991

• The SYSLIN Procedure - SAS Support

support.sas.com

*1:YOB=30) OR (YOB=40

*2:YOB=31) OR (YOB=41

*3:YOB=32) OR (YOB=42

*4:YOB=33) OR (YOB=43

*5:YOB=34) OR (YOB=44

*6:YOB=35) OR (YOB=45

*7:YOB=36) OR (YOB=46

*8:YOB=37) OR (YOB=47

*9:YOB=38) OR (YOB=48

*10:YOB=39) OR (YOB=49