##
## Analysis: Democratic peace theory
##

####################
## Produce tables ##
####################

# only select relevant variables
conflicts <- select(states.active, actor, conflictid, regime, max_prev) 

# Number of wars in which each regime type has used SV (per SV category)
conflict.numbers <- rbind(
  cbind( # create new row with sums
    Regime="democracy",
    #SV0=length(unique(filter(conflicts, regime=="democracy", max_prev==0)$conflictid)),
    SV1=length(unique(filter(conflicts, regime=="democracy", max_prev==1)$conflictid)),
    SV2=length(unique(filter(conflicts, regime=="democracy", max_prev==2)$conflictid)),
    SV3=length(unique(filter(conflicts, regime=="democracy", max_prev==3)$conflictid)),
    Total.withSV=length(unique(filter(conflicts, regime=="democracy", max_prev>0)$conflictid))
  ),
  cbind( # create new row with sums
    Regime="autocracy",
    #SV0=length(unique(filter(conflicts, regime=="autocracy", max_prev==0)$conflictid)),
    SV1=length(unique(filter(conflicts, regime=="autocracy", max_prev==1)$conflictid)),
    SV2=length(unique(filter(conflicts, regime=="autocracy", max_prev==2)$conflictid)),
    SV3=length(unique(filter(conflicts, regime=="autocracy", max_prev==3)$conflictid)),
    Total.withSV=length(unique(filter(conflicts, regime=="autocracy", max_prev>0)$conflictid))
  ),
  cbind( # create new row with sums
    Regime="anocracy",
    #SV0=length(unique(filter(conflicts, regime=="anocracy", max_prev==0)$conflictid)),
    SV1=length(unique(filter(conflicts, regime=="anocracy", max_prev==1)$conflictid)),
    SV2=length(unique(filter(conflicts, regime=="anocracy", max_prev==2)$conflictid)),
    SV3=length(unique(filter(conflicts, regime=="anocracy", max_prev==3)$conflictid)),
    Total.withSV=length(unique(filter(conflicts, regime=="anocracy", max_prev>0)$conflictid))
  )
)

# Number of regimes/actors/governments (per regime), that have used SV (per SV category)
actor.numbers <- rbind(
  cbind( # create new row with sums
    Regime="democracy",
    #SV0=length(unique(filter(conflicts, regime=="democracy", max_prev==0)$conflictid)),
    SV1=length(unique(filter(conflicts, regime=="democracy", max_prev==1)$actor)),
    SV2=length(unique(filter(conflicts, regime=="democracy", max_prev==2)$actor)),
    SV3=length(unique(filter(conflicts, regime=="democracy", max_prev==3)$actor)),
    Total.withsV=length(unique(filter(conflicts, regime=="democracy", max_prev>0)$actor))
  ),
  cbind( # create new row with sums
    Regime="autocracy",
    #SV0=length(unique(filter(conflicts, regime=="autocracy", max_prev==0)$conflictid)),
    SV1=length(unique(filter(conflicts, regime=="autocracy", max_prev==1)$actor)),
    SV2=length(unique(filter(conflicts, regime=="autocracy", max_prev==2)$actor)),
    SV3=length(unique(filter(conflicts, regime=="autocracy", max_prev==3)$actor)),
    Total.withSV=length(unique(filter(conflicts, regime=="autocracy", max_prev>0)$actor))
  ),
  cbind( # create new row with sums
    Regime="anocracy",
    #SV0=length(unique(filter(conflicts, regime=="anocracy", max_prev==0)$conflictid)),
    SV1=length(unique(filter(conflicts, regime=="anocracy", max_prev==1)$actor)),
    SV2=length(unique(filter(conflicts, regime=="anocracy", max_prev==2)$actor)),
    SV3=length(unique(filter(conflicts, regime=="anocracy", max_prev==3)$actor)),
    Total.withSV=length(unique(filter(conflicts, regime=="anocracy", max_prev>0)$actor)),
    Total.all=length(unique(filter(conflicts, regime=="anocracy")$actor))
  )
)

# Democracies that have engaged in wartime sexual violence
View(states[states$xpolity_highest==1 & 
              states$max_prev > 0 & 
              states$conflictyear==1,]
     [c("year", "actor", "max_prev", "xpolity", "location")])



##################################
## ORDINAL LOGISTIC REGRESSIONS ##
##################################

##
## Linear model
##

