replication_Ross2012

『石油の呪い』第3章表3.7の再現

まず、ライブラリを読み込んでおく。

library(lme4)#ロジット・パネル分析で使用
library(stargazer)#分析結果の表の出力に使用

データの読み込み

load("Replication data for The Oil Curse - Ross 2012.RData")#データ読み込み

パネルデータ分析を行う際には、１列目をID、２列目を時間（この場合は年）にしておく。ロスのオリジナルのデータでは、６列目にid2があるのでこれをIDに利用し、7行目のyearを時間（年）に使用する。このため、最初の５列を削除し、x2に格納してこれを使用する。

x2<-x[,-c(1:5)]#データ整形

パネル・ロジット分析

ロスの解説にある通り、この分析では1960年から5年ごとので区切られた期間ダミーが投入される。これはロスのデータセットのYEAR1960_64 YEAR1965_69 YEAR1970_74 YEAR1975_79 YEAR1980_84 YEAR1985_89 YEAR1990_94 YEAR1995_99 YEAR2000_04 YEAR2005_09 に該当する。ただし、最後のYEAR2005_09を投入しても共線性を産むために自動で削除されるので、以下のモデルには含んでいない。

#モデル１（統制変数のみ）
result1<-glmer(dem_transition ~ logGDP_cap2000_sup_1 + log_authyrs_1 + YEAR1960_64 + YEAR1965_69 + YEAR1970_74 + 
                 YEAR1975_79+YEAR1980_84 + YEAR1985_89 + YEAR1990_94 + YEAR1995_99 +YEAR2000_04+(1|id2),data=x2,
               family=binomial,control =glmerControl(optimizer = "bobyqa"),nAGQ = 12)
summary(result1)

Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 12) [glmerMod]
 Family: binomial  ( logit )
Formula: dem_transition ~ logGDP_cap2000_sup_1 + log_authyrs_1 + YEAR1960_64 +  
    YEAR1965_69 + YEAR1970_74 + YEAR1975_79 + YEAR1980_84 + YEAR1985_89 +  
    YEAR1990_94 + YEAR1995_99 + YEAR2000_04 + (1 | id2)
   Data: x2
Control: glmerControl(optimizer = "bobyqa")

     AIC      BIC   logLik deviance df.resid 
   795.8    876.4   -384.9    769.8     3626 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-0.8869 -0.1446 -0.1061 -0.0693 10.1227 

Random effects:
 Groups Name        Variance Std.Dev.
 id2    (Intercept) 2.625    1.62    
Number of obs: 3639, groups:  id2, 125

Fixed effects:
                       Estimate Std. Error z value Pr(>|z|)   
(Intercept)          -4.3222185  1.3325368  -3.244  0.00118 **
logGDP_cap2000_sup_1  0.1392227  0.1505823   0.925  0.35519   
log_authyrs_1        -0.0004891  0.1791622  -0.003  0.99782   
YEAR1960_64          -1.8946762  0.7784805  -2.434  0.01494 * 
YEAR1965_69          -2.1162789  0.7032924  -3.009  0.00262 **
YEAR1970_74          -2.2346520  0.6998821  -3.193  0.00141 **
YEAR1975_79          -1.5107291  0.5851695  -2.582  0.00983 **
YEAR1980_84          -1.4422578  0.5821806  -2.477  0.01324 * 
YEAR1985_89          -0.7425178  0.5211552  -1.425  0.15423   
YEAR1990_94           0.2523970  0.4498043   0.561  0.57471   
YEAR1995_99          -0.9428532  0.5710373  -1.651  0.09871 . 
YEAR2000_04          -0.2964313  0.4959482  -0.598  0.55004   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) lGDP_2 lg_t_1 YEAR1960 YEAR1965 YEAR1970 YEAR1975 YEAR1980
lGDP_2000__ -0.851                                                           
log_thyrs_1 -0.511  0.100                                                    
YEAR1960_64 -0.166  0.056 -0.014                                             
YEAR1965_69 -0.155  0.037 -0.039  0.421                                      
YEAR1970_74 -0.093 -0.006 -0.105  0.405    0.455                             
YEAR1975_79 -0.118 -0.010 -0.106  0.462    0.514    0.517                    
YEAR1980_84 -0.092 -0.022 -0.150  0.460    0.507    0.511    0.600           
YEAR1985_89 -0.132  0.025 -0.195  0.449    0.499    0.504    0.595    0.630  
YEAR1990_94 -0.303  0.106 -0.020  0.436    0.481    0.480    0.573    0.581  
YEAR1995_99 -0.257  0.096  0.006  0.325    0.359    0.359    0.430    0.433  
YEAR2000_04 -0.276  0.082  0.027  0.350    0.389    0.386    0.463    0.466  
            YEAR1985 YEAR1990 YEAR1995
lGDP_2000__                           
log_thyrs_1                           
YEAR1960_64                           
YEAR1965_69                           
YEAR1970_74                           
YEAR1975_79                           
YEAR1980_84                           
YEAR1985_89                           
YEAR1990_94  0.645                    
YEAR1995_99  0.482    0.564           
YEAR2000_04  0.520    0.606    0.489

#モデル２（石油収入を追加）
result2<-glmer(dem_transition ~ logGDP_cap2000_sup_1 + log_authyrs_1 +logoil_gas_valuePOP_1+ YEAR1960_64 +
                 YEAR1965_69 +YEAR1970_74 + YEAR1975_79 + YEAR1980_84 + YEAR1985_89 + YEAR1990_94 + YEAR1995_99 +
                 YEAR2000_04+(1|id2),data=x2, family=binomial,control =glmerControl(optimizer = "bobyqa"),nAGQ = 12)

