antivax-attitudes

Reanalyses of data from Horne, Powell, Hummel & Holyoak (2015)
git clone https://git.eamoncaddigan.net/antivax-attitudes.git
Log | Files | Refs | README | LICENSE

ggPostPlot.R (3237B)

      1 # ggplot2-based plots for some of the things
      2 
      3 library(ggplot2)
      4 library(dplyr)
      5 
      6 # Find the responses and the posterior predictions for the given factor levels
      7 calcPosteriorPredictive <- function(modelData, codaObject,
      8                                     x1Level, x2Level = NA, x3Level = NA) {
      9   mcmcMat <- as.matrix(codaObject)
     10   
     11   # Separate the observed counts for the given levels
     12   # XXX - hard-coded factor names because laziness
     13   cellData <- filter(modelData, as.numeric(as.factor(question)) == x1Level)
     14   if (!is.na(x2Level)) {
     15     cellData <- filter(cellData, as.numeric(as.factor(intervention)) == x2Level)
     16   }
     17   if (!is.na(x3Level)) {
     18     cellData <- filter(cellData, as.numeric(as.factor(interval)) == x3Level)
     19   }
     20   y <- cellData[["response"]]
     21   
     22   # Find the mean, standard deviation, and thresholds
     23   mu <- mcmcMat[, "b0"] + mcmcMat[, paste0("b1[", x1Level, "]")]
     24   if (!is.na(x2Level)) {
     25     mu <- mu + mcmcMat[, paste0("b2[", x2Level, "]")] + 
     26       mcmcMat[, paste0("b1b2[", x1Level, ",", x2Level, "]")]
     27   }
     28   if (!is.na(x3Level)) {
     29     mu <- mu + mcmcMat[, paste0("b3[", x3Level, "]")] + 
     30       mcmcMat[, paste0("b1b3[", x1Level, ",", x3Level, "]")]
     31     if (!is.na(x2Level)) {
     32       mu <- mu + mcmcMat[, paste0("b1b2b3[", x1Level, ",", x2Level, ",", x3Level, "]")]
     33     }
     34   }
     35   sigma <- mcmcMat[, paste0("sigma[", x1Level, "]")]
     36   
     37   # Find the HDI and mean of the posterior probabilities of each of the y levels
     38   # (Stolen from Kruschke)
     39   chainLength <- nrow(mcmcMat)
     40   outProb <- matrix(0, nrow=chainLength, ncol=max(y))
     41   for (stepIdx in 1:chainLength) {
     42     threshCumProb <- pnorm((mcmcMat[stepIdx, 
     43                                     paste0("thresh[", x1Level, ",",1:(max(y)-1),"]")]
     44                             - mu[stepIdx]) / sigma[stepIdx])
     45     outProb[stepIdx,] <- c(threshCumProb, 1) - c(0, threshCumProb)
     46   }
     47   outHdi <- apply(outProb, 2, HDIofMCMC)
     48   outMean <- apply(outProb, 2, mean, na.rm=TRUE)
     49   
     50   # Create a data.frame with the posterior data
     51   plotData <- cellData %>%
     52     count(response) %>% 
     53     left_join(data_frame(response = 1:max(y)), ., by = "response") %>%
     54     mutate(n = ifelse(is.na(n), 0, n),
     55            response_proportion = n / sum(n),
     56            posterior_mean = outMean,
     57            posterior_hdi_low = outHdi[1,],
     58            posterior_hdi_high = outHdi[2,])
     59   return(plotData)
     60 }
     61 
     62 
     63 # Plot the responses and posterior predictions for the given factor levels
     64 ggPosteriorPredictive <- function(modelData, codaObject, 
     65                                   x1Level, x2Level = NA, x3Level = NA) {
     66   plotData <- calcPosteriorPredictive(modelData, codaObject,
     67                                       x1Level, x2Level, x3Level)
     68   numLevels = length(unique(plotData[["response"]]))
     69 
     70   p <- ggplot(plotData, aes(x = response, y = response_proportion)) + 
     71     geom_bar(stat = "identity", binwidth=1, color="white", fill="skyblue") + 
     72     geom_point(aes(y = posterior_mean), 
     73                size=3, color="red") + 
     74     geom_errorbar(aes(ymin = posterior_hdi_low, ymax = posterior_hdi_high), 
     75                   size=2, color="red", width=0) + 
     76     scale_x_continuous(breaks = 1:numLevels) +
     77     scale_y_continuous(limits = c(0, 1)) +
     78     xlab("Response Level") + 
     79     ylab("Proportion")
     80   return(p)
     81 }