R/AdaPowCor.R

In abnormally-distributed/cvreg: Cross Validation and Robust Estimation Utilities

#' Bayesian Linear Regression with Adaptive Powered Correlation Prior
#'
#' @description The adaptive powered correlation prior extends the Zellner-Siow Cauchy g-prior by allowing the crossproduct of the
#' model matrix to be raised to powers other than -1. The power here will be referred to as "lambda". A lambda of 0 results
#' in an identity matrix, which results in a ridge-regression like prior. Under the ridge prior, all coefficients are shrunk towards
#' zero. Positive values of lambda adapt to collinearity by tending to apply shrinkage to entire groups of correlated predictors.
#' Negative values of lambda on the other hand favor sparsity within sets of correlated predictors by tending to apply shrinkage
#' to all but one predictor within a grouping. This can be understood as projecting the information matrix into a new space which
#' leads to a model similar in function to principal components regression (Krishna et al., 2009). Alternatively, this can be viewed
#' as a successor to ridge regression and an alternative to the elastic net. The model is estimated using an optimization routine.
#' Standard errors are computed from the hessian evaluated at the maximum a posteriori estimate.
#'
#'
#' \cr
#' Additional comments: This particular model is worth the monicker "author's pick". This model seems to
#' perform very well, and is relatively straightforward to estimate. Furthermore, the way it adapts to collinearity
#' seems more flexible than the elastic net. \cr \cr
#' I recommend only using numeric covariates here.
#'
#' @param formula a model formula
#' @param data a data frame
#' @param lambda a value for the power parameter. defaults to "auto" which means it is jointly estimated
#' with the other model parameters.
#' @param g a value for the g prior. if left as NULL, it defaults to nrow(X)/ncol(X)
#'
#' @return a roblm object
#' @export
#'
#' @references
#'
#' Zellner, A. & Siow S. (1980). Posterior odds ratio for selected regression hypotheses. In Bayesian statistics. Proc. 1st int. meeting (eds J. M. Bernardo, M. H. DeGroot, D. V. Lindley & A. F. M. Smith), 585–603. University Press, Valencia. \cr
#' \cr
#' Liang, Paulo, Molina, Clyde, & Berger (2008) Mixtures of g Priors for Bayesian Variable Selection, Journal of the American Statistical Association, 103:481, 410-423, DOI: 10.1198/016214507000001337 \cr
#' \cr
#' Krishna, A., Bondell, H. D., & Ghosh, S. K. (2009). Bayesian variable selection using an adaptive powered correlation prior. Journal of statistical planning and inference, 139(8), 2665–2674. doi:10.1016/j.jspi.2008.12.004
#' \cr
apclm <- function(formula, data, lambda="auto", g=NULL){

  if (is.null(g)){
    X <- model.matrix(formula, data)
    g<-nrow(X)/ncol(X) ; rm(X)
  }else{
    g<-g
  }
  if (g=="auto"){fixed.g<-F}else{fixed.g<-T}; if (lambda=="auto"){fixed.lambda<-F}else{fixed.lambda<-T}
  if (g=="auto"){
  if (lambda=="auto"){
    .logposterior <- function(coef,x,y){
      .apcprior <- function(X,lambda=-1){
        Trace = function(mat){sum(diag(mat))}
        xtx <- crossprod(X)
        if (lambda==0){xtxlam <- diag(1, ncol(X), ncol(X))}
        else{xtxlam <- cvreg:::ridgepow(xtx, lambda)}
        xtxinv <- cvreg:::ridgepow(xtx, -1)
        xtxlam <- mpd(symm((xtxlam / Trace(xtxlam)) * Trace(xtxinv)))
        xtxlam
      }
      x<-as.matrix(x); y<-as.vector(y)
      lambda<-coef[1]; sigma <- coef[2]; g <- coef[3]; coef <- coef[-c(1,2,3)]
      xtxlam <- .apcprior(X=x, lambda=lambda)
      yhat <-  as.vector(x%*%coef)
      zeros <- rep(0, length(coef))
      log_g      <- dgamma(1/g, shape = 1.5 , rate = 1.5 * (nrow(x)/(ncol(x))), log = T)
      log_lambda <- dshash(lambda, -2, 0.75, 2, 3,log=T) #dpowexp(lambda, 0, 2, 6, log = T)
      tau <- 1 / var(y); sh <- tau^2 / tau; ra <- tau / tau
      log_invsd2 <- dgamma(1/(sigma^2), sh, ra, log=T)
      log_intcpt <- dst(coef[1],3,mean(y),sd(y),log=T)
      log_coeffs <- sum(mvtnorm::dmvnorm(coef[-1], zeros[-1], sigma^2 * xtxlam[-1,-1] * g, log=T) )
      log_likehd <- sum(dnorm(y, yhat, sigma, log=T))
      log_post <- log_coeffs + log_likehd + log_invsd2 + log_intcpt + log_lambda + log_g
      return(-1 * log_post)
    }

    X = model.matrix(formula, data)
    Y = model.response(model.frame(formula, data))
    init <- as.vector(XtXinv(X)%*%crossprod(X,Y))
    init <- c(0, mad(as.vector(Y)-as.vector(X%*%init)), nrow(X)/ncol(X), init)
    fit = nlminb(init,
                 .logposterior,
                 x=X,
                 y=Y,
                 lower=c(-1500, 1e-9, 0, rep(-1e100, length(init)-1)),
                 upper=c(1500,  1e100, nrow(X), rep(1e100, length(init)-1)),
                 control=list(eval.max=10000,
                              iter.max=10000,
                              abs.tol=0,
                              sing.tol=1e-40)
    )
    gfun <- function(coef) {.logposterior(coef,X,Y) }
    hess <- optimx::gHgen(fit$par, gfun)$Hn
    std.err <- hessian2se(hess, vcov = T)
    vcov <- std.err$vcov
    std.err<-std.err$SE
    coef <- fit$par[-c(1,2,3)];
    lambda <- fit$par[1]; lambda.std.err <- std.err[1];
    sigma <- fit$par[2]; sigma.std.err <- std.err[2];
    g <- fit$par[3]; g.std.err <- std.err[3]; std.err <- std.err[-c(1,2,3)]
    fitted <- as.vector(X%*%coef); residuals <- sigma*((as.vector(Y)-fitted)/sd(as.vector(Y)-fitted))
    names(coef) <- names(std.err) <- colnames(X)
    res <- structure(list(coefficients = coef, std.err = std.err, sigma = sigma, lambda = lambda, g = g, g.se = g.std.err, sigma.se = sigma.std.err, lambda.se = lambda.std.err, fitted = fitted, residuals = residuals, log_post = fit$objective, vcov = vcov),class = "apclm")

  }
  else{

    .apcprior <- function(X,lambda=-1){
      Trace = function(mat){sum(diag(mat))}
      xtx <- crossprod(X)
      if (lambda==0){xtxlam <- diag(1, ncol(X), ncol(X))}
      else{xtxlam <- cvreg:::ridgepow(xtx, lambda)}
      xtxinv <- cvreg:::ridgepow(xtx, -1)
      xtxlam <- mpd(symm((xtxlam / Trace(xtxlam)) * Trace(xtxinv)))
      xtxlam
    }

    .logposterior <- function(coef,x,y,xtxlam){
      x<-as.matrix(x); y<-as.vector(y)
      sigma <- coef[1]; g <- coef[2];coef <- coef[-c(1,2)]
      yhat <-  as.vector(x%*%coef)
      zeros <- rep(0, length(coef))
      log_g      <- dgamma(1/g, shape = 1.5 , rate = 1.5 * ((nrow(x))/ncol(x)), log = T)
      tau <- 1 / var(y); sh <- tau^2 / tau; ra <- tau / tau
      log_invsd2 <- dgamma(1/(sigma^2), sh, ra, log=T)
      log_intcpt <- dst(coef[1],3,mean(y),sd(y),log=T)
      log_coeffs <- sum(mvtnorm::dmvnorm(coef[-1], zeros[-1], sigma^2 * xtxlam[-1,-1] * g,log=T))
      log_likehd <- sum(dnorm(y, yhat, sigma, log=T))
      log_post <- log_coeffs + log_likehd + log_invsd2 + log_intcpt + log_g
      return(-1 * log_post)
    }

    X = model.matrix(formula, data)
    Y = model.response(model.frame(formula, data))
    init <- as.vector(XtXinv(X)%*%crossprod(X,Y))
    init <- c(mad(as.vector(Y)-as.vector(X%*%init)), nrow(X)/ncol(X), init)
    xtxlam <- .apcprior(X, lambda)
    fit = nlminb(init,
                 .logposterior,
                 xtxlam=xtxlam,
                 x=X,
                 y=Y,
                 lower=c(1e-9, 0, rep(-1e100,length(init)-1)),
                 upper=c(1e100, nrow(X), rep(1e100, length(init)-1)),
                 control=list(eval.max=10000,
                              iter.max=10000,
                              abs.tol=0,
                              sing.tol=1e-40)
    )
    gfun <- function(coef) {.logposterior(coef,X,Y,xtxlam) }
    hess <- optimx::gHgen(fit$par, gfun)$Hn
    std.err <- hessian2se(hess, vcov = T)
    vcov <- std.err$vcov
    std.err<-std.err$SE
    coef <- fit$par[-c(1,2)]; sigma <- fit$par[1]; sigma.std.err <- std.err[1];
    g <- fit$par[2]; g.std.err <- std.err[2]; std.err <- std.err[-c(1,2)]
    names(coef) <- names(std.err) <- colnames(X)
    fitted <- as.vector(X%*%coef); residuals <- sigma*((as.vector(Y)-fitted)/sd(as.vector(Y)-fitted))
    names(coef) <- names(std.err) <- colnames(X)
    res <-structure(list(coefficients=coef,std.err=std.err,sigma=sigma,lambda=lambda,g=g,sigma.se=sigma.std.err,g.se=g.std.err,log_post=fit$objective,vcov=vcov),class = "apclm")

  }

}else{
  if (lambda=="auto"){

  .logposterior <- function(coef,x,y,G){
      .apcprior <- function(X,lambda=-1,g=NULL){
        Trace = function(mat){sum(diag(mat))}
        xtx <- crossprod(X)
        if (lambda==0){xtxlam <- diag(1, ncol(X), ncol(X))}
        else{xtxlam <- cvreg:::ridgepow(xtx, lambda)}
        xtxinv <- cvreg:::ridgepow(xtx, -1)
        xtxlam <- mpd(symm((xtxlam / Trace(xtxlam)) * Trace(xtxinv)))
        xtxlam * g
      }
      x<-as.matrix(x); y<-as.vector(y)
      lambda<-coef[1]; sigma <- coef[2]; coef <- coef[-c(1,2)]
      xtxlam <- .apcprior(X=x, lambda=lambda, g=G)
      yhat <-  as.vector(x%*%coef)
      zeros <- rep(0, length(coef))
      log_lambda <- dshash(lambda, -2, 0.75, 2, 3, log = T) #dpowexp(lambda, 0, 2, 6, log = T)
      tau <- 1 / var(y); sh <- tau^2 / tau; ra <- tau / tau
      log_invsd2 <- dgamma(1/(sigma^2), sh, ra, log=T)
      log_intcpt <- dst(coef[1],3,mean(y),sd(y),log=T)
      log_coeffs <- sum(mvtnorm::dmvnorm(coef[-1], zeros[-1], sigma^2 * xtxlam[-1,-1],log=T))
      log_likehd <- sum(dnorm(y, yhat, sigma, log=T))
      log_post <- log_coeffs + log_likehd + log_invsd2 + log_intcpt + log_lambda
      return(-1 * log_post)
    }

    X = model.matrix(formula, data)
    Y = model.response(model.frame(formula, data))
    if (is.null(g)){g<-nrow(X)/ncol(X)}else{g<-g}
    init <- as.vector(XtXinv(X)%*%crossprod(X,Y))
    init <- c(0, mad(as.vector(Y)-as.vector(X%*%init)), init)
    fit = nlminb(init,
               .logposterior,
               x=X,
               y=Y,
               G=g,
               lower=c(-1500, 1e-9, rep(-1e100, length(init)-1)),
               upper=c(1500,  1e100, rep(1e100, length(init)-1)),
               control=list(eval.max=10000,
                            iter.max=10000,
                            abs.tol=0,
                            sing.tol=1e-40)
    )
    gfun <- function(coef) {.logposterior(coef,X,Y,g) }
    hess <- optimx::gHgen(fit$par, gfun)$Hn
    std.err <- hessian2se(hess, vcov = T)
    vcov <- std.err$vcov
    std.err<-std.err$SE
    coef <- fit$par[-c(1,2)];
    lambda <- fit$par[1]; lambda.std.err <- std.err[1];
    sigma <- fit$par[2]; sigma.std.err <- std.err[2]; std.err <- std.err[-c(1,2)]
    fitted <- as.vector(X%*%coef); residuals <- sigma*((as.vector(Y)-fitted)/sd(as.vector(Y)-fitted))
    names(coef) <- names(std.err) <- colnames(X)
    res <-structure(list(coefficients = coef, std.err = std.err, sigma = sigma, lambda = lambda, g = g, sigma.se = sigma.std.err, lambda.se = lambda.std.err, g.se = NULL, fitted = fitted, residuals = residuals, log_post = fit$objective, vcov = vcov),class = "apclm")

  }
  else{

  .apcprior <- function(X,lambda=-1,g=NULL){
    if (is.null(g)){g<-nrow(X)/ncol(X)}else{g<-g}
    Trace = function(mat){sum(diag(mat))}
    xtx <- crossprod(X)
    if (lambda==0){xtxlam <- diag(1, ncol(X), ncol(X))}
    else{xtxlam <- cvreg:::ridgepow(xtx, lambda)}
    xtxinv <- cvreg:::ridgepow(xtx, -1)
    xtxlam <- mpd(symm((xtxlam / Trace(xtxlam)) * Trace(xtxinv)))
    xtxlam * g
  }

  .logposterior <- function(coef,x,y,xtxlam){
    x<-as.matrix(x); y<-as.vector(y)
    sigma <- coef[1]; coef <- coef[-1]
    yhat <-  as.vector(x%*%coef)
    zeros <- rep(0, length(coef))
    tau <- 1 / var(y); sh <- tau^2 / tau; ra <- tau / tau
    log_invsd2 <- dgamma(1/(sigma^2), sh, ra, log=T)
    log_intcpt <- dst(coef[1],3,mean(y),sd(y),log=T)
    log_coeffs <- sum(mvtnorm::dmvnorm(coef[-1], zeros[-1], sigma^2 * xtxlam[-1,-1],log=T))
    log_likehd <- sum(dnorm(y, yhat, sigma, log=T))
    log_post <- log_coeffs + log_likehd + log_invsd2 + log_intcpt
    return(-1 * log_post)
  }

  X = model.matrix(formula, data)
  Y = model.response(model.frame(formula, data))
  if (is.null(g)){g<-nrow(X)/ncol(X)}else{g<-g}
  init <- as.vector(XtXinv(X)%*%crossprod(X,Y))
  init <- c(mad(as.vector(Y)-as.vector(X%*%init)), init)
  xtxlam <- .apcprior(X, lambda, g)
  fit = nlminb(init,
             .logposterior,
             xtxlam=xtxlam,
             x=X,
             y=Y,
             lower=c(1e-9, rep(-1e100,length(init)-1)),
             upper=c(1e100, rep(1e100, length(init)-1)),
             control=list(eval.max=10000,
                          iter.max=10000,
                          abs.tol=0,
                          sing.tol=1e-40)
             )
  gfun <- function(coef) {.logposterior(coef,X,Y,xtxlam) }
  hess <- optimx::gHgen(fit$par, gfun)$Hn
  std.err <- hessian2se(hess, vcov = T)
  vcov <- std.err$vcov
  std.err<-std.err$SE
  coef <- fit$par[-1]; sigma <- fit$par[1]; sigma.std.err <- std.err[1]; std.err <- std.err[-1]
  names(coef) <- names(std.err) <- colnames(X)
  fitted <- as.vector(X%*%coef); residuals <- sigma*((as.vector(Y)-fitted)/sd(as.vector(Y)-fitted))
  names(coef) <- names(std.err) <- colnames(X)
  res <- structure(list(coefficients=coef,std.err=std.err,sigma=sigma,lambda=lambda,g=g,sigma.se=sigma.std.err,lambda.se=NULL,g.se=NULL,log_post=fit$objective,vcov=vcov),class = "apclm")
  }
}
  .apc <- function(X,lambda=-1,g=NULL){
    if (is.null(g)){g<-nrow(X)/ncol(X)}else{g<-g}
    Trace = function(mat){sum(diag(mat))}
    xtx <- crossprod(X)
    if (lambda==0){xtxlam <- diag(1, ncol(X), ncol(X))}
    else{xtxlam <- cvreg:::ridgepow(xtx, lambda)}
    xtxinv <- cvreg:::ridgepow(xtx, -1)
    xtxlam <- mpd(symm((xtxlam / Trace(xtxlam)) * Trace(xtxinv)))
    xtxlam * g
  }
  prior_covmat <- .apc(X, res$lambda, res$G) * (var(Y) * (nrow(X)-1) / nrow(X))
  prior_sigma <- sqrt(diag(prior_covmat))
  res$postOdds <-  exp(dnorm(coef, coef, std.err, T) - dnorm(0, coef, std.err, T))
  res$BF <- exp(dnorm(0, coef, std.err, T) - dnorm(0, 0, prior_sigma, T))
  res$sampler <- make_sampler(fixed.g, fixed.lambda, Y, X, res$coef, res$sigma, res$lambda, res$g, res$vcov, res$sigma.se, res$lambda.se, res$g.se)
  return(res)
}

