#############################################################
## Analysis State capacity / Military coherence (H2a, H2b) ##
#############################################################

model.capacity <- 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,
                      link="logit", data=states.active)

              
# Coefficients and summary statistics
summary(model.capacity)

# predict probabilities (with model w/o centered variables)
model.capacity.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 = unique(states.active$troop_quality_log),
  wdi.fertility = mean(states.active$wdi.fertility, na.rm = T),
  xpolity = mean(states.active$xpolity, na.rm = T),
  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.capacity.uncentred, newdata=newData, interval=TRUE, type="prob")$fit)
tmp.fit <- tmp.fit[c("troop_quality_log", "0", "1", "2", "3")]
tmp.fit.melt <- melt(tmp.fit, id=c("troop_quality_log"))

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

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

ggplot(tmp.fit.melt, aes(x = troop_quality_log, 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="Troop quality (log)", 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))