##
## Impact of Timing on sexual violence (interm, postc, conflictyear)
##

# predict basic model with ordinal
model <- clm(max_prev ~ interm + postc, link="logit", data=states)

interm.newData <- expand.grid(
  max_prev=as.factor(c(0,1,2,3)),
  interm=as.factor(c(0,1)),
  postc=as.factor(c(1)))
                   
postc.newData <- expand.grid(
  max_prev=as.factor(c(0,1,2,3)),
  interm=as.factor(c(0)),
  postc=as.factor(c(0,1)))

conflictyear.newData <- expand.grid(
  max_prev=as.factor(c(0,1,2,3)),
  conflictyear=as.factor(c(0,1)))

# probabilities for interim phase
interm.fit <- cbind(interm.newData, probability=predict(model, newdata=interm.newData, interval=TRUE, type="prob")$fit)
interm.lwr <- cbind(interm.newData, probability=predict(model, newdata=interm.newData, interval=TRUE, type="prob")$lwr)
interm.upr <- cbind(interm.newData, probability=predict(model, newdata=interm.newData, interval=TRUE, type="prob")$upr)

# probabilities for post-conflict phase
postc.fit <- cbind(postc.newData, probability=predict(model, newdata=postc.newData, interval=TRUE, type="prob")$fit)
postc.lwr <- cbind(postc.newData, probability=predict(model, newdata=postc.newData, interval=TRUE, type="prob")$lwr)
postc.upr <- cbind(postc.newData, probability=predict(model, newdata=postc.newData, interval=TRUE, type="prob")$upr)

# probabilities for conflict phase
conflict.model <- clm(max_prev ~ conflictyear, link="logit", data=states)
conflict.fit <- cbind(conflictyear.newData, probability=predict(conflict.model, newdata=conflictyear.newData, interval=TRUE, type="prob")$fit)
conflict.lwr <- cbind(conflictyear.newData, probability=predict(conflict.model, newdata=conflictyear.newData, interval=TRUE, type="prob")$lwr)
conflict.upr <- cbind(conflictyear.newData, probability=predict(conflict.model, newdata=conflictyear.newData, interval=TRUE, type="prob")$upr)


# graphics
g1 <- ggplot(interm.fit, aes(x=max_prev, y=probability, fill=interm)) + 
  geom_bar(stat = "identity", position="dodge") +
  geom_errorbar(aes(ymin=interm.lwr$probability, ymax=interm.upr$probability), stat = "identity", position =position_dodge(1), width = 0.25) +
  scale_fill_discrete(name="Intermediate\nphase") +
  labs(title="", x="Prevalence of SV", y="Probability") +
  theme_MA()

g2 <- ggplot(postc.fit, aes(x=max_prev, y=probability, fill=postc)) + 
  geom_bar(stat = "identity", position="dodge") +
  geom_errorbar(aes(ymin=postc.lwr$probability, ymax=postc.upr$probability), stat = "identity", position =position_dodge(1), width = 0.25) +
  scale_fill_discrete(name="Post-conflict\nphase") +
  labs(title="", x="Prevalence of SV", y="Probability") +
  theme_MA()

g3 <- ggplot(conflict.fit, aes(x=max_prev, y=probability, fill=conflictyear)) + 
  geom_bar(stat = "identity", position="dodge") +
  scale_y_continuous(expand = c(0,0), limits = c(0,1)) +
  geom_errorbar(aes(ymin=conflict.lwr$probability, ymax=conflict.upr$probability), stat = "identity", position = position_dodge(1), width = 0.25) +
  scale_fill_discrete(name="Conflict\nphase") +
  labs(title="", x="Prevalence of SV", y="Probability") +
  theme_MA()


library(gridExtra)
grid.arrange(g1, g2, g3, ncol=3)