#' print method for apclm objects
#'
#' @param fit the model fit
#' @param cred.level the credibility level for calculation of highest density posterior intervals. defaults to 0.95.
#' @param digits the number of significant digits to display. defaults to 3.
#' @method print apclm
#' @export
#' @keywords internal
#'
print.apclm <- function(fit, cred.level = 0.95, digits = 4){
  m <- as.matrix(fit$coefficients)
  m <- cbind(m, fit$std.err,
             lower.ci = m[,1] - (abs(qnorm((1 - cred.level)/2))*fit$std.err),
             upper.ci = m[,1] + (abs(qnorm((1 - cred.level)/2))*fit$std.err))
  #sig <- sapply(1:nrow(m), function(i) ifelse(sign(m[i,3]) == sign(m[i,4]), "*" , " "))
  rownames(m) <- names(fit$coefficients)
  m <- format(round(m, digits), nsmall = digits);
  m <- cbind.data.frame(m, formatNumber(fit$BF), formatNumber(fit$postOdds))
  m <- as.data.frame(m) #; m<-cbind.data.frame(m, sig)
  colnames(m) <- c("estimate", "std.err", paste0("lower ", 100*cred.level, "% CrI" ), paste0("upper ", 100*cred.level, "% CrI"), "BF0", "postOdds")
  print(noquote(m), right=T)

  calculate_ci <- function(ICDF, cred.level, tol=1e-8){
    intervalWidth <- function(lowTailPr, ICDF, cred.level) {
      ICDF(cred.level  + lowTailPr) - ICDF(lowTailPr)
    }
    optInfo <- optimize(intervalWidth, c(0, 1 - cred.level), ICDF = ICDF, cred.level  = cred.level, tol = tol)
    HDIlowTailPr <- optInfo$minimum
    result <- c(lower = ICDF(HDIlowTailPr), upper = ICDF(cred.level+HDIlowTailPr))
    attr(result, "cred.level ") <- cred.level
    return(result)
  }

  qgamma_factory<-function(mean,sd){
    function(p){
      q <- qgamma(p, shape = (mean^2)/(sd^2), rate = (mean)/(sd^2))
      return(q)
    }
  }

  qlambda_factory<-function(mean,sd){
    function(p){q <- qnorm(p, mean, sd); return(q)}
  }

  cifun <- qlambda_factory(fit$sigma, fit$sigma.se)
  sigma_ci <- calculate_ci(cifun, cred.level)
  pars <- cbind(fit$sigma, fit$sigma.se, t(sigma_ci))
  colnames(pars) <- c("estimate", "std.err", paste0("lower ", 100*cred.level, "% CrI" ), paste0("upper ", 100*cred.level, "% CrI"))
  rownames(pars) <- "sigma"
  if (!is.null(fit$lambda.se)){
    cifun <- qlambda_factory(fit$lambda, fit$lambda.se)
    lambda_ci <- calculate_ci(cifun, cred.level)
    pars <- rbind(pars, cbind(fit$lambda, fit$lambda.se, t(lambda_ci)))
    rownames(pars) <- c("sigma", "lambda")
  }
  current.rownames<-rownames(pars)
  if (!is.null(fit$g.se)){
    cifun <- qgamma_factory(fit$g, fit$g.se)
    g_ci <- calculate_ci(cifun, cred.level)
    pars <- rbind(pars, cbind(fit$g, fit$g.se, t(g_ci)))
    rownames(pars) <- c(current.rownames, "g")
  }
  cat("\n", crayon::underline("model parameter estimates:\n\n"))
  pars <- format(round(pars, digits), nsmall = digits)
  print(noquote(pars), right=T)
}


