replication_noah_carl_submission_2016.R

### ANALYTIC REPLICATIONS OF NOAH CARL'S Explaining Terrorism Threat Level Across Western Countries


# libs --------------------------------------------------------------------

library(pacman)
p_load(XLConnect, kirkegaard, psych, weights, magrittr, effsize, lsr, compute.es, MASS)


# data --------------------------------------------------------------------

d_main = readWorksheetFromFile(file = "Terrorism Threat.xlsx", sheet = 2)

rownames(d_main) = d_main$country

# derived variables ------------------------------------------

#dummies
d_main$west = as.logical(d_main$west)
d_main$europe = as.logical(d_main$europe)
d_main$turkey = as.logical(d_main$turkey)

#muslim 2015 estimation
d_main$muslim15 = d_main$muslim10 + .25 * (d_main$muslim30 - d_main$muslim10)

#death count
d_main$deaths2 = log(1 + d_main$deaths)

#log gdp
d_main$gdp_log = log(d_main$gdpcap14)

#factors
d_main$any_num = d_main$any
d_main$any = as.factor(d_main$any)
d_main$part = as.factor(d_main$part)
d_main$terror_f = as.factor(d_main$terror)

#winzor turkey's value
d_main["Turkey", "muslim15"] = 19.075

# subsets -----------------------------------------------------------------

d_west = d_main[d_main$west, ]
d_euro = d_main[d_main$europe, ]

# descriptive -------------------------------------------------------------

#threat level
psych::describe(d_west$terror)

#muslims
psych::describe(d_west$muslim15)

#deaths
psych::describe(d_west$deaths)

#any
table(d_west$any)

#against isis
table(d_west$part)

# cors --------------------------------------------------------------------

wtd.cors(d_west)

cor.test(d_west$terror, d_west$muslim15)
cor.test(d_west$terror, d_west$deaths2)


# simple lm's ---------------------------------------------------------------------
lm("terror ~ poly(muslim15, 1)", data = d_west) %>% summary()
lm("terror ~ poly(muslim15, 2)", data = d_west) %>% summary()

lm("terror ~ poly(deaths2, 1)", data = d_west) %>% summary()
lm("terror ~ poly(deaths2, 2)", data = d_west) %>% summary()

# SMD ---------------------------------------------------------------------

#calculate the d
SMD_matrix(d_west$terror, d_west$any)
cohen.d(d = d_west$terror, f = as.factor(d_west$any))
SMD_matrix(d_west$terror, d_west$part)
cohen.d(d = d_west$terror, f = as.factor(d_west$part))

#try to get a p value
des(d = cohen.d(d = d_west$terror, f = d_west$any)$estimate, n.1 = 21, n.2 = 7)
des(d = cohen.d(d = d_west$terror, f = d_west$part)$estimate, n.1 = 21, n.2 = 7)


# complex lm's ------------------------------------------------------------

#all west
d_betas = MOD_APSLM(dependent = "terror", predictors = c("muslim15", "any", "deaths2", "part", "gdp_log", "unemp14", "ineq0911"), data = d_west)[[1]] %>% round(2)

#the models we want
d_betas[c(8:10, 108, 112, 113), ]

#all oecd
d_betas = MOD_APSLM(dependent = "terror", predictors = c("muslim15", "any", "deaths2", "part", "gdp_log", "unemp14", "ineq0911"), data = d_main)[[1]] %>% round(2)

#the models we want
d_betas[c(8:10, 108, 112, 113), ]

#europe only
d_betas = MOD_APSLM(dependent = "terror", predictors = c("muslim15", "any", "deaths2", "part", "gdp_log", "unemp14", "ineq0911"), data = d_euro)[[1]] %>% round(2)

#the models we want
d_betas[c(8:10, 108, 112, 113), ]

#logistic
polr("terror_f ~ muslim15 + any", data = d_west) %>% summary()
polr("terror_f ~ muslim15 + deaths2", data = d_west) %>% summary()
polr("terror_f ~ muslim15 + part", data = d_west) %>% summary()

#lasso
d_betas_lasso = MOD_repeat_cv_glmnet(df = d_west, dependent = "terror", predictors = c("muslim15", "any", "deaths2", "part", "gdp_log", "unemp14", "ineq0911"), runs = 500)
MOD_summarize_models(d_betas_lasso)