2015年10月19日月曜日

SAS での回帰不連続デザイン (RDD) とビジネスでの応用について

2011年の米南東部 SAS ユーザー会 (SESUG) 発表資料で政策の効果を推定するための回帰不連続デザイン (RDD)を実行するマクロが紹介されている発表 Schoeneberger (2011) "RDPLOT: A SAS® MACRO for Generating Regression Discontinuity Plots" を見つけたので、
  1. RDD マクロの解説
  2. RDD そのものの解説
  3. ビジネスでの応用
をあわせて紹介したいと思います。かなりマニアックなテーマなので、「SAS にこじつけて自分の趣味の話をしたいだけだろ」と思われるかも知れません。おおむねその通りですがどうかご容赦ください。

RDDとは?

回帰不連続デザイン (RDD, Regression Discontinuity Design) というのは回帰分析の一種で、初め Thistlethwaite & Campbell (1960) "Regression-discontinuity analysis: An alternative to the ex post facto experiment" で発表されたものです。社会科学の研究では割りとポピュラーな手法ですが、日本語での情報はあまりありません。そのため, RDD の和訳も定着しておらず、切断回帰デザインとか非連続回帰計画というふうに呼ばれることもあります。例えば日本語ですと
  • アングリスト&ピスケ (2013) 『ほとんど無害な計量経済学 –応用経済学のための実証分析ガイド–』, 大森義明・田中隆一・野口晴子・小原美紀訳, NTT出版. 原著: Angrist, Joshua D. and Jrn-Steffen Pischke (2008) “Mostly Harmless Econometrics: An Empiricist’s companion,” Princeton University Press
  • 津川友介 (2015) 『回帰分断デザインRegression discontinuity design(および分割時系列デザイン)|医療政策学×医療経済学
  • 森田果 (2014) 『実証分析入門 –データから「因果関係」を読み解く作法–』, 日本評論社
などで RDD が言及されています。より専門的な解説をした英語の論文は末尾の参考文献に挙げておきます。私は学生時代に、教育政策の効果を検証する方法として学びました。では、ビジネスの文脈、例えばマーケティングなどに応用されてはいないかと調べたところ、 Marketing Science 誌に掲載された論文で、RDD をマーケティング全般へ応用する方法について論じた Hartmann, Nair, & Narayanan (2011) "Identifying Causal Marketing Mix Effects Using a Regression Discontinuity Design" を見つけました。
上記の文献を読むだけでもなにを意図したものかおおむねわかるかと思いますが、この Hartmann らの論文に沿って簡単に解説したいと思います。
まず、ある政策を実施して、その効果がどれだけあったかを知りたいとします。 民間企業の立場なら、顧客に対するプロモーションを実施して、どの程度売上につながったか、とも解釈できます。どちらの想定でも、何らかの成果を期待してコストを払って施策を行っているわけなので、実際に効果があがっているかの検証は大事です。Hartmann et al. (2011) では、カジノでは少数の上客 (ハイローラー) が大きな利益をもたらしてくれることから、ハイローラーに対する優遇サービスを例に挙げていますが、 他の業種でのプロモーションの施策に置き換えて考えても問題ありません。単純に考えると、ハイローラーに対し優遇サービスを実施し、それ以外の利用客には実施しないとすると、
(効果) = (優遇を受けた利用客の平均) - (優遇を受けなかった利用客の平均)
というふうに優遇を受けた利用客と受けなかった利用客で分割し、それぞれのグループについて1人あたり売上額などを計算し、さらにその差をとれば、効果の大きさとみなせそうです。視覚的にあらわすならば、以下の図1のような棒グラフで表せます。
図1: 施策の有無と1人あたり売上額の差

ところが、実際は効果の大きさを正しく見積もるのは容易ではありません。2つのグループの利用客の傾向がかたよるからです。そもそも、カジノ側は、最初から多く遊んでくれるハイローラーに対してのみ優遇します。そのためハイローラーは、仮に施策を受けなかったとしても多く遊ぶ可能性が高いのです。逆に優遇を受けなかったグループは、優遇を受けたとしても商品を購入する可能性が比較的低いかもしれません。よって、ここで計算した差は、純粋な施策の効果だけでなく、偏った平均によるバイアスも含まれてしまい、過大に見積もってしまいます。同様の理屈で、逆にあまり遊ばない利用客に対して重点的に優遇を行った場合、効果が過少評価されるでしょう。
医薬品の臨床実験や、農作物の品種試験でも、薬品を投与した場合・しない場合や、新しい品種・既存の品種を比較して性能の差を見るのですが、このような偏りが起こらないようにコントロールするのが普通です (いわゆるランダム化)。 しかし、マーケティングに関しては利益に直結するため必ずしもこのコントロールができるとは限りません。
このように、偏りがある場合の対処方法として考案された方法が RDD です。 特に、この施策を受けるか受けないかの判断基準がしきい値 (これをカットオフ点ともいいます) で表現されるとき RDD が利用可能です。 ある数値を超えた人間が施策を受け、超えなかった人間は施策を受けないとします。
カジノの例でいうならば、利用客ごとの1日にカジノで使う金額を見て、一定額以上ならハイローラーに分類し、優待サービスを実施します。事実、ハイローラーに分類された利用客は1日あたりの使用金額が非常に多く、かつ来場回数も増えることが確認されています。利用客のハイローラーへの分類の前と後の利用金額の関係を擬似的に表すと、以下の図2のようになります。
図2: RDDの視覚的な説明