model.dempeace.linear <- clm(max_prev ~
                               # ethnic hatred
                               epr_ethwar + 
                               # state capacity
                               taxgdp.mwz + 
                               # military coherence
                               troop_quality_log.mwz + 
                               # feminist theory
                               wdi.fertility.mwz + 
                               # democratic peace theory
                               xpolity.mwz + 
                               # control variables
                               duration.mwz +
                               type_intrastate + type_internationalized + 
                               wdi.pop.density.log.mwz + mountain_log.mwz +
                               max_prev_lag.dummy,# + population_log.mwz, 
                             link="logit", data=states.active)
                           
# Coefficients and summary statistics
summary(model.dempeace.linear)
summary(model.dempeace.linear)$coefficients
model.dempeace.linear$nobs # Number of observations
exp(model.dempeace.linear$coefficients) %>% data.frame() # Odd Ratios (in a data frame)
#exp(confint(model.dempeace.linear, type = "Wald")) # confidence intervals for OR with Wald statistic
logLik(model.dempeace.linear)*-2        # -2 log-likelihood
AIC(model.dempeace.linear)              # AIC
nagelkerke(model.dempeace.linear$formula, states.active)

# predict probabilities (with model w/o centered variables)
model.dempeace.linear.uncentred <- clm(max_prev ~
                                         epr_ethwar + taxgdp + troop_quality_log + wdi.fertility + xpolity + 
                                         duration + type_intrastate + type_internationalized +
                                         wdi.pop.density.log + mountain_log + max_prev_lag.dummy,
                                       link="logit", data=states.active)
newData <- expand.grid(
  epr_ethwar = as.factor(0),#mean(as.numeric(states.active$pitf_ethwar), na.rm = T),
  taxgdp = mean(states.active$taxgdp, na.rm = T),
  troop_quality_log = mean(states.active$troop_quality_log, na.rm = T),
  wdi.fertility = mean(states.active$wdi.fertility, na.rm = T),
  xpolity = unique(states.active$xpolity),
  duration = mean(states.active$duration, na.rm = T),
  type_intrastate = as.factor(1),#mean(as.numeric(states.active$type_intrastate), na.rm = T),
  type_internationalized = as.factor(0),#mean(as.numeric(states.active$type_internationalized), na.rm = T),
  max_prev_lag.dummy = as.factor(1),
  wdi.pop.density.log = mean(states.active$wdi.pop.density.log, na.rm = T),
  mountain_log = mean(states.active$mountain_log, na.rm = T))
tmp.fit <- cbind(newData, predict(model.dempeace.linear.uncentred, newdata=newData, interval=TRUE, type="prob")$fit)
tmp.fit <- tmp.fit[c("xpolity", "0", "1", "2", "3")]
tmp.fit.melt <- melt(tmp.fit, id=c("xpolity"))

tmp.lwr <- cbind(newData, predict(model.dempeace.linear.uncentred, newdata=newData, interval=TRUE, type="prob")$lwr)
tmp.lwr <- tmp.lwr[c("xpolity", "0", "1", "2", "3")]
tmp.lwr.melt <- melt(tmp.lwr, id=c("xpolity"))

tmp.upr <- cbind(newData, predict(model.dempeace.linear.uncentred, newdata=newData, interval=TRUE, type="prob")$upr)
tmp.upr <- tmp.upr[c("xpolity", "0", "1", "2", "3")]
tmp.upr.melt <- melt(tmp.upr, id=c("xpolity"))

ggplot(tmp.fit.melt, aes(x = xpolity, y= value, colour=variable)) + 
  geom_ribbon(aes(ymin=tmp.lwr.melt$value, ymax=tmp.upr.melt$value, linetype=NA, fill=variable), alpha=.3) +
  labs(title="", x="Democracy (xpolity)", y="Probability") +
  scale_fill_discrete(name="Prevalence of\nSexual Violence",
                      labels=c("None", "Isolated", "Widespread", "Massive")) +
  scale_colour_discrete(name="Prevalence of\nSexual Violence",
                        labels=c("None", "Isolated", "Widespread", "Massive")) +
  geom_line() +
  theme_MA()

# show probabilities
View(unique(tmp.fit))


##
## Parabolic model
##

model.dempeace.parabolic <- clm(max_prev ~
                               # ethnic hatred
                               epr_ethwar + 
                               # state capacity
                               taxgdp.mwz + 
                               # military coherence
                               troop_quality_log.mwz + 
                               # feminist theory
                               wdi.fertility.mwz + 
                               # democratic peace theory
                               xpolity.mwz + xpolity_quadratic.mwz +
                               # control variables
                               duration.mwz +
                               type_intrastate + type_internationalized + 
                               wdi.pop.density.log.mwz + mountain_log.mwz +
                               max_prev_lag.dummy,
                             link="logit", data=states.active)

