replication_Ross2012_table5_4

表5.4　内戦の勃発 1960-2006の再現

準備

この表では、ロジスティック回帰分析が行われている。Rでロジスティック回帰分析を行うときには、一般にはglm関数が使用される。しかし、glmではロスが使用したSTATAの頑健性標準誤差を得られない。この問題に対応するため、ここではfixestパッケージを使用する。必要に応じて、fixestパッケージをインストールしておく。

install.packages("fixest")#パッケージをインストールしていない場合。

library(fixest)

データの読み込み

load("Replication data for The Oil Curse - Ross 2012.RData")

1列目の分析

result1<-feglm(dom_conflict_onset_2_dich ~  logGDP_cap2000_sup_1 + pop_log_1 + peaceyears + X_spline1 + X_spline2 + X_spline3, data = x, family = "logit", cluster = "id")
summary(result1)

GLM estimation, family = binomial(link = "logit"), Dep. Var.: dom_conflict_onset_2_dich
Observations: 6,426 
Standard-errors: Clustered (id) 
                      Estimate Std. Error   t value   Pr(>|t|)    
(Intercept)          -6.090606   1.204936 -5.054716 4.3103e-07 ***
logGDP_cap2000_sup_1 -0.316016   0.060997 -5.180813 2.2092e-07 ***
pop_log_1             0.314434   0.072559  4.333502 1.4676e-05 ***
peaceyears            0.027245   0.088154  0.309059 7.5728e-01    
X_spline1             0.001718   0.001369  1.255137 2.0943e-01    
X_spline2            -0.001418   0.000842 -1.683892 9.2203e-02 .  
X_spline3             0.000630   0.000282  2.232526 2.5580e-02 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Log-Likelihood:  -801.7   Adj. Pseudo R2: 0.071675
           BIC: 1,664.8     Squared Cor.: 0.024154

上記の分析結果がglm関数を用いた場合とどのように異なるのか、確認しておこう。

result1_b<-glm(dom_conflict_onset_2_dich ~ logGDP_cap2000_sup_1 + pop_log_1 + peaceyears + X_spline1 + X_spline2 + X_spline3,
              data = x, family = binomial(link = "logit"))
summary(result1_b)

Call:
glm(formula = dom_conflict_onset_2_dich ~ logGDP_cap2000_sup_1 + 
    pop_log_1 + peaceyears + X_spline1 + X_spline2 + X_spline3, 
    family = binomial(link = "logit"), data = x)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.7930  -0.2762  -0.2094  -0.1453   3.1496  

Coefficients:
                       Estimate Std. Error z value Pr(>|z|)    
(Intercept)          -6.0906064  0.8766580  -6.948 3.72e-12 ***
logGDP_cap2000_sup_1 -0.3160159  0.0552011  -5.725 1.04e-08 ***
pop_log_1             0.3144342  0.0436499   7.204 5.87e-13 ***
peaceyears            0.0272446  0.0645479   0.422   0.6730    
X_spline1             0.0017180  0.0010659   1.612   0.1070    
X_spline2            -0.0014183  0.0006965  -2.036   0.0417 *  
X_spline3             0.0006301  0.0002641   2.386   0.0170 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1740.2  on 6425  degrees of freedom
Residual deviance: 1603.5  on 6419  degrees of freedom
  (2097 observations deleted due to missingness)
AIC: 1617.5

Number of Fisher Scoring iterations: 7

推定された回帰係数はfeglmとglmの間で変わらないが、標準誤差が異なる。feglmで推定した値の方が、Ross(2012)の頑健性標準誤差に近いことがわかる。以下では、feglmの分析結果のみを記載する。

2列目の分析

result2<-feglm(dom_conflict_onset_2_dich ~ logGDP_cap2000_sup_1 + pop_log_1 + logoil_gas_valuePOP_1 + peaceyears + X_spline1 + X_spline2 + X_spline3, data = x, family = "logit",cluster = "id")
summary(result2)

GLM estimation, family = binomial(link = "logit"), Dep. Var.: dom_conflict_onset_2_dich
Observations: 6,426 
Standard-errors: Clustered (id) 
                       Estimate Std. Error   t value   Pr(>|t|)    
(Intercept)           -4.632924   1.321424 -3.506009 4.5488e-04 ***
logGDP_cap2000_sup_1  -0.443890   0.069068 -6.426878 1.3025e-10 ***
pop_log_1              0.257728   0.077615  3.320577 8.9832e-04 ***
logoil_gas_valuePOP_1  0.133092   0.038309  3.474124 5.1252e-04 ***
peaceyears             0.039352   0.089671  0.438853 6.6077e-01    
X_spline1              0.001886   0.001388  1.358449 1.7432e-01    
X_spline2             -0.001501   0.000851 -1.763862 7.7755e-02 .  
X_spline3              0.000640   0.000283  2.262821 2.3647e-02 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Log-Likelihood:  -794.7   Adj. Pseudo R2: 0.078656
           BIC: 1,659.4     Squared Cor.: 0.026772

