R/savage_dickey.R

In LudvigOlsen/LRO.utilities: Utility functions and addins

#   __________________ #< 89d7a9f84fe8faae45b2ec8cd64cbcb0 ># __________________
#   Savage Dickey Bayes Factor                                              ####

#' @title Calculate Bayes Factor using Savage-Dickey density ratio.
#' @description Creates plot of prior and posterior distributions and calculates the
#'  Bayes factor at the point of interest (Q).
#' @details Bayes factors are calculated using the polspline package.
#'
#'  BF10 <- dlogspline(Q, logspline(post)) /
#'  dlogspline(Q,logspline(prior))
#'
#'  BF01 <- dlogspline(Q, logspline(prior)) /
#'  dlogspline(Q,logspline(post))
#' @param post Samples from posterior distribution (Numeric)
#' @param prior Samples from prior distribution (Numeric)
#' @param Q Point on x-axis for calculating Bayes factor.
#' @param xlab Label for x-axis
#' @param ylab Label for y-axis
#' @param plot Create density plot (Logical)
#' @param print_plot (Logical)
#' @return List with ggplot2 object (optional), BF10 and BF01,
#'  Maximum a posteriori and priori at point Q
#' @author Benjamin Hugh Zachariae
#' @author Ludvig Renbo Olsen
#' @export
#' @examples
#' # Sampling from two different gaussian distributions
#' prior <- rnorm(1000, mean=0, sd=1)
#' posterior <- rnorm(1000, mean=2, sd=2)
#'
#' # Calculating BF and generating plot of the given distributions
#' s_d <- savage_dickey(posterior, prior, Q = 0, plot = TRUE)
#' @importFrom polspline dlogspline logspline
#' @importFrom dplyr '%>%'
savage_dickey <- function(post, prior, Q,
                          xlab = 'Parameter values',
                          ylab = 'Density',
                          plot = TRUE,
                          print_plot = FALSE){


##  .................. #< e3616350b525535a2bfdf402f2e6b31e ># ..................
##  Check arguments                                                         ####


  stopifnot(is.logical(plot),
            is.logical(print_plot),
            is.numeric(post) | is.integer(post),
            is.numeric(prior) | is.integer(prior),
            is.numeric(Q) | is.integer(Q))

  if ( !isTRUE(plot) && isTRUE(print_plot)){
    message("'print_plot' has no effect when 'plot' is FALSE")
  }

##  Maximum
  maximum_a_posteriori = dlogspline(Q, logspline(post))
  maximum_a_priori = dlogspline(Q, logspline(prior))

##  .................. #< fc3ac4bea4b5054bdda78ec4479a8730 ># ..................
##  Bayes factor                                                            ####

  BF10 <-
    maximum_a_posteriori /
    maximum_a_priori
  BF01 <-
    maximum_a_priori /
    maximum_a_posteriori


##  .................. #< 597ccd506c1a6e9cb52cd7fca23de8d8 ># ..................
##  Plot                                                                    ####


  if (isTRUE(plot)){

    # Gather posterior and prior
    results <- data.frame("posterior" = post,
                          "prior" = prior )

    plot_ <- results %>%
      tidyr::gather("postprior", "gathered") %>%
      ggplot2::ggplot() +
      ggplot2::geom_density(ggplot2::aes(gathered, fill = postprior), alpha = .5) +
      ggplot2::labs(x=xlab, y=ylab) +
      ggplot2::geom_vline(xintercept=Q, color='steelblue2')

    if(isTRUE(print_plot)){
      print(plot_)
    }


##  .................. #< 4310b4db83cb41a819f09cb31930189a ># ..................
##  Return list                                                             ####

    return(list("post_prior_plot" = plot_, "BF10" = BF10, "BF01" = BF01,
                "Maximum a posteriori at Q" = maximum_a_posteriori,
                "Maximum a priori at Q" = maximum_a_priori))

  } else {

    return(list("BF10" = BF10, "BF01" = BF01,
                "Maximum a posteriori at Q" = maximum_a_posteriori,
                "Maximum a priori at Q" = maximum_a_priori))
  }


}