Stataによる予測値のプロット(1)

最近の論文では、回帰モデル推定後に予測値をプロットすることが多いようです。Stataで分析しているときにそれをやる方法についてまとめておきます。

まず、連続量の説明変数による被説明変数の予測値の変化を見たい場合、Stataで最も簡単に直接グラフを出力するやり方は、predictコマンドを使うことです。架空のデータを使って説明します。（変数colは大卒ダミー、age2は年齢二乗項。）

use http://homepage3.nifty.com/sociology/data/predict, clear
reg income male col age age2
predict pincome

これでpincomeという予測値が各ケースごとに出力されました。これを年齢をX軸にとってプロットしますが、何も考えずにすると以下のようになります。

scatter pincome age, scheme(economist)

これは4つのカテゴリー（性別×学歴）ごとに収入予測値がプロットされるからです。年齢との交互作用が推定されていない場合、何かにカテゴリーを決めてプロットしても良いのですが、age以外の変数について平均値をとったとき、という設定で出力するのが一般的でしょう。

オーソドックスな方法は、個々の推定された係数と説明変数の平均値を使って、マニュアル的に予測値変数をつくってやることです。

まずは共変量の平均値をマクロに保存する。

su male
global m_male=r(mean)
su col
global m_col=r(mean)

次にそれを使って、ageとage2以外については平均値を使って予測値を計算する。

qui reg income male col age age2
gen pincome2=_b[_cons]+_b[male]*$m_male+_b[col]*$m_col+
_b[age]*age+_b[age2]*age2
scatter pincome2 age, scheme(economist)

とはいえ、この方法では共変量が多くなったときはかなり面倒になるので、次のような方法（裏ワザ？）を使います。

qui reg income male col age age2
adjust male col, by(age age2) gen(pincome3)

adjustコマンドは引数にとった変数（ここではmaleとcol）については平均値であると措定し、オプションbyの引数にとった変数の個々の値の予測値を出力してくれるコマンドです。ともあれ、これでpincome2と全く同じ数値がpincome3として出力されます。（adjustコマンドはversion 11からmarginsという新しいコマンドに置き換えられています。ただadjust自体はまだ使えます。）

ロジットきょうだい(*logit)の場合、予測値はprオプションを使った確率予測値を使うのが一般的だと思います。ここでもやり方は同じですが、せっかくだから学歴だけ平均値で、性別と年齢ごとのプロットをしてみます。

logit seiki male col age age2
qui adjust col, by(male age age2) gen(pseiki) pr
scatter pseiki age, scheme(economist)

続きは後日。

補足

あとから気づきました。以下のようにすれば、スムーズな曲線になるのでした。

tw (function seiki=invlogit(_b[_cons]+_b[male]*0+_b[col]*1+_b[age]*x), range(30 60)) (function seiki=invlogit(_b[_cons]+_b[male]*1+_b[col]*1+_b[age]*x), range(30 60)), scheme(economist) legend(off)