3列目の分析

3列目では、サンプルを中所得国（2000年基準価格の一人当たり国民所得が5000ドル未満）に限定する。「2000年基準価格の一人当たり国民所得」はgdpcap2000_supなので、これを使って部分集合を作成する。

x2<-x[x$gdpcap2000_sup<5000,]#部分集合の作成
result3<-feglm(dom_conflict_onset_2_dich ~ logGDP_cap2000_sup_1 + pop_log_1 + logoil_gas_valuePOP_1 + peaceyears + X_spline1 + X_spline2 + X_spline3, data = x2, family = "logit",cluster = "id")
summary(result3)

GLM estimation, family = binomial(link = "logit"), Dep. Var.: dom_conflict_onset_2_dich
Observations: 4,554 
Standard-errors: Clustered (id) 
                       Estimate Std. Error   t value   Pr(>|t|)    
(Intercept)           -5.545599   1.614431 -3.435018 0.00059252 ***
logGDP_cap2000_sup_1  -0.279814   0.097961 -2.856398 0.00428478 ** 
pop_log_1              0.255387   0.087886  2.905904 0.00366194 ** 
logoil_gas_valuePOP_1  0.123603   0.042515  2.907249 0.00364623 ** 
peaceyears             0.032184   0.094380  0.341005 0.73309965    
X_spline1              0.001708   0.001453  1.175587 0.23975996    
X_spline2             -0.001360   0.000900 -1.512101 0.13050815    
X_spline3              0.000574   0.000308  1.864977 0.06218461 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Log-Likelihood:  -709.2   Adj. Pseudo R2: 0.04709 
           BIC: 1,485.7     Squared Cor.: 0.023554

４列目の分析

４列目の分析は、サンプルを高所得国（2000年基準価格の一人当たり国民所得が5000ドル以上）に限定する。以下の分析結果はロスの分析結果と異なる。この違いは、ロスが作成したSTATAのプログラムが間違っているためである。ロスの間違いについては「ロスの分析に含まれる誤り」を参照されたい。STATA上でロスのプログラムの間違いを修正すると、以下のRの分析結果とほぼ一致する。

x3<-x[x$gdpcap2000_sup>=5000,]
result4<-feglm(dom_conflict_onset_2_dich ~ logGDP_cap2000_sup_1 + pop_log_1 + logoil_gas_valuePOP_1 + peaceyears + X_spline1 + X_spline2 + X_spline3, data = x3, family = "logit",cluster = "id")
summary(result4)

GLM estimation, family = binomial(link = "logit"), Dep. Var.: dom_conflict_onset_2_dich
Observations: 1,659 
Standard-errors: Clustered (id) 
                       Estimate Std. Error   t value   Pr(>|t|)    
(Intercept)           -6.940543   3.816240 -1.818686 6.8959e-02 .  
logGDP_cap2000_sup_1  -1.248528   0.357635 -3.491066 4.8110e-04 ***
pop_log_1              0.808567   0.123215  6.562242 5.3005e-11 ***
logoil_gas_valuePOP_1  0.166711   0.116884  1.426295 1.5378e-01    
peaceyears            -0.069441   0.240951 -0.288194 7.7320e-01    
X_spline1              0.000331   0.003007  0.109983 9.1242e-01    
X_spline2             -0.000597   0.001728 -0.345797 7.2950e-01    
X_spline3              0.000367   0.000623  0.588914 5.5592e-01    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Log-Likelihood: -60.0   Adj. Pseudo R2: 0.057093
           BIC: 179.4     Squared Cor.: 0.012176

正しいプログラムの分析結果はロスの推定結果と異なるが、ロスの議論そのものを損なうものではない。

５列目の分析

result5<-feglm(dom_conflict_onset_2_dich ~ logGDP_cap2000_sup_1 + pop_log_1 + logoil_gdpcap2000_sup_Q_1 + peaceyears + X_spline1 + X_spline2 + X_spline3, data = x, family = "logit", cluster = "id")
summary(result5)

GLM estimation, family = binomial(link = "logit"), Dep. Var.: dom_conflict_onset_2_dich
Observations: 6,382 
Standard-errors: Clustered (id) 
                           Estimate Std. Error   t value   Pr(>|t|)    