モデル１のlog_authyrs_1の推定値はロスの分析結果と異なるが、おそらくStataとRの間での計算処理の違いだと思われる。統計ソフトが実行するパネル・ロジット分析には複雑な計算が採用されており、様々なオプションでこれを調整するが、その数学的背景は筆者の能力を超えるので、ここでは取り扱わ（え）ない。ただし、この変数の推定値の符号が負であり、また統計的に有意でないことはロスの分析と一致しているために、大きな問題ではないと思われる。モデル２以降の分析結果はほぼロスの分析と一致するので、summary( )を使って結果を表示していない。

モデル３とモデル４では、サンプルの対象年を限定する。このため、ここではsubsetを使って対象年の異なる部分集合を二つ作成して使用する。なお、サンプルに期間の制限がかけられているため、それに合わせて投入する期間ダミーを調整した。

#モデル３（1960-1979年）
pre1980<-subset(x2,year<1980)
result3<-glmer(dem_transition ~ logGDP_cap2000_sup_1 + log_authyrs_1 +logoil_gas_valuePOP_1+ YEAR1960_64 + YEAR1965_69 +
                 YEAR1970_74 +(1|id2),data=pre1980, family=binomial,control = glmerControl(optimizer = "bobyqa"),nAGQ = 12)
#モデル４（1980-2006)
post1980<-subset(x2,year>1979)
result4<-glmer(dem_transition ~ logGDP_cap2000_sup_1 + log_authyrs_1 +logoil_gas_valuePOP_1+YEAR1980_84 + YEAR1985_89 +
                 YEAR1990_94 + YEAR1995_99 +YEAR2000_04+(1|id2),data=post1980, family=binomial,
               control = glmerControl(optimizer = "bobyqa"), nAGQ = 12)

この二つのモデルの分析結果は、ロスの分析結果とほぼ同じである。

#モデル５（石油収入と南米ダミーの交差項）
result5<-glmer(dem_transition ~ logGDP_cap2000_sup_1 + log_authyrs_1 +logoil_gas_valuePOP_1+ latin_logoil_1+YEAR1960_64 +
                 YEAR1965_69 +YEAR1970_74 + YEAR1975_79 + YEAR1980_84 + YEAR1985_89 + YEAR1990_94 + YEAR1995_99 +
                 YEAR2000_04+(1|id2),data=x2, family=binomial,control = glmerControl(optimizer = "bobyqa"), nAGQ = 12)
#モデル６（過去の民主制、経済成長、ムスリム人口）
result6<-glmer(dem_transition ~ logGDP_cap2000_sup_1 + log_authyrs_1 +logoil_gas_valuePOP_1+latin_logoil_1+
                 dem_prior_update+ gdpgrowth_alt_new_1+ mus_pct_sup_1+YEAR1960_64+YEAR1965_69 +YEAR1970_74 +
                 YEAR1975_79 +YEAR1980_84 + YEAR1985_89 + YEAR1990_94 +YEAR1995_99 +YEAR2000_04+(1|id2),data=x2,
               family=binomial, control = glmerControl(optimizer = "bobyqa"),nAGQ = 12)

Warning in optwrap(optimizer, devfun, start, rho$lower, control = control, :
convergence code 1 from bobyqa: bobyqa -- maximum number of function evaluations
exceeded

この二つのモデルも、ロスの分析結果とほぼ同じである。なお、モデル６には上記のwarningが表示されるが、ここでは無視する。

分析結果の表示

stargazerを用いて分析結果を表で記述する。ロスの表は年次や国家数、観察数、欠損値の割合を含むが、この中で観察数以外はRでは別途算出する必要があるので、ここでは表示していない。

stargazer(result1,result2,result3,result4,result5,result6,keep =c("logGDP_cap2000_sup_1","log_authyrs_1",
                                                                  "logoil_gas_valuePOP_1",
                                                                   "latin_logoil_1","dem_prior_update",
                                                                   "gdpgrowth_alt_new_1","mus_pct_sup_1"),
          keep.stat = "n",type = "html")


	Dependent variable:

	dem_transition
	(1)	(2)	(3)	(4)	(5)	(6)

logGDP_cap2000_sup_1	0.139	0.333^*	0.406	0.124	0.256	0.225
	(0.151)	(0.171)	(0.254)	(0.134)	(0.168)	(0.162)

log_authyrs_1	-0.0005	-0.045	-0.122	-0.438^***	0.010	0.387^*
	(0.179)	(0.182)	(0.303)	(0.141)	(0.182)	(0.229)

logoil_gas_valuePOP_1		-0.179^**	-0.078	-0.129^**	-0.292^***	-0.197^**
		(0.079)	(0.114)	(0.059)	(0.092)	(0.090)

latin_logoil_1					0.673^***	0.414^***
					(0.152)	(0.143)

dem_prior_update						1.915^***
						(0.465)

gdpgrowth_alt_new_1						-0.054^***
						(0.017)

mus_pct_sup_1						-0.721
						(0.551)


Observations	3,639	3,507	1,297	2,210	3,507	3,422

Note:	p<0.1; p<0.05; p<0.01

replication_Ross2012_table3.7

Masaki Matsuo

2/18/2021

『石油の呪い』第3章表3.7の再現

パネル・ロジット分析

分析結果の表示