R/posterior_freq_checks.R

In duncan-clark/Blolog: Bayesian Inference for LOLOG Network Models

#' @export
#function to do frequentist style checks:
#Draws a value from the posterior, then simulates from that value and compares to the observed

posterior_freq_checks <- function(sample,
                             formula_rhs,
                             gof_formula_rhs,
                             net,
                             draws,
                             title = "Goodness of Fit Graph",
                             xlab = NULL,
                             ylab = NULL){
  
  formula <- as.formula(paste("net ~ ",formula_rhs))
  
  #sample
  theta = unlist(sample(sample,1))
  tmp = lolog::createLatentOrderLikelihood(formula,theta = theta)
  #simulate
  test_stats <- lapply(1:draws,function(x){
    sim = tmp$generateNetwork()$network
    sim = lolog::as.network(sim)
    formula_gof <- as.formula(paste("sim ~ ",gof_formula_rhs))
    stats = lolog::calculateStatistics(formula_gof)
    return(stats)
  })
  
  #make into dataframe for plotting
  test_stats <- do.call(rbind,test_stats)
  test_stats <- as.data.frame(test_stats)
  test_stats$sim <- seq(1,draws,1)
  test_stats <- reshape2::melt(test_stats,id.vars = c("sim"))
  
  
  #make the obs stats into a dataframe for plotting:
  formula_gof <- as.formula(paste("net ~ ",gof_formula_rhs))
  obs_stats <- lolog::calculateStatistics(formula_gof)
  obs_stats <- as.data.frame(obs_stats)
  obs_stats$variable <- rownames(obs_stats)
  obs_stats$sim  <- 1
  names(obs_stats) <- c("value","variable","sim")
  
  if(!is.null(ylab)){
    ylab <- "Statistic"
  }
  
  #plot 
  gg_line <- ggplot(data = test_stats,aes(x = variable,y=value,group = sim))+
    geom_line(alpha = lineAlpha, size = lineSize)+
    geom_line(data = obs_stats,aes(x = variable,y=value,group = sim),color = "red", size = lineSize)+
    ggtitle(title)+
    ggplot2::theme_bw() + ggplot2::ylab(ylab) + ggplot2::xlab("") + 
    ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, 
                                                       hjust = 1), panel.grid.major.x = ggplot2::element_blank(), 
                   panel.grid.minor.x = ggplot2::element_blank())
  
  gg_box <- ggplot(data = test_stats,aes(x = variable,y=value))+
    geom_boxplot()+
    geom_point(data = obs_stats,aes(x = variable,y=value),color = "red")+
    ggtitle(title)+
    ggplot2::theme_bw() + ggplot2::ylab(ylab) + ggplot2::xlab("") + 
    ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, 
                                                       hjust = 1), panel.grid.major.x = ggplot2::element_blank(), 
                   panel.grid.minor.x = ggplot2::element_blank())
  
  return(list(gg_line = gg_line,
              gg_box = gg_box))
}