#' Predict method for apclm models
#'
#'
#' @param fit the model fit
#' @param newdata a data frame containing the new data. can be left NULL, in which case the fitted
#' observations are returned.
#'
#' @return a vector of predictions
#' @method predict apclm
#' @export
predict.apclm <- function(fit, newdata = NULL){
  if (is.null(newdata)){
    return(fit$fitted)
  }
  else{
    if (!is.data.frame(newdata)) {newdata <- as.data.frame(newdata)}
    terms <- colnames(fit$model.frame)
    newdata <- newdata[,terms]
    if (any(names(fit$coeffcients) == "(Intercept)")) {
      x <- model.matrix(~., newdata)
    }
    else{
      x <- model.matrix(~0+., newdata)
    }
    return(as.vector(x%*%fit$coefficients))
  }
}





make_sampler <- function(fixed.g, fixed.lambda, y, X, coef, sigma, lambda, g, vcov, sigma.se, lambda.se=NULL, g.se=NULL){

  if (!fixed.g && !fixed.lambda){
    return(function(reps=1000) {
      list <- list()
      pvar <- length(coef)
      list[["prior"]] <- .sample_prior1(y, X, reps)
      list[["posterior"]] <- .sample_posterior1(X, coef, sigma, lambda, g, vcov[1:pvar, 1:pvar], sigma.se, lambda.se, g.se, reps)
      list
    })
  }
  else if (!fixed.g && fixed.lambda){
    return(function(reps=1000) {
      list <- list()
      pvar <- length(coef)
      list[["prior"]] <- .sample_prior2(y, X, lambda, reps)
      list[["posterior"]] <- .sample_posterior2(X, coef, sigma, lambda, g, vcov[1:pvar, 1:pvar], sigma.se, g.se, reps)
      list
    })
  }
  else if(fixed.g && !fixed.lambda){
    return(function(reps=1000) {
      list <- list()
      pvar <- length(coef)
      list[["prior"]] <- .sample_prior3(y, X, g, reps)
      list[["posterior"]] <- .sample_posterior3(X, coef, sigma, lambda, g, vcov[1:pvar, 1:pvar], sigma.se, lambda.se, reps)
      list
    })
  }
  else if(fixed.g && fixed.lambda){
    return(function(reps=1000) {
      list <- list()
      pvar <- length(coef)
      list[["prior"]] <- .sample_prior4(y, X, lambda, g, reps)
      list[["posterior"]] <- .sample_posterior4(X, coef, sigma, lambda, g, vcov[1:pvar, 1:pvar], sigma.se, reps)
      list
    })
  }

}



