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
################Model1#####
summary(m1 <- glm.nb(incidentsusvictims_L1 ~ lnpopterr + logGDPCAP+logUSMilitaryAid+logArea+Polity2+CINC+WECON, data = HD1))
##
## Call:
## glm.nb(formula = incidentsusvictims_L1 ~ lnpopterr + logGDPCAP +
## logUSMilitaryAid + logArea + Polity2 + CINC + WECON, data = HD1,
## init.theta = 0.1215498483, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.9445 -0.5357 -0.4207 -0.3010 4.4536
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -13.36463 1.14589 -11.663 < 2e-16 ***
## lnpopterr 0.76339 0.07460 10.233 < 2e-16 ***
## logGDPCAP 0.23761 0.05735 4.143 3.42e-05 ***
## logUSMilitaryAid 0.03960 0.01010 3.920 8.87e-05 ***
## logArea -0.14153 0.05290 -2.676 0.00746 **
## Polity2 0.02446 0.01183 2.068 0.03868 *
## CINC -13.63610 4.42736 -3.080 0.00207 **
## WECON -0.79809 0.13405 -5.953 2.63e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.1215) family taken to be 1)
##
## Null deviance: 1261.69 on 3154 degrees of freedom
## Residual deviance: 999.11 on 3147 degrees of freedom
## (1845 observations deleted due to missingness)
## AIC: 3144.5
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.1215
## Std. Err.: 0.0102
##
## 2 x log-likelihood: -3126.5110
cov.m1 <- vcovHC(m1, type="HC0")
std.err1 <- sqrt(diag(cov.m1))
r.est1 <- cbind(Estimate= coef(m1), "Robust SE" = std.err1,
"Pr(>|z|)" = 2 * pnorm(abs(coef(m1)/std.err1), lower.tail=FALSE),
LL = coef(m1) - 1.96 * std.err1,
UL = coef(m1) + 1.96 * std.err1)
r.est1
## Estimate Robust SE Pr(>|z|) LL
## (Intercept) -13.36462713 1.34437538 2.756476e-23 -1.599960e+01
## lnpopterr 0.76338969 0.06855539 8.439280e-29 6.290211e-01
## logGDPCAP 0.23760729 0.06370351 1.915611e-04 1.127484e-01
## logUSMilitaryAid 0.03959710 0.01560174 1.114906e-02 9.017678e-03
## logArea -0.14152765 0.04669954 2.440678e-03 -2.330588e-01
## Polity2 0.02446225 0.01291673 5.824591e-02 -8.545379e-04
## CINC -13.63609659 4.09445923 8.672819e-04 -2.166124e+01
## WECON -0.79808625 0.13844794 8.188939e-09 -1.069444e+00
## UL
## (Intercept) -10.72965137
## lnpopterr 0.89775825
## logGDPCAP 0.36246618
## logUSMilitaryAid 0.07017651
## logArea -0.04999656
## Polity2 0.04977904
## CINC -5.61095650
## WECON -0.52672829
round(r.est1,3)
## Estimate Robust SE Pr(>|z|) LL UL
## (Intercept) -13.365 1.344 0.000 -16.000 -10.730
## lnpopterr 0.763 0.069 0.000 0.629 0.898
## logGDPCAP 0.238 0.064 0.000 0.113 0.362
## logUSMilitaryAid 0.040 0.016 0.011 0.009 0.070
## logArea -0.142 0.047 0.002 -0.233 -0.050
## Polity2 0.024 0.013 0.058 -0.001 0.050
## CINC -13.636 4.094 0.001 -21.661 -5.611
## WECON -0.798 0.138 0.000 -1.069 -0.527
with(m1, cbind(res.deviance = deviance, df = df.residual,
p = pchisq(deviance, df.residual, lower.tail=FALSE)))
## res.deviance df p
## [1,] 999.1088 3147 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.
s1 <- deltamethod(list(~ exp(x1), ~ exp(x2), ~ exp(x3), ~ exp(x4),~ exp(x5),~ exp(x6),~ exp(x7),~ exp(x8)), coef(m1), cov.m1)
## exponentiate old estimates dropping the p values
rexp.est1 <- exp(r.est1[, -3])
## replace SEs with estimates for exponentiated coefficients
rexp.est1[, "Robust SE"] <- s1
round(rexp.est1,3)
## Estimate Robust SE LL UL
## (Intercept) 0.000 0.000 0.000 0.000
## lnpopterr 2.146 0.147 1.876 2.454
## logGDPCAP 1.268 0.081 1.119 1.437
## logUSMilitaryAid 1.040 0.016 1.009 1.073
## logArea 0.868 0.041 0.792 0.951
## Polity2 1.025 0.013 0.999 1.051
## CINC 0.000 0.000 0.000 0.004
## WECON 0.450 0.062 0.343 0.591
#########################
#Robustness tests
###########################
#R1: Model 1, plus a dummy for Muslim majority countries (Muslim, column, N)
summary(r1 <- glm.nb(incidentsusvictims_L1~ lnpopterr + logGDPCAP+logUSMilitaryAid+logArea+Polity2+CINC+WECON+Muslim, data = HD1))
##
## Call:
## glm.nb(formula = incidentsusvictims_L1 ~ lnpopterr + logGDPCAP +
## logUSMilitaryAid + logArea + Polity2 + CINC + WECON + Muslim,
## data = HD1, init.theta = 0.1230881645, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.9665 -0.5421 -0.4132 -0.2996 4.1557
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -13.644597 1.162278 -11.740 < 2e-16 ***
## lnpopterr 0.772297 0.074982 10.300 < 2e-16 ***
## logGDPCAP 0.262225 0.058204 4.505 6.63e-06 ***
## logUSMilitaryAid 0.042719 0.010141 4.213 2.52e-05 ***
## logArea -0.125672 0.052808 -2.380 0.01732 *
## Polity2 0.007761 0.013046 0.595 0.55193
## CINC -14.511681 4.471467 -3.245 0.00117 **
## WECON -0.875147 0.136604 -6.406 1.49e-10 ***
## Muslim -0.509376 0.184632 -2.759 0.00580 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.1231) family taken to be 1)
##
## Null deviance: 1270.56 on 3154 degrees of freedom
## Residual deviance: 998.68 on 3146 degrees of freedom
## (1845 observations deleted due to missingness)
## AIC: 3139.4
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.1231
## Std. Err.: 0.0103
##
## 2 x log-likelihood: -3119.3980
cov.r1 <- vcovHC(r1, type="HC0")
std.err1r <- sqrt(diag(cov.r1))
r.est1r <- cbind(Estimate= coef(r1), "Robust SE" = std.err1r,
"Pr(>|z|)" = 2 * pnorm(abs(coef(r1)/std.err1r), lower.tail=FALSE),
LL = coef(r1) - 1.96 * std.err1r,
UL = coef(r1) + 1.96 * std.err1r)
round(r.est1r,3)
## Estimate Robust SE Pr(>|z|) LL UL
## (Intercept) -13.645 1.358 0.000 -16.306 -10.983
## lnpopterr 0.772 0.071 0.000 0.634 0.911
## logGDPCAP 0.262 0.066 0.000 0.133 0.391
## logUSMilitaryAid 0.043 0.015 0.005 0.013 0.073
## logArea -0.126 0.045 0.005 -0.214 -0.037
## Polity2 0.008 0.014 0.570 -0.019 0.035
## CINC -14.512 4.411 0.001 -23.156 -5.867
## WECON -0.875 0.143 0.000 -1.156 -0.595
## Muslim -0.509 0.209 0.015 -0.920 -0.099
#R2: Model 1, plus the measure for physical integrity (PhysicalIntregity, column O)
summary(r2 <- glm.nb(incidentsusvictims_L1 ~ lnpopterr + logGDPCAP+logUSMilitaryAid+logArea+Polity2+CINC+WECON+PhysicalIntegrity, data = HD1))
##
## Call:
## glm.nb(formula = incidentsusvictims_L1 ~ lnpopterr + logGDPCAP +
## logUSMilitaryAid + logArea + Polity2 + CINC + WECON + PhysicalIntegrity,
## data = HD1, init.theta = 0.1527578932, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0538 -0.5292 -0.3880 -0.2565 5.0295
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -10.25355 1.18389 -8.661 < 2e-16 ***
## lnpopterr 0.51706 0.07646 6.762 1.36e-11 ***
## logGDPCAP 0.47476 0.05999 7.915 2.48e-15 ***
## logUSMilitaryAid 0.03130 0.00956 3.274 0.001062 **
## logArea -0.12337 0.05139 -2.401 0.016370 *
## Polity2 0.04108 0.01196 3.434 0.000595 ***
## CINC -5.88005 3.99877 -1.470 0.141436
## WECON -0.47824 0.13065 -3.660 0.000252 ***
## PhysicalIntegrity -0.36661 0.03777 -9.705 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.1528) family taken to be 1)
##
## Null deviance: 1424.4 on 3141 degrees of freedom
## Residual deviance: 1011.6 on 3133 degrees of freedom
## (1858 observations deleted due to missingness)
## AIC: 3040.3
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.1528
## Std. Err.: 0.0134
##
## 2 x log-likelihood: -3020.3080
cov.r2 <- vcovHC(r2, type="HC0")
std.err2r <- sqrt(diag(cov.r2))
r.est2r <- cbind(Estimate= coef(r2), "Robust SE" = std.err2r,
"Pr(>|z|)" = 2 * pnorm(abs(coef(r2)/std.err2r), lower.tail=FALSE),
LL = coef(r2) - 1.96 * std.err2r,
UL = coef(r2) + 1.96 * std.err2r)
round(r.est2r,3)
## Estimate Robust SE Pr(>|z|) LL UL
## (Intercept) -10.254 1.574 0.000 -13.338 -7.169
## lnpopterr 0.517 0.105 0.000 0.312 0.722
## logGDPCAP 0.475 0.067 0.000 0.344 0.605
## logUSMilitaryAid 0.031 0.014 0.029 0.003 0.059
## logArea -0.123 0.045 0.006 -0.212 -0.035
## Polity2 0.041 0.013 0.001 0.017 0.066
## CINC -5.880 3.196 0.066 -12.144 0.384
## WECON -0.478 0.133 0.000 -0.738 -0.218
## PhysicalIntegrity -0.367 0.056 0.000 -0.477 -0.256
#R3: Model 1, plus the measure of American alliances (Alliances, column P)
summary(r3 <- glm.nb(incidentsusvictims_L1 ~ lnpopterr + logGDPCAP+logUSMilitaryAid+logArea+Polity2+CINC+WECON+Alliances, data = HD1))
##
## Call:
## glm.nb(formula = incidentsusvictims_L1 ~ lnpopterr + logGDPCAP +
## logUSMilitaryAid + logArea + Polity2 + CINC + WECON + Alliances,
## data = HD1, init.theta = 0.1547452841, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.1514 -0.5462 -0.3626 -0.2336 4.0158
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -11.287085 1.127878 -10.007 < 2e-16 ***
## lnpopterr 0.629805 0.073150 8.610 < 2e-16 ***
## logGDPCAP 0.118604 0.059379 1.997 0.045780 *
## logUSMilitaryAid 0.027593 0.009899 2.787 0.005313 **
## logArea -0.112046 0.052557 -2.132 0.033017 *
## Polity2 -0.043343 0.013023 -3.328 0.000875 ***
## CINC -1.376650 4.161613 -0.331 0.740798
## WECON -0.895107 0.132726 -6.744 1.54e-11 ***
## Alliances 1.952493 0.175759 11.109 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.1547) family taken to be 1)
##
## Null deviance: 1438.48 on 3154 degrees of freedom
## Residual deviance: 999.83 on 3146 degrees of freedom
## (1845 observations deleted due to missingness)
## AIC: 3022
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.1547
## Std. Err.: 0.0133
##
## 2 x log-likelihood: -3001.9570
cov.r3 <- vcovHC(r3, type="HC0")
std.err3r <- sqrt(diag(cov.r3))
r.est3r <- cbind(Estimate= coef(r3), "Robust SE" = std.err3r,
"Pr(>|z|)" = 2 * pnorm(abs(coef(r3)/std.err3r), lower.tail=FALSE),
LL = coef(r3) - 1.96 * std.err3r,
UL = coef(r3) + 1.96 * std.err3r)
round(r.est3r,3)
## Estimate Robust SE Pr(>|z|) LL UL
## (Intercept) -11.287 1.092 0.000 -13.428 -9.147
## lnpopterr 0.630 0.065 0.000 0.503 0.756
## logGDPCAP 0.119 0.068 0.080 -0.014 0.251
## logUSMilitaryAid 0.028 0.014 0.043 0.001 0.054
## logArea -0.112 0.046 0.016 -0.203 -0.021
## Polity2 -0.043 0.012 0.000 -0.066 -0.021
## CINC -1.377 4.237 0.745 -9.682 6.928
## WECON -0.895 0.129 0.000 -1.147 -0.643
## Alliances 1.952 0.175 0.000 1.609 2.296
# R4: Model 1, plus the count of terrorist incidents in the preceding year (incidentsusvictims, column Q).
oldw <- getOption("warn")
options(warn = -1)
summary(r4 <- glm.nb(incidentsusvictims_L1 ~ lnpopterr + logGDPCAP+logUSMilitaryAid+logArea+Polity2+CINC+WECON+incidentsusvictims, data = HD1))
##
## Call:
## glm.nb(formula = incidentsusvictims_L1 ~ lnpopterr + logGDPCAP +
## logUSMilitaryAid + logArea + Polity2 + CINC + WECON + incidentsusvictims,
## data = HD1, init.theta = 0.1703013083, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.0327 -0.4963 -0.3943 -0.3030 4.1536
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -11.472717 1.098513 -10.444 < 2e-16 ***
## lnpopterr 0.627117 0.070887 8.847 < 2e-16 ***
## logGDPCAP 0.187020 0.054423 3.436 0.000589 ***
## logUSMilitaryAid 0.017165 0.009477 1.811 0.070110 .
## logArea -0.130900 0.050579 -2.588 0.009653 **
## Polity2 0.007262 0.011284 0.644 0.519855
## CINC -8.270395 3.914170 -2.113 0.034606 *
## WECON -0.482507 0.126114 -3.826 0.000130 ***
## incidentsusvictims 0.572814 0.022605 25.340 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.1703) family taken to be 1)
##
## Null deviance: 1511.2 on 3151 degrees of freedom
## Residual deviance: 1038.0 on 3143 degrees of freedom
## (1848 observations deleted due to missingness)
## AIC: 3012.2
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.1703
## Std. Err.: 0.0152
## Warning while fitting theta: alternation limit reached
##
## 2 x log-likelihood: -2992.2320
cov.r4 <- vcovHC(r4, type="HC0")
std.err4r <- sqrt(diag(cov.r4))
r.est4r <- cbind(Estimate= coef(r4), "Robust SE" = std.err4r,
"Pr(>|z|)" = 2 * pnorm(abs(coef(r4)/std.err4r), lower.tail=FALSE),
LL = coef(r4) - 1.96 * std.err4r,
UL = coef(r4) + 1.96 * std.err4r)
round(r.est4r,3)
## Estimate Robust SE Pr(>|z|) LL UL
## (Intercept) -11.473 1.214 0.000 -13.853 -9.092
## lnpopterr 0.627 0.075 0.000 0.481 0.774
## logGDPCAP 0.187 0.061 0.002 0.067 0.307
## logUSMilitaryAid 0.017 0.014 0.229 -0.011 0.045
## logArea -0.131 0.046 0.004 -0.221 -0.041
## Polity2 0.007 0.012 0.538 -0.016 0.030
## CINC -8.270 3.278 0.012 -14.696 -1.845
## WECON -0.483 0.137 0.000 -0.751 -0.214
## incidentsusvictims 0.573 0.012 0.000 0.550 0.596
options(warn = oldw)