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)