『石油の呪い』図6.9
Masaki Matsuo
はじめに
『石油の呪い』には有益で魅力的な図が豊富に収録されている。しかし、その全てにおいてデータと作成方法は公表されていない。もちろん、全ての図にはデータの出所が示されているため、自分でデータを収集してなるべく近いものを作成することは可能である。しかし、自分でデータを収集するのはあまりにも手間がかかる。
しかし幸運なことに、『石油の呪い』の図6.9に関しては、公表されているデータセットをほんのわずか修正するだけで作成することが可能である。また、図6.9は(1)散布図を作成し、(2)観察対象の分類に従って凡例を変え（産油国が黒丸、非産油国が白抜き丸）、(3)産油国にラベルをつけ（アルファベット３文字の国コードをつける）、(4)回帰直線を描くという、計量分析の基本的な技術を習得するのに適している。
以上の点を踏まえ、ここでは図6.9を作成する手順を通して、散布図と回帰直線を使って、分析をビジュアルで表現する方法を学ぶ。
散布図作成を作成する
まず、データを読み込む。作業ディレクトリは各自で設定すること。
load("Replication data for The Oil Curse - Ross 2012.RData")図6.9は2005年のデータを元に作成されている（日本語訳には示されていないが、原文では2005年のデータであることが示されている）。このため、2005年のデータだけを抽出して部分集合を作成する。部分集合の作成にはsubset()を使用する。
y2005<-subset(x,year==2005)縦軸は「政府の効率性」であり、これはロスのデータセットの中ではgoveffectに該当する。横軸の「一人あたり国民所得（対数）」は同じくlogGDPcapである。この二つの変数を使うだけで散布図を描くことはできる。まずは難しいことを考えず、散布図を作成してみよう。使用する関数はplot()である。
plot()の書式には、x軸とy軸の値を引数で指定する方法と、モデル式で表す方法がある。ここでは、モデル式で表す方法を採用している。
plot(goveffect~logGDPcap,data = y2005)
グループ分け
図6.9では、人口一人当たり1000ドル以上の石油収入を有する国にラベルをつけ、黒丸で描画している。しかし、ロスのデータセットには1000ドル以上を抜き出すダミー変数はない。また、ロスのデータセットに収録されている一人当たりの石油収入は対数である。このため、1000の対数を計算し、この値を元に条件判断を行い、1000の対数を超える対象のみ点を黒く表示し、それ以外を白抜きで表示する。点の表示はpch()で指定する。引数に16を指定すると黒丸、21を指定すると白抜きの丸になる。
条件判断にはifelse()を用いる。書式はifelse(判断基準,Yes,No)の形式である。判断基準に合致する場合、Yesの値が与えられ、合致しない場合はNoの値が与えられる。判断基準は「人口一人当たり1000ドル以上の石油収入を有するかどうか」である。ロスのデータセットにあるlogoil_gas_valuePOPが人口一人当たりの石油・ガス収入の対数値なので、これを元に1000の対数を超える場合がYes、超えない場合がNoになる。超えた場合に黒丸で示すようにしたいので、ifelse()の条件判断の結果をpch()の引数指定と組み合わせて用いると良い。
plot(goveffect~logGDPcap,data = y2005,pch=ifelse(y2005$logoil_gas_valuePOP>=log(1000),16,21))
黒丸と白抜き丸で産油国と非産油国を分けて表示することができた。
産油国にラベルをつける
ロスのデータセットには、アルファベット３文字の国コードはidに格納されている。各点にラベルをつけるときには、text()で座標を指定し、引数のlabelsでidを指定することで、アルファベット３文字のコードをつけることができる。しかし、単にこれを実行すると、白抜き丸（非産油国）にもラベルがつけられてしまう。
plot(goveffect~logGDPcap,data = y2005,
pch=ifelse(y2005$logoil_gas_valuePOP>=log(1000),16,21))
text(x=y2005$logGDPcap,y=y2005$goveffect+0.1,labels = y2005$id,cex=0.5)
#点の上にラベルを表示させるため、yの値に0.1を足した値をy座標に指定している。
#ラベルの文字の大きさをcexで小さくしている。ラベルを産油国に限定するにはどうしたらいいだろうか。ここでは、産油国にのみアルファベット３文字の国コードが格納され、それ以外にはデータが格納されていない変数id3を新たに作成する(id2は既にロスのデータセットで使用されている)。これをlabelsの引数に指定することで、非産油国には何もないデータが表示され（つまり何も表示されず）、産油国にのみ国コードが表示されるようなプログラムを作成することで、この問題を解決しよう。産油国と非産油国の識別方法は、上記のifelse()と同じ方法を採用する。これがYesの場合、元々のロスのデータセットの国コード(つまり変数id)が転記され、Noの場合は何も記載されないというプログラムを組み、これをy2005に新たに追加するid3に格納する。
y2005$id3<-ifelse(y2005$logoil_gas_valuePOP>=log(1000),y2005$id,"")では、こうして出来上がったid3を用いて、散布図にラベルを付けよう。
plot(goveffect~logGDPcap,data = y2005,
pch=ifelse(y2005$logoil_gas_valuePOP>=log(1000),16,21))
text(x=y2005$logGDPcap,y=y2005$goveffect+0.1,labels = y2005$id3,cex=0.5)
これで産油国のみにラベルを貼ることができた。
回帰直線の記入
次に、散布図の上から回帰直線を描画しよう。最も簡単な方法は、直線を描く関数であるabline()の引数に回帰分析の結果を指定する方法である。ここでは、回帰分析は従属変数がgoveffect、独立変数がlogGDPcapである。分析結果をresultに格納することにしよう。
result<-lm(goveffect~logGDPcap,data = y2005)この分析結果を上記の散布図の上から描くと、以下のようになる。
plot(goveffect~logGDPcap,data = y2005,
pch=ifelse(y2005$logoil_gas_valuePOP>=log(1000),16,21))
text(x=y2005$logGDPcap,y=y2005$goveffect+0.1,labels = y2005$id3,cex=0.5)
abline(result,lty="dashed")#ltyで破線を指定した。
それっぽくなってきた。しかし、再現したい図6.9では、回帰直線はデータが存在しない部分は描かれていないのに対して、ここで描いた回帰直線はグラフの端から端まで描かれている。
回帰直線の書き方にはいくつかの考え方がある。その一つに、現実を分析した結果を示す図に、現実には存在しないトレンドを描くことは不適切だ、という考え方がある。図6.9はそのような考え方に基づいて描かれたものだろう。この考えに基づいて、ここでも現実のデータが存在しない領域には、回帰直線を描かないことにしたい。これは、どうすればいいだろうか。
予測結果を描く
lines()を使うと、線の両端の座標を指定して線を描くことができる。この座標指定において、現実のデータ範囲を超えないように指定すれば、図6.9の回帰直線を再現することができる。具体的には、以下の手順でおこなう。
線の左端の座標を(x1,y1)、右端を(x2,y2)と置いた場合、x1はx軸データであるlogGDPcapの最小値に該当し、同じくx2は最大値に該当する。y1とy2は、我々が推定した回帰式(つまりresult)にx1とx2を代入すれば求めることができる。推定した回帰式に値を代入して結果を得る行為を「予測」と呼ぶ。Rではpredict()を用いてこれを行う。
predict()は、引数に回帰式と代入するデータを指定して用いる。回帰式はresultであり、代入するデータはx1とx2があれば良い。x1とx2を数値で指定するのは煩雑なので、最小値を得るmin()と、最大値を得るmax()を利用し、その結果をpdl.rangeに格納しておこう。その際、データにNAが含まれていると最小値と最大値が得られないので、NAがある場合にはそれを参照しないようにオプションを設定しておこう。これは具体的には、\(na.rm=TRUE\)である。
pdl.range<-c(min(y2005$logGDPcap,na.rm = T),max(y2005$logGDPcap,na.rm = T))
pdl.range[1] 4.510586 10.862510これでx1とx2の値を得ることができた。これをpredict()の引数に使用し、y1とy2を得る。
pdl<-predict(result,newdata =data.frame(logGDPcap=pdl.range))
pdl 1 2
-1.67931 1.67614 これらの値をlines()の引数に設定して、散布図の上に回帰直線を描く。
plot(goveffect~logGDPcap,data = y2005,
pch=ifelse(y2005$logoil_gas_valuePOP>=log(1000),16,21))
text(x=y2005$logGDPcap,y=y2005$goveffect+0.1,labels = y2005$id3,cex=0.5)
lines(pdl.range,pdl,lty="dashed")
軸ラベルをつける
あとは、x軸とy軸のラベルを設定すれば良い。軸ラベルの設定はplot()の引数にxlabとylabを用いて行う。Windowsの場合は問題ないが、Macで日本語を表示させるには、plot()を指定する前にpar()でフォントを指定する必要がある。なお、ここでは点の大きさをcexでやや小さめに設定した。
par(family="HiraKakuPro-W3")#Macの場合のみ設定。
plot(xlim=c(4,12),ylim=c(-2,2.5),goveffect~logGDPcap,data = y2005,
pch=ifelse(y2005$logoil_gas_valuePOP>=log(1000),16,21),cex=0.8,
xlab="一人当たり国民所得（対数)",ylab="政府の効率性")
text(x=y2005$logGDPcap,y=y2005$goveffect+0.1,labels = y2005$id3,cex=0.5)
lines(pdl.range,pdl,lty="dashed")