commit 9d6f834859cd7c3fa6ebf4057372609c5d71f2ee
parent f2ec3e6d5188be0057b4a8928f6f8dfaa6e2727e
Author: eamoncaddigan <eamon.caddigan@gmail.com>
Date: Sat, 29 Aug 2015 22:58:37 -0400
Seperated post. calculation from plotting
Diffstat:
1 file changed, 23 insertions(+), 12 deletions(-)
diff --git a/ggPostPlot.R b/ggPostPlot.R
@@ -3,8 +3,9 @@
library(ggplot2)
library(dplyr)
-ggPosteriorPredictive <- function(modelData, codaObject,
- x1Level, x2Level = NA, x3Level = NA) {
+# Find the responses and the posterior predictions for the given factor levels
+calcPosteriorPredictive <- function(modelData, codaObject,
+ x1Level, x2Level = NA, x3Level = NA) {
mcmcMat <- as.matrix(codaObject)
# Separate the observed counts for the given levels
@@ -32,7 +33,7 @@ ggPosteriorPredictive <- function(modelData, codaObject,
}
}
sigma <- mcmcMat[, paste0("sigma[", x1Level, "]")]
-
+
# Find the HDI and mean of the posterior probabilities of each of the y levels
# (Stolen from Kruschke)
chainLength <- nrow(mcmcMat)
@@ -46,25 +47,35 @@ ggPosteriorPredictive <- function(modelData, codaObject,
outHdi <- apply(outProb, 2, HDIofMCMC)
outMean <- apply(outProb, 2, mean, na.rm=TRUE)
- # Create the plot
+ # Create a data.frame with the posterior data
plotData <- cellData %>%
count(response) %>%
- left_join(data_frame(response = 1:max(y)), ., by = "response")
- plotData <- plotData %>%
+ left_join(data_frame(response = 1:max(y)), ., by = "response") %>%
mutate(n = ifelse(is.na(n), 0, n),
response_proportion = n / sum(n),
posterior_mean = outMean,
posterior_hdi_low = outHdi[1,],
posterior_hdi_high = outHdi[2,])
+ return(plotData)
+}
+
+
+# Plot the responses and posterior predictions for the given factor levels
+ggPosteriorPredictive <- function(modelData, codaObject,
+ x1Level, x2Level = NA, x3Level = NA) {
+ plotData <- calcPosteriorPredictive(modelData, codaObject,
+ x1Level, x2Level, x3Level)
+ numLevels = length(unique(plotData[["response"]]))
+
p <- ggplot(plotData, aes(x = response, y = response_proportion)) +
- geom_bar(stat = "identity", binwidth=1, color="white") +
+ geom_bar(stat = "identity", binwidth=1, color="white", fill="skyblue") +
geom_point(aes(y = posterior_mean),
- size=3, color="red") +
+ size=5, color="red") +
geom_errorbar(aes(ymin = posterior_hdi_low, ymax = posterior_hdi_high),
- size=2, color="red", width=0) +
- scale_x_continuous(breaks = 1:6) + # XXX hardcoded 6
+ size=3, color="red", width=0) +
+ scale_x_continuous(breaks = 1:numLevels) +
scale_y_continuous(limits = c(0, 1)) +
- xlab("Response") +
- ylab("Proportion of Responses")
+ xlab("Response Level") +
+ ylab("Proportion")
return(p)
}