R/Bayesian-function.R

In lcw68/G3proj: BIOS 735 Project

#' Bayesian Logistic Regression Model Prediction on MCMC Beta Samples
#'
#' This function computes predicted probabilities based on MCMC Beta samples
#' generated by BLRM_fit_mwg function.
#'
#' @param Y.test a vector of binary data
#' @param X.test a design matrix that has the same number of columns with the dimension of Beta
#' @param MCsamples a matrix of MCMC sampled Betas stacked in rows
#'
#' @return a data frame of input binary data and corresponding predicted probabilities of being 1.
#'
#' @examples
#'
#' ## simulate data
#'
#' set.seed(1);
#' N  = 1000;
#' p  = 10;
#'
#' X.test    = matrix(data = rnorm(N*p), nrow=N, ncol=p)
#' beta_true = c(rep(1,p/2),rep(0,p/2))
#' eta       = X.test %*% beta_true
#' pi        = exp(eta) / (1 + exp(eta))
#' Y.test    = rbinom(N,1,pi)
#'
#' ## simulate MCMC beta samples (beta_true + N(0, 0.5) random errors)
#'
#' M = 1000
#' MCsamples = matrix(data = rep(beta_true, M), nrow = M, ncol = p, byrow = T) + matrix(data = rnorm(M*p, sd = 0.5), nrow=M, ncol=p)
#'
#' ## predict based on MCMC beta samples
#'
#' prediction = predict_BLRM(Y.test, X.test, MCsamples)
#'
#' ## Classification Metrics with cutoff = 0.5
#' caret::confusionMatrix(data = factor(prediction$data), reference = factor(ifelse(prediction$pred.prob > 0.5, 1, 0)), positive = "1")
#'
#' @export
predict_BLRM = function(Y.test, X.test, MCsamples){

  N = nrow(X.test)
  M = nrow(MCsamples)

  eta_mat = X.test %*% t(MCsamples)
  mu_mat = exp(eta_mat) / (1 + exp(eta_mat));

  z_mat = matrix(rbinom(N*M, 1, mu_mat), N, M)
  pred_z = rowMeans(z_mat)

  return(data.frame(data = Y.test, pred.prob = pred_z))
}