model <- model.dempeace.parabolic
summary(model)
model$nobs # Number of observations
logLik(model)*-2        # -2 log-likelihood
nagelkerke(model$formula, states.active)



##
## Threshold effect model
##

model.dempeace.threshold <- clm(max_prev ~
                               # ethnic hatred
                               epr_ethwar + 
                               # state capacity
                               taxgdp.mwz + 
                               # military coherence
                               troop_quality_log.mwz + 
                               # feminist theory
                               wdi.fertility.mwz + 
                               # democratic peace theory
                               xpolity_highest + xpolity_lowest +
                               # control variables
                               duration.mwz +
                               type_intrastate + type_internationalized + 
                               wdi.pop.density.log.mwz + mountain_log.mwz +
                               max_prev_lag.dummy,
                             link="logit", data=states.active)

model <- model.dempeace.threshold
summary(model)
model$nobs # Number of observations
logLik(model)*-2        # -2 log-likelihood
nagelkerke(model$formula, states.active)


# predict probabilities (with model w/o centered variables)
model.dempeace.threshold.uncentred <- clm(max_prev ~
                                         epr_ethwar + taxgdp + troop_quality_log + wdi.fertility + 
                                         xpolity_highest + xpolity_lowest +
                                         duration + type_intrastate + type_internationalized +
                                         wdi.pop.density.log + mountain_log + max_prev_lag.dummy,
                                       link="logit", data=states.active)
newData <- expand.grid(
  epr_ethwar = as.factor(0),#mean(as.numeric(states.active$pitf_ethwar), na.rm = T),
  taxgdp = mean(states.active$taxgdp, na.rm = T),
  troop_quality_log = mean(states.active$troop_quality_log, na.rm = T),
  wdi.fertility = mean(states.active$wdi.fertility, na.rm = T),
  xpolity_highest = as.factor(c(0, 1)),
  xpolity_lowest = as.factor(0),
  duration = mean(states.active$duration, na.rm = T),
  type_intrastate = as.factor(1),#mean(as.numeric(states.active$type_intrastate), na.rm = T),
  type_internationalized = as.factor(0),#mean(as.numeric(states.active$type_internationalized), na.rm = T),
  max_prev_lag.dummy = as.factor(1),
  wdi.pop.density.log = mean(states.active$wdi.pop.density.log, na.rm = T),
  mountain_log = mean(states.active$mountain_log, na.rm = T))
tmp.fit <- cbind(newData, predict(model.dempeace.threshold.uncentred, newdata=newData, interval=TRUE, type="prob")$fit)
tmp.fit <- tmp.fit[c("xpolity_highest", "0", "1", "2", "3")]
tmp.fit.melt <- melt(tmp.fit, id=c("xpolity_highest"))

tmp.lwr <- cbind(newData, predict(model.dempeace.threshold.uncentred, newdata=newData, interval=TRUE, type="prob")$lwr)
tmp.lwr <- tmp.lwr[c("xpolity_highest", "0", "1", "2", "3")]
tmp.lwr.melt <- melt(tmp.lwr, id=c("xpolity_highest"))

tmp.upr <- cbind(newData, predict(model.dempeace.threshold.uncentred, newdata=newData, interval=TRUE, type="prob")$upr)
tmp.upr <- tmp.upr[c("xpolity_highest", "0", "1", "2", "3")]
tmp.upr.melt <- melt(tmp.upr, id=c("xpolity_highest"))

ggplot(tmp.fit.melt, aes(x = variable, y=value, fill=xpolity_highest)) + 
  geom_bar(stat = "identity", position = "dodge") +
  labs(title="", x="Prevalence of Sexual Violence", y="Probability") +
  scale_fill_discrete(name="Democracy\n(xpolity)",
                      labels=c("Non-Democracy", "Full Democracy")) +
  scale_colour_discrete(name="Democracy\n(xpolity)",
                        labels=c("Non-Democracy", "Full Democracy")) +
  geom_errorbar(aes(ymin=tmp.lwr.melt$value, ymax=tmp.upr.melt$value), stat = "identity", position =position_dodge(1), width = 0.25) +
  theme_MA()

# show probabilities
View(unique(tmp.fit))
