lgraph: Stataで任意のエラーバー付きのグラフを描く

はじめに

分析をしているときに、いろいろな場面でモデルから得られる予測値と信頼区間をプロットして、視覚的に把握したいというときがあります。Stataであればmarginsmarginsplotというコマンドがあり、これを使うときれいな予測値のグラフが書けます。

しかし、Stata内で行った回帰分析から得られた結果であればいいのですが、場合によっては、予測値と標準誤差のデータだけしか手元にないとか、3重交互作用など複雑なパターンであるとか、様々な事情のためにmarginsの範囲内ではいい感じに図示できないことがあります。

こういうとき、Excelでは、信頼区間を手で入力してグラフを作成するということができます。たとえば以下の記事のような感じです。

→ Excelで棒グラフに個別のエラーバーをつける方法 

→ Excelでエラーバー付きグラフを作成する

これと似たようなことをStataでできないかと探していたところTimothy Makさんという方が作成したlgraphという関数を使うとできるようです。使うためにはこの関数をStataにインストールする必要があります。

findit lgraph

今回はこの関数の使い方メモです。

どう使う?

たとえば以下のように、車の重量、および車種(国内車か外車か)によって燃費の予測値がどの程度異なるか、およびその予測値の標準誤差に関するデータを持っているとします(このデータの作り方自体は後述します)。

weightというのが重さ、foreignというのが車種で、これらの変数の値ごとに、予測値(predicted)とその標準誤差(stderr)の値が並んでいる、という状態です。

こういうデータが手元にあるときに、予測値と信頼区間を簡単に図示できる関数がlgraphということになります。

まず、上側信頼区間と下側信頼区間の値を計算します。

gen upperbound = predicted + stderr * 1.96
gen lowerbound = predicted - stderr * 1.96

今回は95%信頼区間を図示しますが、1.96の部分の値を変えれば任意の信頼区間を求めることができます。

あとはlgraphを使って次のように記述するだけです。一応、y軸とx軸のタイトルもつけておきます。

lgraph predicted weight ///
, by(foreign) errortype(collapse(upperbound lowerbound)) ///
ytitle("Mileage") xtitle("Weight (lb)")

このような図ができます。

これはデフォルトの図で、これだけでも十分綺麗だと思いますが、いろいろオプションを指定することで見た目も変えられます。

lgraph predicted weight ///
, by(foreign) errortype(collapse(upperbound lowerbound)) ///
ytitle("Mileage") xtitle("Weight (lb)") ///
loptions(recast(scatter) ; 0 color(navy) ; 1 color(cranberry) ) ///
eoptions(lp(solid) ; 0 color(navy) ; 1 color(cranberry) ) ///
 scheme(s1mono) ylabel(0(10)30, grid) scale(1.2) title("Predicted values", box bexpand size(medium))

X軸が連続変数でない場合は、点をつながないグラフのほうがよいと思います。そのほか、くわしいオプションや使い方はlgraphのhelpを参照してください。

(参考)今回のデータの作り方

適当なデータを使って、予測値と標準誤差だけがあるという状況を作ります。ここではestoutパッケージ群に入っているesttabという関数を使います。入っていない人はあらかじめインストールしておいてください(このパッケージはめちゃめちゃ便利なのでぜひ入れておくのをおすすめします)。

ssc install estout

使用するデータはStataヘルプページを頻繁に見ている人にとってはおなじみ(?)のautoデータです。

use "http://www.stata-press.com/data/r13/auto", clear

regress mpg c.weight##i.foreign 
margins foreign, at(weight = (2000(500)4500)) post
est sto margins1

esttab margins1 using margins1.csv, b se nogaps wide nostar plain noobs replace

** ディレクトリの名前は架空です。
import delimited "/Users/RyotaMugiyama/Desktop/margins1.csv", clear 

browse

読み込んだデータは次のような状態です。v2列(b)というのが各カテゴリの予測値で、v3列(se)というのが標準誤差です。このデータをもう少しちゃんとした状態にします。

gen row = _n 
drop if row <= 2

rename v2 predicted
rename v3 stderr

gen weight = 1
replace weight = weight + weight[_n-2] if weight[_n-2] ~= .
replace weight = (weight - 1) * 500 + 2000
bysort weight: gen foreign = _n - 1

drop v1
destring, replace

lab def foreignlab 0"Domestic" 1"Foreign"
lab val foreign foreignlab

これで先ほどのデータができます。見るときはこんな感じ。

browse weight foreign predicted stderr

(参考)そのほかの作り方

このようなウェブページもありました。

How can I make a bar graph with error bars? | Stata FAQ

Say that you were looking at writing scores broken down by race and ses. You might want to graph the mean and confidence interval for each group using a bar chart with error bars as illustrated below. This FAQ shows how you can make a graph like this, building it up step by step.

個人的に棒グラフにエラーバーつけるのはあまり見た目が良くないように思うのでやらないのですが、Stataだとグラフの種類自体はオプションでrecast()などを使うとわりと簡単に変えられるので、このやり方も参考になると思います。

あとこのページも見つけました。ひょっとするとこちらのほうがシンプルで柔軟性が高いかもしれません。

割合と信頼区間のグラフを描く - researchmap

割合の95%信頼区間を算出する にて,割合と信頼区間の算出を紹介しました。 これらをStataでをグラフにするには, 割合と信頼区間それぞれが変数となった新たなデータセットが必要です。 以下,バ...

土屋政雄先生のStataメモのページはかゆいところに手が届く感がすごく、おすすめです。

(参考)今回のdofile