0

Below is the code I am using to generate predicted probabilities after running my glmer regression model. The code below works fine at generating a graph that shows the predicted probabilities along all the values of intra-EU trade share in the data frame.

#GRAPH OF SITC 1 intra-eu trade share 
tmpdat1 <- multi.sanctions.bust1.full@frame[, c("caseid", "lndist", "lnlaggdpt", "duration", "partnercode", "lnlaggdpp", "lagtradeopenP", "lageutradeshare100", "lagsitc1100","EUmem", "nobust1", "nobust1sq", "nobust1cb", "colonial", "lagtradecontrol1", "SEM")]

jvalues1eu <- with(multi.sanctions.bust1.full, seq(from = min(multi.sanctions.bust1.full@frame[["lageutradeshare100"]]), to = max(multi.sanctions.bust1.full@frame[["lageutradeshare100"]]), length.out = 100))

pp1eu <- lapply(jvalues1eu, function(j) {
  tmpdat1$lageutradeshare100 <- j
  predict(multi.sanctions.bust1.full, newdata = tmpdat1, type = "response", re.form = NULL)
})

#subset(multi.sanctions.bust1.full@frame,allbuster1==0)

#sapply(pp1[c(0, 1, 2, 3, 4, 5)], mean)

plotdat1eu <- t(sapply(pp1eu, function(x) {
    c(M = mean(x), Med = median(x), quantile(x, c(0.25, 0.75), na.rm = TRUE), (mean(x)-(2*sd(x))),
      (mean(x)+(2*sd(x))))
}))

plotdat1eu <- as.data.frame(cbind(plotdat1eu, jvalues1eu))
colnames(plotdat1eu) <- c("PredictedProbabilityMean", "PredProbMedian", "quartile1", "quartile3", "lowersd", "uppersd", "lageutradeshare100")
head(plotdat1eu)
tail(plotdat1eu)

sitc1eu <- ggplot() + geom_line(data=plotdat1eu, aes(x = lageutradeshare100, y = PredictedProbabilityMean), size = 2, color="blue") + 
  geom_ribbon(data=plotdat1eu, aes(x = lageutradeshare100, ymin = lowersd, ymax = uppersd),
              fill = "grey50", alpha=.5) +
  ylim(c(-0.25, 0.55)) + 
  geom_hline(yintercept=0) +
    geom_rug(data=subset(multi.sanctions.bust1.full@frame,allbuster1==0), aes(x=lageutradeshare100), color="black", size=1.0, sides="b", alpha= 3/4, length = unit(0.05, "npc")) +
    geom_rug(data=subset(multi.sanctions.bust1.full@frame,allbuster1==1), aes(x=lageutradeshare100), color="red", size=1.0, sides="b", alpha = 1) +  
 theme(panel.grid.major = element_line(colour = "gray", linetype = "dotted"), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.y = element_text(size=14, face="bold"), axis.title.x = element_text(size=14, face="bold")) +
  xlab("Intra-EU Trade Share") + 
  ylab("Probability of Sanctions Busting") + facet_wrap(~partnercode,nrow=4)

sitc1eu

What I would like to do is get the facet_wrap(~partnercode,nrow=4) to show different predicted probabilities depending on the 27 countries in my data set. The rug plot changes, but the line does not (nor did I expect it to change). I am not sure where or how to alter the code. I am guessing that I need to generate predicted probabilities for each country or partnercode?

Any help greatly appreciated.

Keith
  • 193
  • 1
  • 10
  • 1
    Maybe you should provide a reproducible example of your data you are trying to plot in order people can assist you. (see: https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) – dc37 Jan 04 '20 at 06:48
  • You need to include partnercode as a variable in plotdat1eu and the facet_wrap should work out. I cannot see where partnercode being used anywhere in the code above, and it will really if you provide some example data and the relevant code which highlights your problem – StupidWolf Jan 04 '20 at 09:39
  • @StupidWolf partnercode is in the first line of code in tmpdat1. I'll work on editing my question. – Keith Jan 04 '20 at 23:02

0 Answers0