実際に効果があるのなら、中央のカットオフ点を境に、直線が断層のようにずれるはずです。この直線の、不連続な地点の上下のずれの大きさが、施策による効果の大きさとみなせます。
ここで、このカットオフ点の前後、たとえば前後1万円とすると、ハイローラーに分類されるかどうかに関係なく、適用前の利用客の使用金額にほとんど違いがありません。よって、この限られた範囲ならば、先に述べた偏りがほとんどないと考えられます。ここから、データから狭い範囲に属する利用客の情報だけを取り出して分析すれば、効果の大きさをうまく見積もることができる、というのが RDD の基本理念です。RDD は回帰分析を応用することでこのずれの大きさを計算し、そこから効果の大きさを推定します。

SAS での実行方法

冒頭で紹介した Schoeneberger (2011) で、 SAS マクロのすべてが公開されています。これは SAS/STAT SAS/GRAPH があれば実行可能なのでたいていの環境で実行できるのですが、 SAS University Edition では SAS/GRAPH がないため gplot プロシジャが使用できません。 そこで、 sgplot プロシジャで でグラフを描画するように書き換えたコードを使用します。このプログラムは記事の末尾に掲載しておきます。 元の regdis マクロはインプットデータを物理名で参照しますが、このままだと使い勝手が悪いので、 そのまま使う場合でも、元のプログラムも9ページ右上の最初のデータステップの
set "%bquote(&path)%bquote(&data_in)";
のところを
set &data_in;
に変えたほうが使いやすいと思います。
次に、以下のコードでインプットデータを作成します。
/*------------------------------
 RDD デモコード
 SAS/STAT が必要
-------------------------------*/

data param;
input _TYPE_$ _NAME_$ x_prior x_post;
cards;
COV x_prior 2 1 
COV x_post  1 3
;
run;

proc simnormal seed=1 numreal=100 data=param(type=cov) out=sample ;
    var x_prior x_post;
run;

data sample;
    call streaminit(1);
    set sample;
    x_prior+10;
    x_post +10;
    is_high=0;
    if x_prior >10 then do;
        x_post+rand('normal',4,.5);
        is_high=1;
    end;
run;
simnormal プロシジャというのは、与えられたパラメータから乱数を生成するプロシジャで、 SAS/STAT があれば使用可能です。これで2次元正規乱数を生成し、さらに加工することで擬似データを作成しています。乱数の種を固定しているので、毎回同じ結果になるはずです。データは各レコードが利用客の1人分の記録を表し、変数の x_prior が施策前の利用客の使用金額、 x_post が施策後の金額、 is_high が適用対象かどうかを表す変数です。
このデータを regdis マクロに入力します。私がデータセットの1行だけ改変したものですと、
%regdis(%str(),work.sample,x_prior,x_post,is_high,10,above,Y);
後で掲載する University Edition 対応版ですと
%regdis(work.sample,x_prior,x_post,is_high,10,above,Y);
を実行してください。以下は引数の解説です。
  1. path : ファイルのあるフォルダパスです。UniversityEdition版は廃止したため、次の data_in が第一引数になります。
  2. data_in: データセットのファイル名です。マクロを改変したので今回はデータセット名を書いています。
  3. pretest: 施策前の状態を表す変数です。X軸に当たります。
  4. posttest: 施策後の状態を表す変数です。Y軸に当たります。
  5. group: 各個人に対して施策がなされたかどうかの数値変数。適用対象なら 1、そうでなければ 0 です。
  6. cutscore: 適用対象を決めるしきい値 (カットオフ点)。今回は 10 とします。
  7. treat_loc: カットオフ点のどちら側が施策適用対象か。今回は 10 を越える利用客が対象なので、 above と入力します。逆は below です。
  8. plot_data: グラフにレコードの散布図を描くかどうか。 Y/N どちらかを入力します。 基本的に Y のほうが良いでしょう。
