model.ethwar <- 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)
summary(model.ethwar)
model.ethwar$nobs # Number of observations
logLik(model.ethwar)*-2        # -2 log-likelihood
nagelkerke(model.ethwar$formula, states.active)


model.secess <- clm(max_prev ~
                      # ethnic hatred
                      epr_secesswar + 
                      # 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)
summary(model.secess)
model.ethwar$nobs # Number of observations
logLik(model.secess)*-2        # -2 log-likelihood
nagelkerke(model.secess$formula, states.active)

model.secess <- clm(max_prev ~
                      # ethnic hatred
                      epr_secesswar + 
                      # 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)
summary(model.secess)
model.ethwar$nobs # Number of observations
logLik(model.secess)*-2        # -2 log-likelihood
nagelkerke(model.secess$formula, states.active)


model.genoc <- clm(max_prev ~
                      # ethnic hatred
                      pitf_genocwar + 
                      # 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)
summary(model.genoc)
model.ethwar$nobs # Number of observations
logLik(model.genoc)*-2        # -2 log-likelihood
nagelkerke(model.genoc$formula, states.active)


model.genoc.mag <- clm(max_prev ~
                     # ethnic hatred
                     pitf_genocwar_deathmag.mwz + 
                     # 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)

model <- model.genoc.mag
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.capacity.uncentred <- clm(max_prev ~
                                  pitf_genocwar_deathmag + 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(
  pitf_genocwar_deathmag = unique(states$pitf_genocwar_deathmag),
  taxgdp = mean(states.active$taxgdp, na.rm = T),
  troop_quality_log = mean(states$troop_quality_log, na.rm = T),
  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("pitf_genocwar_deathmag", "0", "1", "2", "3")]
tmp.fit.melt <- melt(tmp.fit, id=c("pitf_genocwar_deathmag"))

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

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

ggplot(tmp.fit.melt, aes(x = pitf_genocwar_deathmag, y= value, colour=variable)) + 
  geom_line() +
  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_ribbon(aes(ymin=tmp.lwr.melt$value, ymax=tmp.upr.melt$value, linetype=NA, fill=variable), alpha=.3) +
  labs(title="", x="Magnitude of Genocide", y="Probability") +
  theme_MA()
