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
#MODEL3
summary(m3 <- glm.nb(incidentsusvictims_L1 ~ lnpopterr + logGDPCAP+logUSMilitaryAid+logArea+Polity2+CINC+WOSOC, data = HD1))
##
## Call:
## glm.nb(formula = incidentsusvictims_L1 ~ lnpopterr + logGDPCAP +
## logUSMilitaryAid + logArea + Polity2 + CINC + WOSOC, data = HD1,
## init.theta = 0.1232447969, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.8725 -0.5396 -0.4233 -0.3024 4.4541
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -14.39150 1.18498 -12.145 < 2e-16 ***
## lnpopterr 0.80680 0.07593 10.626 < 2e-16 ***
## logGDPCAP 0.23817 0.05814 4.097 4.19e-05 ***
## logUSMilitaryAid 0.03580 0.01032 3.470 0.000520 ***
## logArea -0.13973 0.05351 -2.611 0.009022 **
## Polity2 0.03172 0.01238 2.562 0.010421 *
## CINC -15.98454 4.54251 -3.519 0.000433 ***
## WOSOC -0.51623 0.10679 -4.834 1.34e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.1232) family taken to be 1)
##
## Null deviance: 1212.90 on 2963 degrees of freedom
## Residual deviance: 961.51 on 2956 degrees of freedom
## (2036 observations deleted due to missingness)
## AIC: 3045.4
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.1232
## Std. Err.: 0.0104
##
## 2 x log-likelihood: -3027.4040
cov.m3 <- vcovHC(m3, type="HC0")
std.err3 <- sqrt(diag(cov.m3))
r.est3 <- cbind(Estimate= coef(m3), "Robust SE" = std.err3,
"Pr(>|z|)" = 2 * pnorm(abs(coef(m3)/std.err3), lower.tail=FALSE),
LL = coef(m3) - 1.96 * std.err3,
UL = coef(m3) + 1.96 * std.err3)
round(r.est3,3)
## Estimate Robust SE Pr(>|z|) LL UL
## (Intercept) -14.392 1.474 0.000 -17.280 -11.503
## lnpopterr 0.807 0.071 0.000 0.668 0.946
## logGDPCAP 0.238 0.060 0.000 0.120 0.356
## logUSMilitaryAid 0.036 0.016 0.022 0.005 0.066
## logArea -0.140 0.050 0.005 -0.237 -0.042
## Polity2 0.032 0.014 0.023 0.004 0.059
## CINC -15.985 4.196 0.000 -24.208 -7.761
## WOSOC -0.516 0.128 0.000 -0.768 -0.265
with(m3, cbind(res.deviance = deviance, df = df.residual,
p = pchisq(deviance, df.residual, lower.tail=FALSE)))
## res.deviance df p
## [1,] 961.5055 2956 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.
s3 <- deltamethod(list(~ exp(x1), ~ exp(x2), ~ exp(x3), ~ exp(x4),~ exp(x5),~ exp(x6),~ exp(x7),~ exp(x8)), coef(m3), cov.m3)
## exponentiate old estimates dropping the p values
rexp.est3 <- exp(r.est3[, -3])
## replace SEs with estimates for exponentiated coefficients
rexp.est3[, "Robust SE"] <- s3
rexp.est3
## Estimate Robust SE LL UL
## (Intercept) 5.621469e-07 8.285794e-07 3.127473e-08 0.0000101043
## lnpopterr 2.240728e+00 1.588374e-01 1.950066e+00 2.5747136423
## logGDPCAP 1.268920e+00 7.644994e-02 1.127587e+00 1.4279678931
## logUSMilitaryAid 1.036453e+00 1.619765e-02 1.005187e+00 1.0686912448
## logArea 8.695955e-01 4.328924e-02 7.887564e-01 0.9587196281
## Polity2 1.032229e+00 1.435666e-02 1.004470e+00 1.0607554143
## CINC 1.142888e-07 4.795272e-07 3.065563e-11 0.0004260859
## WOSOC 5.967646e-01 7.659034e-02 4.640402e-01 0.7674506850
round(rexp.est3,3)
## Estimate Robust SE LL UL
## (Intercept) 0.000 0.000 0.000 0.000
## lnpopterr 2.241 0.159 1.950 2.575
## logGDPCAP 1.269 0.076 1.128 1.428
## logUSMilitaryAid 1.036 0.016 1.005 1.069
## logArea 0.870 0.043 0.789 0.959
## Polity2 1.032 0.014 1.004 1.061
## CINC 0.000 0.000 0.000 0.000
## WOSOC 0.597 0.077 0.464 0.767
#########################
#Robustness tests
###########################
#R9: Model 3, plus a dummy for Muslim majority countries (Muslim, column, N)
summary(r9 <- glm.nb(incidentsusvictims_L1~ lnpopterr + logGDPCAP+logUSMilitaryAid+logArea+Polity2+CINC+WOSOC+Muslim, data = HD1))
##
## Call:
## glm.nb(formula = incidentsusvictims_L1 ~ lnpopterr + logGDPCAP +
## logUSMilitaryAid + logArea + Polity2 + CINC + WOSOC + Muslim,
## data = HD1, init.theta = 0.1249031411, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.8868 -0.5466 -0.4164 -0.2930 4.0935
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -14.84388 1.20785 -12.290 < 2e-16 ***
## lnpopterr 0.82423 0.07643 10.784 < 2e-16 ***
## logGDPCAP 0.27533 0.05935 4.639 3.50e-06 ***
## logUSMilitaryAid 0.03755 0.01035 3.628 0.000286 ***
## logArea -0.12207 0.05349 -2.282 0.022475 *
## Polity2 0.01451 0.01337 1.086 0.277657
## CINC -17.29179 4.58532 -3.771 0.000163 ***
## WOSOC -0.63091 0.11087 -5.690 1.27e-08 ***
## Muslim -0.61675 0.19162 -3.219 0.001288 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.1249) family taken to be 1)
##
## Null deviance: 1222.0 on 2963 degrees of freedom
## Residual deviance: 959.1 on 2955 degrees of freedom
## (2036 observations deleted due to missingness)
## AIC: 3038.1
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.1249
## Std. Err.: 0.0106
##
## 2 x log-likelihood: -3018.1120
cov.r9 <- vcovHC(r9, type="HC0")
std.err9r <- sqrt(diag(cov.r9))
r.est9r <- cbind(Estimate= coef(r9), "Robust SE" = std.err9r,
"Pr(>|z|)" = 2 * pnorm(abs(coef(r9)/std.err9r), lower.tail=FALSE),
LL = coef(r9) - 1.96 * std.err9r,
UL = coef(r9) + 1.96 * std.err9r)
round(r.est9r,3)
## Estimate Robust SE Pr(>|z|) LL UL
## (Intercept) -14.844 1.478 0.000 -17.740 -11.948
## lnpopterr 0.824 0.073 0.000 0.682 0.967
## logGDPCAP 0.275 0.062 0.000 0.154 0.397
## logUSMilitaryAid 0.038 0.015 0.014 0.008 0.067
## logArea -0.122 0.048 0.011 -0.217 -0.027
## Polity2 0.015 0.015 0.320 -0.014 0.043
## CINC -17.292 4.463 0.000 -26.039 -8.545
## WOSOC -0.631 0.135 0.000 -0.896 -0.366
## Muslim -0.617 0.216 0.004 -1.040 -0.194
#R10: Model 3, plus the measure for physical integrity (PhysicalIntregity, column O)
summary(r10 <- glm.nb(incidentsusvictims_L1 ~ lnpopterr + logGDPCAP+logUSMilitaryAid+logArea+Polity2+CINC+WOSOC+PhysicalIntegrity, data = HD1))
##
## Call:
## glm.nb(formula = incidentsusvictims_L1 ~ lnpopterr + logGDPCAP +
## logUSMilitaryAid + logArea + Polity2 + CINC + WOSOC + PhysicalIntegrity,
## data = HD1, init.theta = 0.1586658595, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0698 -0.5356 -0.3895 -0.2475 5.2138
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -10.843110 1.212527 -8.943 < 2e-16 ***
## lnpopterr 0.538670 0.077734 6.930 4.22e-12 ***
## logGDPCAP 0.493331 0.061319 8.045 8.60e-16 ***
## logUSMilitaryAid 0.031168 0.009763 3.192 0.001411 **
## logArea -0.127041 0.052013 -2.442 0.014586 *
## Polity2 0.048210 0.012426 3.880 0.000105 ***
## CINC -6.979140 4.101545 -1.702 0.088833 .
## WOSOC -0.268154 0.105850 -2.533 0.011298 *
## PhysicalIntegrity -0.388710 0.038534 -10.087 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.1587) family taken to be 1)
##
## Null deviance: 1388.9 on 2953 degrees of freedom
## Residual deviance: 971.3 on 2945 degrees of freedom
## (2046 observations deleted due to missingness)
## AIC: 2930.3
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.1587
## Std. Err.: 0.0140
##
## 2 x log-likelihood: -2910.3010
cov.r10 <- vcovHC(r10, type="HC0")
std.err10r <- sqrt(diag(cov.r10))
r.est10r <- cbind(Estimate= coef(r10), "Robust SE" = std.err10r,
"Pr(>|z|)" = 2 * pnorm(abs(coef(r10)/std.err10r), lower.tail=FALSE),
LL = coef(r10) - 1.96 * std.err10r,
UL = coef(r10) + 1.96 * std.err10r)
round(r.est10r,3)
## Estimate Robust SE Pr(>|z|) LL UL
## (Intercept) -10.843 1.631 0.000 -14.039 -7.647
## lnpopterr 0.539 0.108 0.000 0.327 0.750
## logGDPCAP 0.493 0.069 0.000 0.359 0.628
## logUSMilitaryAid 0.031 0.014 0.025 0.004 0.058
## logArea -0.127 0.046 0.006 -0.218 -0.036
## Polity2 0.048 0.013 0.000 0.022 0.074
## CINC -6.979 3.239 0.031 -13.328 -0.630
## WOSOC -0.268 0.155 0.083 -0.571 0.035
## PhysicalIntegrity -0.389 0.054 0.000 -0.494 -0.284
#R11: Model 3 plus the measure of American alliances (Alliances, column P)
summary(r11 <- glm.nb(incidentsusvictims_L1 ~ lnpopterr + logGDPCAP+logUSMilitaryAid+logArea+Polity2+CINC+WOSOC+Alliances, data = HD1))
##
## Call:
## glm.nb(formula = incidentsusvictims_L1 ~ lnpopterr + logGDPCAP +
## logUSMilitaryAid + logArea + Polity2 + CINC + WOSOC + Alliances,
## data = HD1, init.theta = 0.1607404402, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.1419 -0.5442 -0.3670 -0.2315 3.9855
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -12.22448 1.16162 -10.524 < 2e-16 ***
## lnpopterr 0.67436 0.07452 9.050 < 2e-16 ***
## logGDPCAP 0.13148 0.05936 2.215 0.0268 *
## logUSMilitaryAid 0.01920 0.01004 1.912 0.0558 .
## logArea -0.11896 0.05323 -2.235 0.0254 *
## Polity2 -0.03248 0.01338 -2.427 0.0152 *
## CINC -2.61641 4.22387 -0.619 0.5356
## WOSOC -0.75195 0.10834 -6.940 3.91e-12 ***
## Alliances 2.10000 0.18025 11.650 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.1607) family taken to be 1)
##
## Null deviance: 1401.9 on 2963 degrees of freedom
## Residual deviance: 960.5 on 2955 degrees of freedom
## (2036 observations deleted due to missingness)
## AIC: 2913.2
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.1607
## Std. Err.: 0.0139
##
## 2 x log-likelihood: -2893.1630
cov.r11 <- vcovHC(r11, type="HC0")
std.err11r <- sqrt(diag(cov.r11))
r.est11r <- cbind(Estimate= coef(r11), "Robust SE" = std.err11r,
"Pr(>|z|)" = 2 * pnorm(abs(coef(r11)/std.err11r), lower.tail=FALSE),
LL = coef(r11) - 1.96 * std.err11r,
UL = coef(r11) + 1.96 * std.err11r)
round(r.est11r,3)
## Estimate Robust SE Pr(>|z|) LL UL
## (Intercept) -12.224 1.170 0.000 -14.519 -9.930
## lnpopterr 0.674 0.067 0.000 0.544 0.805
## logGDPCAP 0.131 0.061 0.032 0.011 0.251
## logUSMilitaryAid 0.019 0.014 0.155 -0.007 0.046
## logArea -0.119 0.049 0.015 -0.215 -0.023
## Polity2 -0.032 0.012 0.009 -0.057 -0.008
## CINC -2.616 4.128 0.526 -10.708 5.475
## WOSOC -0.752 0.128 0.000 -1.003 -0.501
## Alliances 2.100 0.175 0.000 1.756 2.444
# R12: Model 1, plus the count of terrorist incidents in the preceding year (incidentsusvictims, column Q).
oldw <- getOption("warn")
options(warn = -1)
summary(r12 <- glm.nb(incidentsusvictims_L1 ~ lnpopterr + logGDPCAP+logUSMilitaryAid+logArea+Polity2+CINC+WOSOC+incidentsusvictims, data = HD1))
##
## Call:
## glm.nb(formula = incidentsusvictims_L1 ~ lnpopterr + logGDPCAP +
## logUSMilitaryAid + logArea + Polity2 + CINC + WOSOC + incidentsusvictims,
## data = HD1, init.theta = 0.1752908864, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.0963 -0.5059 -0.3961 -0.3038 4.5269
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -11.799132 1.123271 -10.504 < 2e-16 ***
## lnpopterr 0.649995 0.071983 9.030 < 2e-16 ***
## logGDPCAP 0.191951 0.055076 3.485 0.000492 ***
## logUSMilitaryAid 0.013618 0.009665 1.409 0.158825
## logArea -0.145071 0.051242 -2.831 0.004638 **
## Polity2 0.019543 0.011821 1.653 0.098292 .
## CINC -9.020898 3.999775 -2.255 0.024111 *
## WOSOC -0.398292 0.101743 -3.915 9.05e-05 ***
## incidentsusvictims 0.573930 0.022290 25.748 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.1753) family taken to be 1)
##
## Null deviance: 1467.00 on 2963 degrees of freedom
## Residual deviance: 996.99 on 2955 degrees of freedom
## (2036 observations deleted due to missingness)
## AIC: 2907.9
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.1753
## Std. Err.: 0.0158
## Warning while fitting theta: alternation limit reached
##
## 2 x log-likelihood: -2887.8560
cov.r12 <- vcovHC(r12, type="HC0")
std.err12r <- sqrt(diag(cov.r12))
r.est12r <- cbind(Estimate= coef(r12), "Robust SE" = std.err12r,
"Pr(>|z|)" = 2 * pnorm(abs(coef(r12)/std.err12r), lower.tail=FALSE),
LL = coef(r12) - 1.96 * std.err12r,
UL = coef(r12) + 1.96 * std.err12r)
round(r.est12r,3)
## Estimate Robust SE Pr(>|z|) LL UL
## (Intercept) -11.799 1.292 0.000 -14.332 -9.267
## lnpopterr 0.650 0.078 0.000 0.497 0.803
## logGDPCAP 0.192 0.058 0.001 0.078 0.306
## logUSMilitaryAid 0.014 0.014 0.340 -0.014 0.042
## logArea -0.145 0.046 0.002 -0.235 -0.055
## Polity2 0.020 0.012 0.108 -0.004 0.043
## CINC -9.021 3.338 0.007 -15.563 -2.479
## WOSOC -0.398 0.141 0.005 -0.675 -0.122
## incidentsusvictims 0.574 0.012 0.000 0.551 0.597
options(warn = oldw)