これを実施すると、 glm プロシジャで直線の傾きと切片が推定されるほか、結果のグラフが以下の図3のように表示されます (SAS University Edition の場合、使用しているスタイルによって結果が若干変わります)。
図3: 施策前後の使用金額の散布図と RDD による回帰直線の当てはめ

このマクロでは、直線回帰のみですが、仮説検定から2つのグループでそれぞれ傾きが異なると判定されたため、それぞれで傾きが異なっています。入力データによっては、両者の傾きが同じになる場合もあります。実際に RDD を使用する場合には、他のモデルと比較したり、使用するデータの範囲を決めたりといったことが必要です。

応用の可能性

Hartmann et al. (2011) ではさらに、いろいろな分野での応用の可能性を指摘しています. Table 6. にその一覧表があるので、和訳したものを掲載します。ここでのスコア変数というのは、上記のグラフの横軸に対応する変数です。一般には RDD は縦と横で異なる変数をとることができます。いずれの場合でも、スコア変数があるカットオフ点を超えた時に施策が行われるという例を示しているので、 RDD のフレームワークを用いて施策の有無と結果変数の差を見ることで、施策の効果の大きさがわかるということです。

SAS University Edition 対応版マクロ

/*---------------------------------------------------------------------
RDPlot macro program modified to use it in SAS University Edition
S-Katagiri, Nature Insight co., 24/Jun/2015
http://blog.n-insight.co.jp/

The original macro program is from
Schoeneberger, J. A. (2011).
RDPLOT: A SASR Macro for Generating Regression Discontinuity Plots. 
Proceedings of the annual Southeastern SAS Users Group Conference.
http://analytics.ncsu.edu/?page_id=3375
---------------------------------------------------------------------*/

**assumes data_in file is one line per unit-of-analysis record;
**assumes group is dichotomous, numeric variable coded 0=control, 1=treatment;
*data_in=the SAS dataset name for analysis;
*pretest=the pre-test score by which a cut-score being applied creates treatment-control groups;
*posttest=the posttest, dependent variable of interest;
*group=variable denoting treatment (1) or control (0) group membership;
*cutscore=numeric value on the pre-test that serves as the cut-score;
*treat_loc=denotes whether the treatment is above or below cut-score: insert 'above' or 'below' for this value;
*plot_data=denotes whether the user wants to plot actual data 'Y' or just RD lines 'N';

