Stataの予測値コマンド(3)
※本記事に一部加筆修正しました。（2012.7.4）
Stata v11から、予測値の計算に新たにmarginsコマンドが導入されました。marginsはほぼすべての推定モデルに対応しており、確率や限界値など、柔軟なかたちでの予測値出力を可能にしてくれます。とはいえ、扱い方に少しクセがあるので、ここで基本的な使い方を説明しておきます。以下の説明はv12のマニュアルを参考にしていますが、v11でも基本的には大丈夫だと思います。
Stataでは、予測値関連のコマンドとしてpredict、adjust、marginsの3つが用意されていますが、このうちadjustはmarginsと機能が被っており、v12ではマニュアルからも消えているようです。adjustに慣れたユーザは、v12ではmarginsを使うしかありません。
基本編
まずはデータを読み込んで回帰分析します。（マニュアルでも使われているデータです。）
. u http://www.stata-press.com/data/r12/margex, clear
. reg y i.sex i.group age, noheader
------------------------------------------------------------------------------
y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
1.sex | 18.44069 .8819175 20.91 0.000 16.71146 20.16991
|
group |
2 | 5.178636 .9584485 5.40 0.000 3.299352 7.057919
3 | 13.45907 1.286127 10.46 0.000 10.93729 15.98085
|
age | -.3298831 .0373191 -8.84 0.000 -.4030567 -.2567094
_cons | 68.63586 1.962901 34.97 0.000 64.78709 72.48463
------------------------------------------------------------------------------
ここから、
y = 68.6 + 18.4*sex + 5.2*group(2) + 13.5*group(3) + -.33*age
という関数が推定されます。Stataではpredictコマンドで個々のケースの予測値を出力してくれますが、その予測値は以上の式に、個々のケースの説明変数の値をあてはめて計算されたものです。たとえば一番上のケースだと、
y = 68.6 + 18.4*1 + 5.2*0 + 13.5*0 -.33*55
で計算されます。結果は68.9（丸めのため上記式から計算するとズレます）。同じ値は、predictコマンドで出力することができます。つまりpredictは、説明変数が個々の値をとったときの被説明変数の予測値(y-e)です。参考のため、予測値の変数を作成しておきましょう。
. predict py
marginsはこれに対して、より自由に値を固定した上で予測値を計算し、さらに予測値の平均値を計算することができます。とはいえ、引数やオプションを何もとらないと、個々のケースごとに計算された予測値を全てのケースについて平均した値が計算されます。
. margins
Predictive margins Number of obs = 3000
Model VCE : OLS
Expression : Linear prediction, predict()
------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_cons | 69.73357 .3619129 192.68 0.000 69.02423 70.4429
------------------------------------------------------------------------------同じ結果は、先ほどの個々のケースの予測値を利用して、
. su py
としても得ることができます。（この場合、回帰分析の性質から、yの予測値の平均値はyの観察値の平均値と一致します。）
では引数を取るときはどうでしょうか。marginsでは引数を
- 直接
- オプション：at( )
- オプション：over()
という3つの方法でとることができるのですが、これらの使い分けには少々説明を要します。
まず、基本的に直接引数とat( )は同じように機能します。直接引数はカテゴリカル変数しかとりません（連続変数だとエラーがでる）。さらに、任意の値を指定するのではなく、カテゴリーの値それぞれについて予測値を計算します。これに対してat( )はカテゴリカルも連続も両方、任意の値を指定できます。
ここらへん分かりにくいと思いますので、具体例で説明します。
まずはsex（性別）を直接引数にとった場合。
. margins sex
Predictive margins Number of obs = 3000
Model VCE : OLS
Expression : Linear prediction, predict()
------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
sex |
0 | 60.50093 .5709154 105.97 0.000 59.38196 61.6199
1 | 78.94162 .5700065 138.49 0.000 77.82442 80.05881
------------------------------------------------------------------------------60.5は、「sexが0（女性とします）、その他の値は個々のケースの観察値としたときに計算した予測値を、すべてのケースで平均した」値です。つまり、sexが1（男性）のケースでも、それを0（女性）として予測値を計算した上で、全ケースで平均値を出しているわけです。また固定されているのはsexのみで、その他の変数（groupとage）は固定されていません。つまり個々のケースの値が代入されています。
同じ数値は
. margins, at(sex=(0 1))
でも得ることができます。
これに対してover()ですが、これは個々のケースの観察値ごとに予測値の平均値を計算するオプションです。
. margins, over(sex)
Predictive margins Number of obs = 3000
Model VCE : OLS
Expression : Linear prediction, predict()
over : sex
------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
sex |
0 | 64.5721 .5121636 126.08 0.000 63.56827 65.57592
1 | 74.88129 .5114812 146.40 0.000 73.87881 75.88378
------------------------------------------------------------------------------だと、個々のケースに個々の値を入れて計算した予測値を、sexについて平均したものです。したがって
. predict py
. tabstat py, by(sex)
のようにしても同じ結果を得ることができます。
over()とat( )は全く異なった働きをするのですが、その違いは少し理解しにくいと思います。簡単に言うと、at( )は個々のケースの属性に関係なく（強制的に）固定値を代入・指定するためのオプションです。たとえば男性のケースでも女性として扱って予測値を計算します。これに対してover()は、予測値の平均値を計算する際にケースの属性に応じて標本を分割して計算するためのオプションです。たとえば男性を男性として扱っていても女性として扱っていても、それとは独立に（観察された）性別ごとに予測値の平均値が出力されます。
「at( )は代入、over( )は標本分割」と覚えておくとよいでしょう。
counterfactualな予測値
では直接引数あるいはat( )と、over()を組み合わせるとどうなるでしょうか。たとえば
. margins, at(sex=(0 1)) over(sex)
Predictive margins Number of obs = 3000
Model VCE : OLS
Expression : Linear prediction, predict()
over : sex
1._at : 0.sex
sex = 0
1.sex
sex = 0
2._at : 0.sex
sex = 1
1.sex
sex = 1
------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_at#sex |
1 0 | 64.5721 .5121636 126.08 0.000 63.56827 65.57592
1 1 | 56.44061 .7184466 78.56 0.000 55.03248 57.84874
2 0 | 83.01278 .7179603 115.62 0.000 81.60561 84.41996
2 1 | 74.88129 .5114812 146.40 0.000 73.87881 75.88378
------------------------------------------------------------------------------とした場合です。ここで出力中_atとあるのはat( )オプションで固定した個々の数値です。1._atおよび2._atとありますが、それぞれsex=0およびsex=1を表しています。sexの値は0と1ですが、出力中では順番に1,2,3...と表示されるので注意。
出力された個々の値を解釈してみましょう。まずは一番上(64.57)は、_atが1でsexが0、つまり女性のケースについて、sex変数に女性(=0)を代入して個々のケースの予測値を計算し、その平均をとったものです。同様に一番下(74.88)は、_atが2でsexが1、つまり男性のケースについて、sex変数に男性(=1)を代入して同様に計算したものになります。これは要するにatの引数をとらないでそのまま予測値を計算し、それを性別ごとに平均した値と同じで、margins, over(sex)によって出力される値と同じです。
これに対して二番目(56.44)と三番目(83.01)は、（モデルが正しく特定化されている場合には）counterfactualな予測値になります。二番目は「男性のケースについて、sex変数に女性を、その他の変数(group, age)はケースの個々の値を代入した値」、三番目はその逆です。二番目についてもう少し説明すると、データ中の男性について、性別だけを女性にして、その他の条件を男性のままにした場合の値、ということになります。
over()とat( )を組み合わせることで、このように「男性のケース（=男性の共変量分布を持つケース）を女性として扱ったときの予測値」を知ることができるわけです。モデルが正しく特定化されている必要はありますが、以下のようなパターンの数値を得ることができます（順番は上記の出力に合わせています）。
- 女性の共変量を持つケース（要するに女性）に女性の効果を適用したときの予測値
- 男性の共変量を持つケース（要するに男性）に女性の効果を適用したときの予測値
- 女性の共変量を持つケース（要するに女性）に男性の効果を適用したときの予測値
- 男性の共変量を持つケース（要するに男性）に男性の効果を適用したときの予測値
モデルの特定化が不完全でも、調査データの場合は説明変数にもある程度の代表性があることが多いので、これらの「反事実的」な予測値は場合によっては参考になると思います。たとえば被説明変数に賃金、説明変数が性別、共変量として学歴・年齢・勤続年数をとった式を推定したとします。このとき上記のような予測値を計算すれば、「学歴・年齢・勤続年数が男性のまま、性別だけを女性にしたときの予測値」（およびその逆）を計算できるわけです。
atの異なる値の予測の際は、上記の場合にはもちろん回帰係数と一致します(83.01-64.57=18.44)。したがってここでみるべきなのは、むしろ1番目と2番目、あるいは3番目と4番目の差です（通常の回帰モデルでは両者の差は同じです）。これは、性別変数の効果に回収されない間接的な性別の効果とも言えますし、場合によっては「女性が置かれた環境にもし男性が置かれたらどうなるか」という反事実的想定のもとでの予測値であるとも言えます。
at( )のその他の指定法
at( )はかなり柔軟な固定値の指定を許容していますが、一定のルールがあります。まず、
margins sex, atmeans
とすれば、sex以外の共変量を平均値に固定した上でsexごとの予測値（2つ）が計算されます。この場合はすべての共変量が固定されていますので、得られるのは説明変数の分布に依存しない予測値です。同じ結果は
. margins sex, at( (means) _all)
でも得ることができますが、
. margins sex, at( (means) _all) at(age=(30(10)50))
とすると、まずat( (means) _all)で予測値（2個）が計算され、次にat(age=(30(10)50))で予測値（6個）が計算されます。つまり2つ別々の作業結果が出力されます（合計8個）。その他の変数は観察値のままで、ケースの予測値の平均値が出力されます。これに対して
. margins sex, atmeans at(age=(30(10)50))
だと、性別と年齢（30歳、40歳、50歳）の組み合わせ（6個）の予測値が、その他変数（group）を平均値に固定した上で計算されます。しかし
. margins sex, at( (means) _all age=(30(10)50))
とすると、「atmeans at(age=(30(10)50))」としたとき（直前）と同じ結果になります。
また、at((means) _all)の部分はたとえばat((median) _all)のように、様々な統計量を用いることができます。
このあたりの文法は煩雑なので、慣れるしかないかもしれません。ひとつだけ特に気にすべきだとすれば、出力されたものが「予測値の平均値」なのかどうか、という点です。直接引数やat( )オプションで説明変数を網羅して固定しているときは、予測値は説明変数の観察値の分布には依存しません。すべて固定されているからです。これに対して、固定していない変数が残っているときは、その変数の観察値の分布に応じた予測値の平均値が計算されます。（勝手に「固定値が指定されていない変数については固定値を平均値とみなす」というルールで計算されているわけではありません。）
共変量の値を固定して予測値を計算するのではなく、共変量の観察分布をもとに予測値を平均する方法は、後で述べるような限界値の計算で必要になる手続です。
交互作用モデルの予測値
以下のようにします（出力は省略）。また、直後にmarginsplotをして図を書いています。
. reg y i.sex#i.group age
. margins sex#group, atmeans
. marginsplot
カテゴリカル変数で交差項を作ると、セル度数が小さい場合には信頼区間が広がります。この場合、女性(sex=1）のgroup=3は該当人数が小さいので信頼区間が広くなっています。交差項を投入してもモデルが有意に改善されないときは（F検定や尤度比検定をする）、投入しないほうがよいでしょう。
また、以上のように、marginsコマンドの直後にmarginsplotをコマンドすることで、marginsで出力された内容を簡単に図にすることができます。これがv12の大きな変更点になっています。
センサーデータの予測値
tobit推定のあとは、センサーされていない予測値とセンサーされた予測値の2つを出力できます。たとえば
. tobit ycn i.sex#i.group age, ul(90)
のように、ycnが上限90でセンサーされている場合、
. margins sex#group
でセンサーされていない予測値を、
. margins sex#group, predict(ystar(.,90))
でセンサーされた予測値を推定できます。センサーされていない予測値は、上限の90を超えることがあります。図は、左がセンサー予測値、右がセンサーされていない予測値です。
限界効果の計算
限界効果とはある変数が一単位変化したときの被説明変数の変化ですが、ロジット、トービット、プロビットなどのモデルでは個々の説明変数の係数をそのまま限界効果として見ることができません。古典的回帰モデルでも、交互作用モデル、二乗項を含んだモデルで同様です。
. logit outcome i.treatment#(i.group c.age) c.age#c.age
（出力省略）
. margins, dydx(*)
Average marginal effects Number of obs = 3000
Model VCE : OIM
Expression : Pr(outcome), predict()
dy/dx w.r.t. : 1.treatment 2.group 3.group age
------------------------------------------------------------------------------
| Delta-method
| dy/dx Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
1.treatment | .0385625 .0162848 2.37 0.018 .0066449 .0704801
|
group |
2 | -.0776906 .0181584 -4.28 0.000 -.1132805 -.0421007
3 | -.1505652 .0400882 -3.76 0.000 -.2291366 -.0719937
|
age | .0095868 .0007796 12.30 0.000 .0080589 .0111148
------------------------------------------------------------------------------
Note: dy/dx for factor levels is the discrete change from the base level.上記のように、アステリスクで投入されたすべての説明変数の限界効果を出力します（個別の効果を知りたい時はその変数を括弧のなかに入れます）。モデルはかなり複雑ですが、限界値にするとすっきり解釈できます。この場合引数やatオプションがないので、共変量は観察値に設定され、各ケースの限界値の平均値が出力されます。
処置(treatment)の限界効果を、その他の共変量の値を固定してみたい時は、次のように入力します。
. margins, dydx(treatment) at(group=(1 2) age=(20(10)60))
. marginsplot, x(age)
dydxのところをeyexのようにすると弾性値が出力されます。
限界効果と回帰係数が一致する場合には、共変量を固定した上で関心のある変数の効果を単純に見るやり方が得策で、他の共変量を固定しないで予測値の平均値を取る意味はあまりないかもしれませんが、共変量に応じて関心のある変数の限界効果が変化するときは、予測値の平均値を見たほうがよいこともあるでしょう（固定値の選択の妥当性の問題）。ケースバイケースで使い分けする必要があります。
差の検定
いわゆるlinear combination（Stataではlincomでも検定可能）、あるいはcontrast statisticsの出力です。複数の予測値の95%信頼区間のみをみて、予測値間の差の検定とすることができない場合もあります。marginsコマンドでは、差の検定のときにはリファレンス値を指定して以下のように命令します。
. margins r.group, atmeans at(age=(1 2))
Contrasts of adjusted predictions
Model VCE : OIM
Expression : Pr(outcome), predict()
1._at : 0.sex = .4993333 (mean)
1.sex = .5006667 (mean)
group = 1
age = 39.799 (mean)
2._at : 0.sex = .4993333 (mean)
1.sex = .5006667 (mean)
group = 2
age = 39.799 (mean)
------------------------------------------------
| df chi2 P>chi2
-------------+----------------------------------
group@_at |
(2 vs 1) 1 | 1 17.52 0.0000
(2 vs 1) 2 | 1 17.52 0.0000
(3 vs 1) 1 | 1 35.24 0.0000
(3 vs 1) 2 | 1 35.24 0.0000
Joint | 2 35.46 0.0000
------------------------------------------------
--------------------------------------------------------------
| Delta-method
| Contrast Std. Err. [95% Conf. Interval]
-------------+------------------------------------------------
group@_at |
(2 vs 1) 1 | -.0615973 .0147146 -.0904374 -.0327573
(2 vs 1) 2 | -.0615973 .0147146 -.0904374 -.0327573
(3 vs 1) 1 | -.1086756 .018307 -.1445566 -.0727946
(3 vs 1) 2 | -.1086756 .018307 -.1445566 -.0727946
--------------------------------------------------------------
. marginsplot
以上は、男女それぞれについて、group=1とgroup=2および3の差の検定をしたものです。
とりあえず以上です。
- 作者: 筒井淳也,水落正明,秋吉美都,坂本和靖,平井裕久,福田亘孝
- 出版社/メーカー: ミネルヴァ書房
- 発売日: 2011/05/15
- メディア: 単行本
- 購入: 2人 クリック: 4回
- この商品を含むブログ (8件) を見る