.apc <- function(X,lambda=-1,g=NULL){
  if (is.null(g)){g<-nrow(X)/ncol(X)}else{g<-g}
  Trace = function(mat){sum(diag(mat))}
  xtx <- crossprod(X)
  if (lambda==0){xtxlam <- diag(1, ncol(X), ncol(X))}
  else{xtxlam <- cvreg:::ridgepow(xtx, lambda)}
  xtxinv <- cvreg:::ridgepow(xtx, -1)
  xtxlam <- mpd(symm((xtxlam / Trace(xtxlam)) * Trace(xtxinv)))
  xtxlam * g
}

## Varying lambda, varying g

.sample_prior1 <- function(y, X, reps = 1000){
  n = nrow(X); p = ncol(X)-1; ymu <- mean(y); tau <- 1 / var(y)
  sampler <- function(rep){
    rep <- 1
    g  <- 1 / rgamma(1, shape = 1.5, rate = 1.5 * (n/p))
    lambda <- rshash(1, -2, 0.75, 2, 3)
    sigma <- sqrt(1 / rgamma(1, tau^2 / tau, tau / tau))
    intcpt <- rnorm(1, ymu, 1)
    xtxlam <- .apc(X, lambda, g)
    coeffs <- as.vector(mvtnorm::rmvnorm(1, rep(0, p), sigma^2 * xtxlam[-1,-1] * g))
    ppd <- rnorm(1, sum(X %*% c(intcpt, coeffs)), sigma)
    return(c(Intercept = intcpt, as.vector(coeffs), sigma = sigma, lambda = lambda, g = g, ppd = ppd))
  }
  out <- t(sapply(1:reps, sampler))
  colnames(out) <- c(colnames(X), "sigma", "lambda", "g", "ppd")
  return(out)
}