%macro regdis(data_in,pretest,posttest,group,cutscore,treat_loc,plotdata);
    data regdis;
        set &data_in.;
        ** homogeneity の検定のために相互効果の変数作成;
        interact=&group*&pretest;
    run;
    ** x, yの最小値最大値をマクロ変数に取得;
    proc sql noprint;
        select min(&pretest), max(&pretest), min(&posttest),
        max(&posttest)
        into: x_min, : x_max, : y_min, : y_max
        from regdis;
    quit;
    ** teart_loc に応じて補間点を並べる;
    %if "%upcase(&treat_loc)"="BELOW" %then %do;
        %let x_0=&x_max;
        %let x_1=&x_min;
    %end;
    %if "%upcase(&treat_loc)"="ABOVE" %then %do;
        %let x_0=&x_min;
        %let x_1=&x_max;
    %end;
    ** glmによる相互作用テスト (homogeneity assumption);
    ods output parameterestimates=paramsint;
    proc glm data=regdis namelen=32;
        model &posttest=&group interact &pretest/ss1
        solution;
        run;
    quit;
    ** homogeneity なら傾きを一定とする。そうでないなら個別に傾きを計算;
    data paramsint;
        set paramsint;
        if (parameter="interact") and (probt gt .05) then
            call symput('hetero','N');
        if (parameter="interact") and (probt le .05) then
            call symput('hetero','Y');
    run;
    ** 以下は傾き一定の場合 (hetero='N');
    %if "&hetero"="N" %then %do;
        ** 傾きを計算。;
        ods output parameterestimates=params;
        proc glm data=regdis namelen=32;
            model &posttest=&group &pretest/ss1 solution;
            run;
        quit;
        ** パラメータをマクロ変数に取得;
        proc sql noprint;
            select estimate into: intercept_g0
            from params
            where parameter="Intercept";
            select estimate into: grp_slp
            from params
            where parameter=lowcase("&group");
            select estimate into: cut_slp_g0
            from params
            where parameter=lowcase("&pretest");
        quit;
        %let intercept_g1=%sysevalf(&intercept_g0+&grp_slp);
        %let cut_slp_g1  =&cut_slp_g0;

    %end;
    ** 以下 homogeneity ではない場合 (hetero='Y');
    %if "&hetero"="Y" %then %do;
        ** 以下の glm で2つのグループの回帰係数を別々に推定;
        ods output parameterestimates=paramsg1;
        ** treatment グループの係数;
        proc glm data=regdis;
            model &posttest=&pretest/ss1 solution;
            where (&group=1);
            run;
            quit;
        ods output parameterestimates=paramsg0;
        ** control グループの係数;
        proc glm data=regdis;
            model &posttest=&pretest/ss1 solution;
            where (&group=0);
            run;
        quit;
        ** sql プロシジャでパラメータを取得;
        proc sql noprint;
            select estimate into: intercept_g1
            from paramsg1
            where parameter="Intercept";
            select estimate into: cut_slp_g1
            from paramsg1
            where parameter=lowcase("&pretest");
            select estimate into: intercept_g0
            from paramsg0
            where parameter="Intercept";
            select estimate into: cut_slp_g0
            from paramsg0
            where parameter=lowcase("&pretest");
        quit;
    %end;
    ** 以下、共通処理;

    ** 回帰直線の点を追加;
    data regdis2;
        set regdis end=eof;
        output;
        if eof then do;
            %do i=0% to 1;
                x=&&x_&i;
                y=&&intercept_g&i + &&cut_slp_g&i* &&x_&i;
                &group=&i;
                output;
                x=&cutscore;
                y=&&intercept_g&i + &&cut_slp_g&i* &cutscore;
                &group=&i;
                output;
            %end;
        end;
    run;
    data regdis2;
        set regdis2;
        if &group=1 then group='treatment';
        else group='control';
    run;
    ** スタイルの定義;
    proc template;
    define style Styles.rdplot;
        parent = Styles.Listing;
        style GraphData1 from GraphData1 / 
            markersymbol = "circle"
            color = blue 
            ;
        style GraphData2 from GraphData2 / 
            markersymbol = "triangle"
            color = red 
            ;
        end;
    run;
    ** グラフ出力;
    ods html path='/folders/myfolders' file='rdplot.html' style=Styles.Rdplot;    
    ** グラフ描画;
    proc sgplot data=regdis2;
            refline &cutscore / axis=x;
            scatter x=&pretest y=&posttest /name='plot' group=group %if "%upcase(&plotdata)"="N" %then markerattrs=(size=0);;
            series x=x y=y /group=group;
            xaxis label="pretest"  minor;
            yaxis label="posttest" minor;
            keylegend 'plot' /noborder across=2 title='Intervention Group';
    run;
    quit;
    ods html close;

%mend regdis;

参考文献

Hahn, J., Todd, P., & Van der Klaauw, W. (2001). Identification and Estimation of Treatment Effects with a Regression-Discontinuity Design. Econometrica, 69(1), 201–209. doi:10.1111/1468-0262.00183
Hartmann, W., Nair, H. S., & Narayanan, S. (2011). Identifying Causal Marketing Mix Effects Using a Regression Discontinuity Design. Marketing Science. doi:10.1287/mksc.1110.0670
Imbens, G. W., & Kalyanaraman, K. (2012). Optimal Bandwidth Choice for the Regression Discontinuity Estimator. Review of Economic Studies, 79(3), 933–959. doi:10.1093/restud/rdr043
Imbens, G. W., & Lemieux, T. (2008). Regression discontinuity designs: A guide to practice. Journal of Econometrics, 142(2), 615–635. doi:10.1016/j.jeconom.2007.05.001
Lee, D. S. (2008). Randomized experiments from non-random selection in U.S. House elections. Journal of Econometrics, 142(2), 675–697. doi:10.1016/j.jeconom.2007.05.004
McCrary, J. (2008). Manipulation of the running variable in the regression discontinuity design: A density test. Journal of Econometrics, 142(2), 698–714. doi:10.1016/j.jeconom.2007.05.005
Schoeneberger, J. (2011). RDPLOT: A SAS® MACRO for Generating Regression Discontinuity Plots. In SESUG Proceedings 2011. URL: http://analytics.ncsu.edu/?page_id=3375
Thistlethwaite, D. L., & Campbell, D. T. (1960). Regression-discontinuity analysis: An alternative to the ex post facto experiment. Journal of Educational Psychology. doi:10.1037/h0044319
2015/7/3 片桐

0 件のコメント:

コメントを投稿

ツイート数からみる"バーチャルYouTuber"ブーム

今や YouTuber の話題の半分を占めるほどのクチコミ数に 当社が提供するソーシャルビッグデータ検索ツールの「 beInsight (ビーインサイト)」を使って、話題の「バーチャル YouTuber 」について調べてみました。 「バーチャル YouTuber...