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)