.sample_posterior1 <- function(X, coef, sigma, lambda, g, vcov, sigma.se, lambda.se, g.se, reps = 1000){
  sampler <- function(rep){
    rep <- 1
    g  <- rgamma(1, shape = g^2 / g.se^2, rate = g / g.se^2)
    lambda <- rnorm(1, lambda, lambda.se)
    sigma <- rgamma(1, sigma^2 / sigma.se^2, sigma / sigma.se^2)
    coeffs <- as.vector(MASS::mvrnorm(1, coef, mpd(vcov), empirical = F))
    mu <- mean(X %*% coeffs)
    yhat <- list()
    for (i in 1:nrow(X)){
      yhat[[i]] <- rnorm(1, X[i,]%*%coeffs, sigma)
    }
    return(c(as.vector(coeffs), sigma = sigma, lambda = lambda, g = g, mu = mu, ypred = as.vector(unlist(yhat))))
  }
  out <- t(sapply(1:reps, sampler))
  colnames(out) <- c(colnames(X), "sigma", "lambda", "g", "mu", paste0("ypred", 1:nrow(X)))
  return(out)
}


## Fixed lambda , Varying g

.sample_prior2 <- function(y, X, lambda, reps = 1000){
  n = nrow(X); p = ncol(X)-1; ymu <- mean(y); tau <- 1 / var(y)
  sampler <- function(rep){
    rep <- 1
    g  <- 1 / rgamma(1, shape = 1.5, rate = 1.5 * (n/p))
    sigma <- sqrt(1 / rgamma(1, tau^2 / tau, tau / tau))
    intcpt <- rnorm(1, ymu, 1)
    xtxlam <- .apc(X, lambda, g)
    coeffs <- as.vector(mvtnorm::rmvnorm(1, rep(0, p), sigma^2 * xtxlam[-1,-1] * g))
    ppd <- rnorm(1, sum(X %*% c(intcpt, coeffs)), sigma)
    return(c(Intercept = intcpt, as.vector(coeffs), sigma = sigma, g = g, ppd = ppd))
  }
  out <- t(sapply(1:reps, sampler))
  colnames(out) <- c(colnames(X), "sigma", "g", "ppd")
  return(out)
}


