## Auto-Regressive

model.autoregressive <- 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 = 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(c(0,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.autoregressive, newdata=newData, interval=TRUE, type="prob")$fit)
tmp.fit <- tmp.fit[c("max_prev_lag.dummy", "0", "1", "2", "3")]
tmp.fit.melt <- melt(tmp.fit, id=c("max_prev_lag.dummy"))

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

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

View(unique(tmp.fit))

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()