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)