(Intercept)               -4.664662   1.320452 -3.532625 4.1146e-04 ***
logGDP_cap2000_sup_1      -0.410497   0.065377 -6.278954 3.4086e-10 ***
pop_log_1                  0.247017   0.077826  3.173974 1.5037e-03 ** 
logoil_gdpcap2000_sup_Q_1  0.108154   0.031597  3.422885 6.1960e-04 ***
peaceyears                 0.034017   0.087958  0.386737 6.9895e-01    
X_spline1                  0.001809   0.001369  1.320790 1.8657e-01    
X_spline2                 -0.001458   0.000843 -1.729764 8.3672e-02 .  
X_spline3                  0.000630   0.000283  2.229863 2.5757e-02 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Log-Likelihood:  -791.4   Adj. Pseudo R2: 0.077283
           BIC: 1,652.9     Squared Cor.: 0.026479

分析結果の一覧表示

feglmによって得られる分析結果のフォーマットはstargazerに対応していないので、ここではtexregパッケージのhtmlregを使って表示する。必要に応じてtexregパッケージをインストールする。

install.packages("texreg")#インストールしていない場合

library(texreg)
htmlreg(list(result1,result2,result3,result4,result5),
          custom.coef.map = list("logGDP_cap2000_sup_1" = "国民所得（対数）",
                                 "pop_log_1" = "人口（対数）",
                                 "logoil_gas_valuePOP_1" = "石油収入（対数）",
                                 "logoil_gdpcap2000_sup_Q_1" =  "国民所得*石油収入（対数）"), digits = 3)

Statistical models

	Model 1	Model 2	Model 3	Model 4	Model 5
国民所得（対数）	-0.316***	-0.444***	-0.280**	-1.249***	-0.410***
	(0.061)	(0.069)	(0.098)	(0.358)	(0.065)
人口（対数）	0.314***	0.258***	0.255**	0.809***	0.247**
	(0.073)	(0.078)	(0.088)	(0.123)	(0.078)
石油収入（対数）		0.133***	0.124**	0.167
		(0.038)	(0.043)	(0.117)
国民所得*石油収入（対数）					0.108***
					(0.032)
Num. obs.	6426	6426	4554	1659	6382
Deviance	1603.452	1589.305	1418.314	120.091	1582.801
Log Likelihood	-801.726	-794.652	-709.157	-60.046	-791.400
Pseudo R2	0.072	0.079	0.047	0.056	0.077
***p < 0.001; **p < 0.01; *p < 0.05

なお、ここで使用しているfixestパッケージのetable関数が一覧表示機能を有しているので、コンソールに表示するだけであればこれを使用するとよい（ここではブラウザに表示させているので見づらくなっている）。

etable(result1, result2, result3, result4, result5, cluster = "id",keep = c("logGDP_cap2000_sup_1", "pop_log_1", "logoil_gas_valuePOP_1", "logoil_gdpcap2000_sup_Q_1"))

                                            result1                   result2
Dependent Var.:           dom_conflict_onset_2_dich dom_conflict_onset_2_dich
                                                                             
logGDP_cap2000_sup_1            -0.3160*** (0.0610)       -0.4439*** (0.0691)
pop_log_1                        0.3144*** (0.0726)        0.2577*** (0.0776)
logoil_gas_valuePOP_1                                      0.1331*** (0.0383)
logoil_gdpcap2000_sup_Q_1                                                    
_________________________ _________________________ _________________________
S.E.: Clustered                              by: id                    by: id
Observations                                  6,426                     6,426
Squared Cor.                                0.02415                   0.02677
Pseudo R2                                   0.07857                   0.08670
BIC                                         1,664.8                   1,659.4

                                            result3                   result4
Dependent Var.:           dom_conflict_onset_2_dich dom_conflict_onset_2_dich
                                                                             
logGDP_cap2000_sup_1             -0.2798** (0.0980)        -1.249*** (0.3576)
pop_log_1                         0.2554** (0.0879)        0.8086*** (0.1232)
logoil_gas_valuePOP_1             0.1236** (0.0425)           0.1667 (0.1169)
logoil_gdpcap2000_sup_Q_1                                                    
_________________________ _________________________ _________________________
S.E.: Clustered                              by: id                    by: id
Observations                                  4,554                     1,659
Squared Cor.                                0.02355                   0.01218
Pseudo R2                                   0.05640                   0.15554
BIC                                         1,485.7                    179.40

                                            result5
Dependent Var.:           dom_conflict_onset_2_dich
                                                   
logGDP_cap2000_sup_1            -0.4105*** (0.0654)
pop_log_1                         0.2470** (0.0778)
logoil_gas_valuePOP_1                              
logoil_gdpcap2000_sup_Q_1        0.1082*** (0.0316)
_________________________ _________________________
S.E.: Clustered                              by: id
Observations                                  6,382
Squared Cor.                                0.02648
Pseudo R2                                   0.08537
BIC                                         1,652.9
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1