require(ggplot2)
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.2.5
require(sandwich)
## Loading required package: sandwich
require(msm)
## Loading required package: msm
## Warning: package 'msm' was built under R version 3.2.5
require(foreign)
## Loading required package: foreign
## Warning: package 'foreign' was built under R version 3.2.5
require(MASS)
## Loading required package: MASS
require(gdata)
## Loading required package: gdata
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'gdata'
library(dplyr)
## Warning: package 'dplyr' was built under R version 3.2.5
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
##
## select
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
#getwd()
HD<-read.csv("HD.csv")
#head(HD)
#attributes(HD)
dim(HD)
## [1] 6200 17
#subsetting for our time frame 1981-2006
HD1<- subset(HD, YEAR >= 1981 & YEAR <2006)
#head(HD1)
dim(HD1)
## [1] 5000 17
##### Model 2 #########
summary(m2 <- glm.nb(incidentsusvictims_L1 ~ lnpopterr + logGDPCAP+logUSMilitaryAid+logArea+Polity2+CINC+WOPOL, data = HD1))
##
## Call:
## glm.nb(formula = incidentsusvictims_L1 ~ lnpopterr + logGDPCAP +
## logUSMilitaryAid + logArea + Polity2 + CINC + WOPOL, data = HD1,
## init.theta = 0.1200766362, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.8963 -0.5380 -0.4133 -0.2891 5.0351
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -13.40181 1.15842 -11.569 < 2e-16 ***
## lnpopterr 0.83001 0.07467 11.116 < 2e-16 ***
## logGDPCAP 0.07921 0.05383 1.472 0.141154
## logUSMilitaryAid 0.04206 0.01013 4.153 3.28e-05 ***
## logArea -0.11864 0.05245 -2.262 0.023709 *
## Polity2 0.03433 0.01231 2.789 0.005285 **
## CINC -15.68286 4.47304 -3.506 0.000455 ***
## WOPOL -0.68883 0.12077 -5.704 1.17e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.1201) family taken to be 1)
##
## Null deviance: 1266.6 on 3190 degrees of freedom
## Residual deviance: 1004.9 on 3183 degrees of freedom
## (1809 observations deleted due to missingness)
## AIC: 3175.2
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.12008
## Std. Err.: 0.00997
##
## 2 x log-likelihood: -3157.22700
cov.m2 <- vcovHC(m2, type="HC0")
std.err2 <- sqrt(diag(cov.m2))
r.est2 <- cbind(Estimate= coef(m2), "Robust SE" = std.err2,
"Pr(>|z|)" = 2 * pnorm(abs(coef(m2)/std.err2), lower.tail=FALSE),
LL = coef(m2) - 1.96 * std.err,
UL = coef(m2) + 1.96 * std.err)
round(r.est2,3)
## Estimate Robust SE Pr(>|z|) LL UL
## (Intercept) -13.402 1.520 0.000 -16.069 -10.734
## lnpopterr 0.830 0.072 0.000 0.697 0.963
## logGDPCAP 0.079 0.055 0.147 -0.043 0.202
## logUSMilitaryAid 0.042 0.016 0.009 0.010 0.075
## logArea -0.119 0.051 0.021 -0.217 -0.020
## Polity2 0.034 0.013 0.008 0.006 0.062
## CINC -15.683 4.421 0.000 -23.650 -7.716
## WOPOL -0.689 0.117 0.000 -0.977 -0.400
with(m2, cbind(res.deviance = deviance, df = df.residual,
p = pchisq(deviance, df.residual, lower.tail=FALSE)))
## res.deviance df p
## [1,] 1004.851 3183 1
#To present the regression results as incident rate ratios and
#their standard errors, together with the confidence interval.
#To compute the standard error for the incident rate ratios,
#use the Delta method. To this end, we make use the function deltamethod implemented in R package msm.
s2 <- deltamethod(list(~ exp(x1), ~ exp(x2), ~ exp(x3), ~ exp(x4),~ exp(x5),~ exp(x6),~ exp(x7),~ exp(x8)), coef(m2), cov.m2)
## exponentiate old estimates dropping the p values
rexp.est2 <- exp(r.est2[, -3])
## replace SEs with estimates for exponentiated coefficients
rexp.est2[, "Robust SE"] <- s2
rexp.est2
## Estimate Robust SE LL UL
## (Intercept) 1.512398e-06 2.299234e-06 1.050096e-07 2.178227e-05
## lnpopterr 2.293335e+00 1.647801e-01 2.007472e+00 2.619904e+00
## logGDPCAP 1.082431e+00 5.913806e-02 9.578109e-01 1.223266e+00
## logUSMilitaryAid 1.042956e+00 1.675618e-02 1.009559e+00 1.077458e+00
## logArea 8.881317e-01 4.554399e-02 8.048948e-01 9.799763e-01
## Polity2 1.034929e+00 1.339285e-02 1.006482e+00 1.064180e+00
## CINC 1.545326e-07 6.831190e-07 5.358287e-11 4.456708e-04
## WOPOL 5.021639e-01 5.858028e-02 3.762870e-01 6.701496e-01
round(rexp.est2,3)
## Estimate Robust SE LL UL
## (Intercept) 0.000 0.000 0.000 0.000
## lnpopterr 2.293 0.165 2.007 2.620
## logGDPCAP 1.082 0.059 0.958 1.223
## logUSMilitaryAid 1.043 0.017 1.010 1.077
## logArea 0.888 0.046 0.805 0.980
## Polity2 1.035 0.013 1.006 1.064
## CINC 0.000 0.000 0.000 0.000
## WOPOL 0.502 0.059 0.376 0.670
#########################
#Robustness tests
###########################
#R5: Model 2, plus a dummy for Muslim majority countries (Muslim, column, N)
summary(r5 <- glm.nb(incidentsusvictims_L1~ lnpopterr + logGDPCAP+logUSMilitaryAid+logArea+Polity2+CINC+WOPOL+Muslim, data = HD1))
##
## Call:
## glm.nb(formula = incidentsusvictims_L1 ~ lnpopterr + logGDPCAP +
## logUSMilitaryAid + logArea + Polity2 + CINC + WOPOL + Muslim,
## data = HD1, init.theta = 0.1217510624, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.8910 -0.5414 -0.4078 -0.2845 4.7509
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -13.65018 1.17280 -11.639 < 2e-16 ***
## lnpopterr 0.84773 0.07520 11.273 < 2e-16 ***
## logGDPCAP 0.08786 0.05439 1.615 0.106232
## logUSMilitaryAid 0.04573 0.01018 4.492 7.06e-06 ***
## logArea -0.09892 0.05239 -1.888 0.059014 .
## Polity2 0.01610 0.01329 1.212 0.225514
## CINC -16.87832 4.50725 -3.745 0.000181 ***
## WOPOL -0.79303 0.12472 -6.359 2.04e-10 ***
## Muslim -0.58924 0.18723 -3.147 0.001649 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.1218) family taken to be 1)
##
## Null deviance: 1276.4 on 3190 degrees of freedom
## Residual deviance: 1003.3 on 3182 degrees of freedom
## (1809 observations deleted due to missingness)
## AIC: 3168.2
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.1218
## Std. Err.: 0.0101
##
## 2 x log-likelihood: -3148.2380
cov.r5 <- vcovHC(r5, type="HC0")
std.err5r <- sqrt(diag(cov.r5))
r.est5r <- cbind(Estimate= coef(r5), "Robust SE" = std.err5r,
"Pr(>|z|)" = 2 * pnorm(abs(coef(r5)/std.err5r), lower.tail=FALSE),
LL = coef(r5) - 1.96 * std.err5r,
UL = coef(r5) + 1.96 * std.err5r)
round(r.est5r,3)
## Estimate Robust SE Pr(>|z|) LL UL
## (Intercept) -13.650 1.529 0.000 -16.648 -10.653
## lnpopterr 0.848 0.075 0.000 0.700 0.995
## logGDPCAP 0.088 0.056 0.115 -0.021 0.197
## logUSMilitaryAid 0.046 0.016 0.003 0.015 0.076
## logArea -0.099 0.050 0.047 -0.197 -0.001
## Polity2 0.016 0.014 0.234 -0.010 0.043
## CINC -16.878 4.699 0.000 -26.088 -7.668
## WOPOL -0.793 0.111 0.000 -1.011 -0.575
## Muslim -0.589 0.200 0.003 -0.981 -0.198
#R6: Model 2, plus the measure for physical integrity (PhysicalIntregity, column O)
summary(r6 <- glm.nb(incidentsusvictims_L1 ~ lnpopterr + logGDPCAP+logUSMilitaryAid+logArea+Polity2+CINC+WOPOL+PhysicalIntegrity, data = HD1))
##
## Call:
## glm.nb(formula = incidentsusvictims_L1 ~ lnpopterr + logGDPCAP +
## logUSMilitaryAid + logArea + Polity2 + CINC + WOPOL + PhysicalIntegrity,
## data = HD1, init.theta = 0.1554427257, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0709 -0.5294 -0.3797 -0.2523 4.8271
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -9.689757 1.182380 -8.195 2.50e-16 ***
## lnpopterr 0.545868 0.076304 7.154 8.44e-13 ***
## logGDPCAP 0.378425 0.057240 6.611 3.81e-11 ***
## logUSMilitaryAid 0.031358 0.009531 3.290 0.0010 **
## logArea -0.112510 0.050855 -2.212 0.0269 *
## Polity2 0.057515 0.012418 4.632 3.63e-06 ***
## CINC -6.370060 3.984075 -1.599 0.1098
## WOPOL -0.603430 0.116202 -5.193 2.07e-07 ***
## PhysicalIntegrity -0.381811 0.037069 -10.300 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.1554) family taken to be 1)
##
## Null deviance: 1452.4 on 3175 degrees of freedom
## Residual deviance: 1018.2 on 3167 degrees of freedom
## (1824 observations deleted due to missingness)
## AIC: 3055.4
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.1554
## Std. Err.: 0.0135
##
## 2 x log-likelihood: -3035.4170
cov.r6 <- vcovHC(r6, type="HC0")
std.err6r <- sqrt(diag(cov.r6))
r.est6r <- cbind(Estimate= coef(r6), "Robust SE" = std.err6r,
"Pr(>|z|)" = 2 * pnorm(abs(coef(r6)/std.err6r), lower.tail=FALSE),
LL = coef(r6) - 1.96 * std.err6r,
UL = coef(r6) + 1.96 * std.err6r)
round(r.est6r,3)
## Estimate Robust SE Pr(>|z|) LL UL
## (Intercept) -9.690 1.513 0.000 -12.655 -6.725
## lnpopterr 0.546 0.103 0.000 0.344 0.748
## logGDPCAP 0.378 0.058 0.000 0.265 0.492
## logUSMilitaryAid 0.031 0.014 0.024 0.004 0.059
## logArea -0.113 0.047 0.017 -0.205 -0.020
## Polity2 0.058 0.012 0.000 0.033 0.082
## CINC -6.370 3.150 0.043 -12.544 -0.196
## WOPOL -0.603 0.103 0.000 -0.805 -0.402
## PhysicalIntegrity -0.382 0.056 0.000 -0.492 -0.272
#R7: Model 2, plus the measure of American alliances (Alliances, column P)
summary(r7 <- glm.nb(incidentsusvictims_L1 ~ lnpopterr + logGDPCAP+logUSMilitaryAid+logArea+Polity2+CINC+WOPOL+Alliances, data = HD1))
##
## Call:
## glm.nb(formula = incidentsusvictims_L1 ~ lnpopterr + logGDPCAP +
## logUSMilitaryAid + logArea + Polity2 + CINC + WOPOL + Alliances,
## data = HD1, init.theta = 0.1531532113, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0857 -0.5474 -0.3642 -0.2315 4.7014
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -11.319728 1.139087 -9.938 < 2e-16 ***
## lnpopterr 0.682184 0.072796 9.371 < 2e-16 ***
## logGDPCAP -0.043792 0.056550 -0.774 0.43871
## logUSMilitaryAid 0.029011 0.009853 2.944 0.00324 **
## logArea -0.074716 0.051764 -1.443 0.14891
## Polity2 -0.035418 0.013313 -2.660 0.00781 **
## CINC -4.078672 4.181093 -0.976 0.32931
## WOPOL -0.687654 0.118783 -5.789 7.07e-09 ***
## Alliances 1.905105 0.176295 10.806 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.1532) family taken to be 1)
##
## Null deviance: 1445.9 on 3190 degrees of freedom
## Residual deviance: 1015.1 on 3182 degrees of freedom
## (1809 observations deleted due to missingness)
## AIC: 3059.9
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.1532
## Std. Err.: 0.0132
##
## 2 x log-likelihood: -3039.9030
cov.r7 <- vcovHC(r7, type="HC0")
std.err7r <- sqrt(diag(cov.r7))
r.est7r <- cbind(Estimate= coef(r7), "Robust SE" = std.err7r,
"Pr(>|z|)" = 2 * pnorm(abs(coef(r7)/std.err7r), lower.tail=FALSE),
LL = coef(r7) - 1.96 * std.err7r,
UL = coef(r7) + 1.96 * std.err7r)
round(r.est7r,3)
## Estimate Robust SE Pr(>|z|) LL UL
## (Intercept) -11.320 1.222 0.000 -13.715 -8.924
## lnpopterr 0.682 0.065 0.000 0.554 0.810
## logGDPCAP -0.044 0.059 0.460 -0.160 0.073
## logUSMilitaryAid 0.029 0.014 0.039 0.002 0.057
## logArea -0.075 0.051 0.141 -0.174 0.025
## Polity2 -0.035 0.012 0.003 -0.059 -0.012
## CINC -4.079 4.366 0.350 -12.637 4.479
## WOPOL -0.688 0.118 0.000 -0.920 -0.456
## Alliances 1.905 0.190 0.000 1.533 2.277
# R4: Model 2, plus the count of terrorist incidents in the preceding year (incidentsusvictims, column Q).
oldw <- getOption("warn")
options(warn = -1)
summary(r8 <- glm.nb(incidentsusvictims_L1 ~ lnpopterr + logGDPCAP+logUSMilitaryAid+logArea+Polity2+CINC+WOPOL+incidentsusvictims, data = HD1))
##
## Call:
## glm.nb(formula = incidentsusvictims_L1 ~ lnpopterr + logGDPCAP +
## logUSMilitaryAid + logArea + Polity2 + CINC + WOPOL + incidentsusvictims,
## data = HD1, init.theta = 0.172158494, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.0750 -0.4962 -0.3835 -0.2911 3.9132
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -10.960285 1.099651 -9.967 < 2e-16 ***
## lnpopterr 0.659393 0.070444 9.360 < 2e-16 ***
## logGDPCAP 0.083624 0.051066 1.638 0.1015
## logUSMilitaryAid 0.017221 0.009467 1.819 0.0689 .
## logArea -0.124084 0.050001 -2.482 0.0131 *
## Polity2 0.020219 0.011711 1.726 0.0843 .
## CINC -8.722123 3.919141 -2.226 0.0260 *
## WOPOL -0.586888 0.113355 -5.177 2.25e-07 ***
## incidentsusvictims 0.580635 0.022024 26.363 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.1722) family taken to be 1)
##
## Null deviance: 1535.9 on 3187 degrees of freedom
## Residual deviance: 1043.7 on 3179 degrees of freedom
## (1812 observations deleted due to missingness)
## AIC: 3029.5
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.1722
## Std. Err.: 0.0152
## Warning while fitting theta: alternation limit reached
##
## 2 x log-likelihood: -3009.4900
cov.r8 <- vcovHC(r8, type="HC0")
std.err8r <- sqrt(diag(cov.r8))
r.est8r <- cbind(Estimate= coef(r8), "Robust SE" = std.err8r,
"Pr(>|z|)" = 2 * pnorm(abs(coef(r8)/std.err8r), lower.tail=FALSE),
LL = coef(r8) - 1.96 * std.err8r,
UL = coef(r8) + 1.96 * std.err8r)
round(r.est8r,3)
## Estimate Robust SE Pr(>|z|) LL UL
## (Intercept) -10.960 1.184 0.000 -13.281 -8.639
## lnpopterr 0.659 0.073 0.000 0.516 0.803
## logGDPCAP 0.084 0.053 0.114 -0.020 0.187
## logUSMilitaryAid 0.017 0.014 0.203 -0.009 0.044
## logArea -0.124 0.045 0.006 -0.213 -0.035
## Polity2 0.020 0.013 0.122 -0.005 0.046
## CINC -8.722 3.295 0.008 -15.180 -2.264
## WOPOL -0.587 0.117 0.000 -0.815 -0.358
## incidentsusvictims 0.581 0.011 0.000 0.558 0.603
options(warn = oldw)