.sample_posterior2 <- function(X, coef, sigma, lambda, g, vcov, sigma.se, g.se, reps = 1000){
  sampler <- function(rep){
    rep <- 1
    g  <- rgamma(1, shape = g^2 / g.se^2, rate = g / g.se^2)
    sigma <- rgamma(1, sigma^2 / sigma.se^2, sigma / sigma.se^2)
    coeffs <- as.vector(MASS::mvrnorm(1, coef, mpd(vcov), empirical = F))
    yhat <- list()
    for (i in 1:nrow(X)){
      yhat[[i]] <- rnorm(1, X[i,]%*%coeffs, sigma)
    }
    return(c(as.vector(coeffs), sigma = sigma, g = g, mu = mu, ypred = as.vector(unlist(yhat))))
  }
  out <- t(sapply(1:reps, sampler))
  colnames(out) <- c(colnames(X), "sigma", "g", paste0("ypred", 1:nrow(X)))
  return(out)
}



## Fixed g, varying lambda

.sample_prior3 <- function(y, X, g, reps = 1000){
  n = nrow(X); p = ncol(X)-1; ymu <- mean(y); tau <- 1 / var(y)
  sampler <- function(rep){
    rep <- 1
    lambda <- rshash(1, -2, 0.75, 2, 3)
    sigma <- sqrt(1 / rgamma(1, tau^2 / tau, tau / tau))
    intcpt <- rnorm(1, ymu, 1)
    xtxlam <- .apc(X, lambda, g)
    coeffs <- as.vector(mvtnorm::rmvnorm(1, rep(0, p), sigma^2 * xtxlam[-1,-1] * g))
    ppd <- rnorm(1, sum(X %*% c(intcpt, coeffs)), sigma)
    return(c(Intercept = intcpt, as.vector(coeffs), sigma = sigma, lambda = lambda, ppd = ppd))
  }
  out <- t(sapply(1:reps, sampler))
  colnames(out) <- c(colnames(X), "sigma", "lambda", "ppd")
  return(out)
}


