R/bridgeDC.R

In abnormally-distributed/Bayezilla: Bayesian Models With JAGS and Several Utilities for Data Exploration and Model Evaluation.

#' Bayesian Bridge Regression with unpenalized design covariates
#'
#' @description The Bayesian Bridge model of Mallick & Yi (2018), but with the allowance for a set of covariates that are not penalized. 
#' For example, you may wish to include variables such as age and gender in all models so that 
#' the coefficients for the other variables are penalized while controlling for these. This
#' is a common need in research. Bridge regression allows you to utilize different Lp norms for the shape 
#' of the prior through the shape parameter kappa of the power exponential distribution (also known as generalized Gaussian). 
#' Norms of 1 and 2 give the Laplace and Gaussian distributions respectively (corresponding to the LASSO and Ridge Regression). 
#' Norms smaller than 1 are very difficult to estimate directly, but have very tall modes at zero and very long, cauchy like tails. 
#' Values greater than 2 become increasingly platykurtic, with the uniform distribution arising as it approaches infinity. #' Using kappa = 1 yields the New Bayesian LASSO, which is a re-parameterization of the Bayesian LASSO utilizing a scale mixture of
#' uniform distributions to obtain the Laplacian priors (Mallick & Yi, 2014).  \cr
#' \cr 
#' \cr
#' JAGS has no built in power exponential distribution, so the distribution is parameterized as a uniform-gamma mixture just as in Mallick & Yi (2018). 
#' The parameterization is given below. For generalized linear models plug-in pseudovariances are used. \cr
#' \cr
#' Model Specification:
#' \cr
#' \if{html}{\figure{bridgeDC.png}{}}
#' \if{latex}{\figure{bridgeDC.png}{}}
#' \cr
#' \cr
#' Plugin Pseudo-Variances: \cr
#' \if{html}{\figure{pseudovar.png}{}}
#' \if{latex}{\figure{pseudovar.png}{}}
#'
#'
#' @param formula the model formula
#' @param design.formula formula for the design covariates.
#' @param data a data frame.
#' @param family one of "gaussian", "binomial", or "poisson".
#' @param kappa the Lp norm you wish to utilize. Default is 1.4.
#' @param log_lik Should the log likelihood be monitored? The default is FALSE.
#' @param iter How many post-warmup samples? Defaults to 10000.
#' @param warmup How many warmup samples? Defaults to 1000.
#' @param adapt How many adaptation steps? Defaults to 2000.
#' @param chains How many chains? Defaults to 4.
#' @param thin Thinning interval. Defaults to 1.
#' @param method Defaults to "parallel". For an alternative parallel option, choose "rjparallel" or. Otherwise, "rjags" (single core run).
#' @param cl Use parallel::makeCluster(# clusters) to specify clusters for the parallel methods. Defaults to two cores.
#' @param ... Other arguments to run.jags.
#'
#' @references Mallick, H. & Yi, N. (2018) Bayesian bridge regression, Journal of Applied Statistics, 45:6, 988-1008, DOI: 10.1080/02664763.2017.1324565 \cr
#' \cr
#' Mallick, H., & Yi, N. (2014). A New Bayesian Lasso. Statistics and its interface, 7(4), 571–582. doi:10.4310/SII.2014.v7.n4.a12 \cr
#' 
#' @return
#' a runjags object
#' @export
#'
#' @examples
#' bridgeDC()
#' 
bridgeDC = function(formula, design.formula, data,  family = "gaussian", kappa = 1.4, log_lik = FALSE, iter=10000, warmup=1000, adapt=2000, chains=4, thin=1, method = "parallel", cl = makeCluster(2), ...){
  
  X = as.matrix(model.matrix(formula, data)[,-1])
  y = model.frame(formula, data)[,1]
  FX <- as.matrix(model.matrix(design.formula, data)[, -1])
  
  if (family == "gaussian"){
    
    jags_bridge = "model{
  
  tau ~ dgamma(.01, .01) 
  sigma <- sqrt(1/tau)
  lambda ~ dgamma(0.50 , 0.20)

  for (i in 1:P){
    u[i] ~ dgamma( (1/kappa) + 1  , lambda )
    beta[i] ~ dunif(-1 * pow(sigma * u[i], 1/kappa), pow(sigma * u[i], 1/kappa))
  }
  
  for (f in 1:FP){
    design_beta[f] ~ dnorm(0, 1e-200)
  }
  
  Intercept ~ dnorm(0, 1e-10)
  
  for (i in 1:N){
    y[i] ~ dnorm(Intercept + sum(beta[1:P] * X[i,1:P]) + sum(design_beta[1:FP] * FX[i,1:FP]) , tau)
    log_lik[i] <- logdensity.norm(y[i], Intercept + sum(beta[1:P] * X[i,1:P]) + sum(design_beta[1:FP] * FX[i,1:FP]) , tau)
    ySim[i] ~ dnorm(Intercept + sum(beta[1:P] * X[i,1:P]) + sum(design_beta[1:FP] * FX[i,1:FP]) , tau)
  }
  
  Deviance <- -2 * sum(log_lik[1:N])
}"
    FP <- ncol(FX)
    P <- ncol(X)
    write_lines(jags_bridge, "jags_bridge.txt")
    jagsdata <- list(X = X, y = y,  FP = FP, FX = FX, N = length(y), P = ncol(X), kappa = kappa)
    monitor <- c("Intercept", "beta", "sigma", "lambda", "Deviance", "ySim", "log_lik")
    if (log_lik == FALSE){
      monitor = monitor[-(length(monitor))]
    }
    inits <- lapply(1:chains, function(z) list("Intercept" = lmSolve(formula, data)[1], 
                                               "beta" = rep(0, P), 
                                               "design_beta" =  lmSolve(design.formula, data)[-1], 
                                               "u" = rgamma(P, (1 / kappa) + 1, 1), 
                                               "lambda" = 1, 
                                               "tau" = 1, 
                                               "ySim" = sample(y, length(y)),
                                               .RNG.name= "lecuyer::RngStream", 
                                               .RNG.seed= sample(1:10000, 1)))
    
    out = run.jags(model = "jags_bridge.txt", modules = c("bugs on", "glm on", "dic off"), monitor = monitor, data = jagsdata, inits = inits, burnin = warmup, sample = iter, thin = thin, adapt = adapt, method = method, cl = cl, summarise = FALSE, ...)
    file.remove("jags_bridge.txt")
    if (!is.null(cl)) {
      parallel::stopCluster(cl = cl)
    }
    return(out)
  }
  
  if (family == "binomial" || family == "logistic"){
    
    jags_bridge = "model{


  lambda ~ dgamma(0.50 , 0.20)
  
  for (i in 1:P){
    u[i] ~ dgamma( (1/kappa) + 1  , lambda )
    beta[i] ~ dunif(-1 * pow(sigma * u[i], 1/kappa), pow(sigma * u[i], 1/kappa))
  }
  
  for (f in 1:FP){
    design_beta[f] ~ dnorm(0, 1e-200)
  }
  
  Intercept ~ dnorm(0, 1e-10)
  
  for (i in 1:N){
      logit(psi[i]) <- Intercept + sum(beta[1:P] * X[i,1:P]) + sum(design_beta[1:FP] * FX[i,1:FP]) 
      y[i] ~ dbern(psi[i])
      log_lik[i] <- logdensity.bern(y[i], psi[i])
      ySim[i] ~ dbern(psi[i])
    }
  
  Deviance <- -2 * sum(log_lik[1:N])
}"
    FP <- ncol(FX)
    P <- ncol(X)
    write_lines(jags_bridge, "jags_bridge.txt")
    jagsdata <- list(X = X, y = y, FP = FP, FX = FX,  N = length(y), P = ncol(X), kappa = kappa, sigma = sqrt(pow(mean(y), -1) * pow(1 - mean(y), -1)))
    monitor <- c("Intercept", "beta", "kappa", "lambda", "Deviance", "ySim", "log_lik")
    if (log_lik == FALSE){
      monitor = monitor[-(length(monitor))]
    }
    inits <- lapply(1:chains, function(z) list("Intercept" = as.vector(coef(glmnet::glmnet(x = X, y = y, family = "binomial", lambda = runif(1, .01, .15), alpha = runif(1, .2, .5), standardize = FALSE))[1,1]), 
                                               "beta" = rep(0, P), 
                                               "design_beta" =  as.vector(coef(glmnet::glmnet(x = FX, y = y, family = "binomial", lambda = runif(1, .01, .15), alpha = runif(1, .2, .5), standardize = FALSE))[-1,1]), 
                                               "u" = rgamma(P, (1 / kappa) + 1, 1), 
                                               "lambda" = 1, 
                                               "ySim" = sample(y, length(y)),
                                               .RNG.name= "lecuyer::RngStream", 
                                               .RNG.seed= sample(1:10000, 1)))
    
    out = run.jags(model = "jags_bridge.txt", modules = c("bugs on", "glm on", "dic off"), monitor = monitor, data = jagsdata, inits = inits, burnin = warmup, sample = iter, thin = thin, adapt = adapt, method = method, cl = cl, summarise = FALSE, ...)
    file.remove("jags_bridge.txt")
    if (!is.null(cl)) {
      parallel::stopCluster(cl = cl)
    }
    return(out)
  }
  
  
  if (family == "poisson"){
    
    jags_bridge = "model{


  lambda ~ dgamma(0.50 , 0.20)
  
    
  for (f in 1:FP){
    design_beta[f] ~ dnorm(0, 1e-200)
  }
  
  
  for (i in 1:P){
    u[i] ~ dgamma( (1/kappa) + 1  , lambda )
    beta[i] ~ dunif(-1 * pow(sigma * u[i], 1/kappa), pow(sigma * u[i], 1/kappa))
  }
  
  Intercept ~ dnorm(0, 1e-10)
  
  for (i in 1:N){
    log(psi[i]) <- Intercept + sum(beta[1:P] * X[i,1:P]) + sum(design_beta[1:FP] * FX[i,1:FP]) 
    y[i] ~ dpois(psi[i])
    log_lik[i] <- logdensity.pois(y[i], psi[i])
    ySim[i] ~ dpois(psi[i])
}
              
  Deviance <- -2 * sum(log_lik[1:N])
}"
  FP <- ncol(FX)
  P <- ncol(X)
  write_lines(jags_bridge, "jags_bridge.txt")
  jagsdata <- list(X = X, y = y, FP = FP, FX = FX, N = length(y), P = ncol(X), kappa = kappa, sigma = sqrt(pow(mean(y) , -1)))
  monitor <- c("Intercept", "beta", "kappa", "lambda", "Deviance", "ySim", "log_lik")
  if (log_lik == FALSE){
    monitor = monitor[-(length(monitor))]
  }
  inits <- lapply(1:chains, function(z) list("Intercept" = as.vector(coef(glmnet::glmnet(x = X, y = y, family = "poisson", lambda = runif(1, .01, .15), alpha = runif(1, .2, .5), standardize = FALSE))[1,1]), 
                                             "beta" = rep(0, P), 
                                             "design_beta" =  as.vector(coef(glmnet::glmnet(x = FX, y = y, family = "poisson", lambda = runif(1, .01, .15), alpha = runif(1, .2, .5), standardize = FALSE))[-1,1]), 
                                             "u" = rgamma(P, (1 / kappa) + 1, 1), 
                                             "lambda" = 1, 
                                             "ySim" = sample(y, length(y)),
                                             .RNG.name= "lecuyer::RngStream", 
                                             .RNG.seed= sample(1:10000, 1)))
  
  out = run.jags(model = "jags_bridge.txt", modules = c("bugs on", "glm on", "dic off"), monitor = monitor, data = jagsdata, inits = inits, burnin = warmup, sample = iter, thin = thin, adapt = adapt, method = method, cl = cl, summarise = FALSE, ...)
  file.remove("jags_bridge.txt")
  if (!is.null(cl)) {
    parallel::stopCluster(cl = cl)
  }
  return(out)
  }

}