Functions for model/data simulations to learn about a binomial proportion

model_data_simulation.R

model_data_simulation <- function(P_vector,
                                  mystat=sum,
                                  N=100,
                                  N_sim=10000){
  require(ggplot2)
  TH <- theme(
    plot.title = element_text(
      colour = "red",
      size = 14,
      hjust = 0.5,
      vjust = 0.8,
      angle = 0
    )
  )
  name_mystat <- deparse(substitute(mystat))
  one_sim <- function(mystat){
    P <- sample(P_vector, size=1)
    S <- mystat(rbinom(N, size=1, prob=P))
    c(P, S)
  }
  S <- data.frame(t(replicate(N_sim, one_sim(mystat))))
  names(S) <- c("P", "Statistic")
  p <- ggplot(S, aes(P, Statistic)) +
    geom_jitter(size=0.2) +
    geom_smooth() +
    ggtitle(paste("Model/Data Simulation \n Statistic:",
                  name_mystat, ", N =", N)) + TH
  print(p)
  S
}

inference_plot <- function(S, value_statistic){
  require(ggplot2)
  require(dplyr)
  TH <- theme(
    plot.title = element_text(
      colour = "red",
      size = 14,
      hjust = 0.5,
      vjust = 0.8,
      angle = 0
    )
  )
  S_new <- filter(S, Statistic==value_statistic)
  limits <- quantile(S_new$P, c(0.25, 0.75))
  p <- ggplot(S_new, aes(P)) +
         geom_density(adjust=2) +
         ggtitle(paste("Posterior of P when Statistic =",
                       value_statistic,
                       "\n 50% Interval Estimate")) + TH +
       geom_vline(xintercept = limits, color="blue")
  print(p)
}