.sample_posterior3 <- function(X, coef, sigma, lambda, g, vcov, sigma.se, lambda.se, reps = 1000){
  sampler <- function(rep){
    rep <- 1
    lambda <- rnorm(1, lambda, lambda.se)
    sigma <- rgamma(1, sigma^2 / sigma.se^2, sigma / sigma.se^2)
    coeffs <- as.vector(MASS::mvrnorm(1, coef[-1], mpd(vcov[-1,-1]), empirical = F, tol = 1e-100))
    coeffs <- c(rnorm(1, coef[1], sqrt(diag(vcov))[1]), coeffs)
    yhat <- list()
    for (i in 1:nrow(X)){
      yhat[[i]] <- rnorm(1, X[i,]%*%coeffs, sigma)
    }
    return(c(as.vector(coeffs), sigma = sigma, lambda = lambda, ypred = as.vector(unlist(yhat))))
  }
  out <- t(sapply(1:reps, sampler))
  colnames(out) <- c(colnames(X), "sigma", "lambda", paste0("ypred", 1:nrow(X)))
  return(out)
}


## Fixed lambda, fixed g

.sample_prior4 <- function(y, X, lambda, g, reps = 1000){
  n = nrow(X); p = ncol(X)-1; ymu <- mean(y); tau <- 1 / var(y)
  sampler <- function(rep){
    rep <- 1
    sigma <- sqrt(1 / rgamma(1, tau^2 / tau, tau / tau))
    intcpt <- rnorm(1, ymu, 1)
    xtxlam <- .apc(X, lambda, g)
    coeffs <- as.vector(mvtnorm::rmvnorm(1, rep(0, p), sigma^2 * xtxlam[-1,-1] * g))
    ppd <- rnorm(1, sum(X %*% c(intcpt, coeffs)), sigma)
    return(c(Intercept = intcpt, as.vector(coeffs), sigma = sigma, ppd = ppd))
  }
  out <- t(sapply(1:reps, sampler))
  colnames(out) <- c(colnames(X), "sigma", "ppd")
  return(out)
}


.sample_posterior4 <- function(X, coef, sigma, lambda, g, vcov, sigma.se, reps = 1000){
  sampler <- function(rep){
    rep <- 1
    sigma <- rgamma(1, sigma^2 / sigma.se^2, sigma / sigma.se^2)
    coeffs <- as.vector(MASS::mvrnorm(1, coef, mpd(vcov), empirical = F))
    yhat <- list()
    for (i in 1:nrow(X)){
      yhat[[i]] <- rnorm(1, X[i,]%*%coeffs, sigma)
    }
    return(c(as.vector(coeffs), sigma = sigma, ypred = as.vector(unlist(yhat))))
  }
  out <- t(sapply(1:reps, sampler))
  colnames(out) <- c(colnames(X), "sigma", paste0("ypred", 1:nrow(X)))
  return(out)
}