Bayesian Penalized Quantile Regression
The quantile varying coefficient model is robust to data heterogeneity, outliers and heavy-tailed distributions in the response variable. In addition, it can flexibly model dynamic patterns of regression coefficients through nonparametric varying coefficient functions. In this package, we have implemented the Gibbs samplers of the penalized Bayesian quantile varying coefficient model with spike-and-slab priors [Zhou et al.(2023)]doi:10.1016/j.csda.2023.107808 for efficient Bayesian shrinkage estimation, variable selection and statistical inference. In particular, valid Bayesian inferences on sparse quantile varying coefficient functions can be validated on finite samples. The Markov Chain Monte Carlo (MCMC) algorithms of the proposed and alternative models can be efficiently performed by using the package.
install.packages("devtools")
devtools::install_github("cenwu/pqrBayes")
install.packages("pqrBayes")
Data <- function(n,p,quant){
sig1 = matrix(0,p,p)
diag(sig1)=1
for (i in 1: p)
{
for (j in 1: p)
{
sig1[i,j]=0.5^abs(i-j)
}
}
x = mvrnorm(n,rep(0,p),sig1)
x = cbind(1,x)
error=rt(n,2) -quantile(rt(n,2),probs = quant)
u = runif(n,0.01,0.99)
gamma0 = 2+2*sin(u*2*pi)
gamma2 = -6*u*(1-u)
gamma1 = 2*exp(2*u-1)
gamma3= -4*u^3
y = gamma1*x[,2] + gamma2*x[,3] + gamma3*x[,4] + gamma0 + error
dat = list(y=y, u=u, x=x, gamma=cbind(gamma0,gamma1,gamma2,gamma3))
return(dat)
}
n=250; p=100; rep=200;
quant = 0.5; # focus on median for Bayesian inference
CI_RBGLSS = CI_RBGL = CI_BGLSS = CI_BGL= c()
for (h in 1:rep) {
dat = Data(n,p,quant)
y = dat$y
u = dat$u
x = dat$x
g = x[,-1]
kn=2
degree=2
u.grid = (1:n)*0.005
gamma_0_grid = 2+2*sin(2*u.grid*pi)
gamma_1_grid = 2*exp(2*u.grid-1)
gamma_2_grid = -6*u.grid*(1-u.grid)
gamma_3_grid = -4*u.grid^3
coefficient = cbind(gamma_0_grid,gamma_1_grid,gamma_2_grid,gamma_3_grid)
fit = pqrBayes(g, y, u, e=NULL,quant=quant, iterations=10000, kn=2, degree=2, robust = TRUE, sparse=TRUE, hyper=NULL,debugging=FALSE)
posterior=fit$posterior
coverage = coverage(posterior,coefficient,u.grid,kn,degree)
fit1 = pqrBayes(g, y, u, e=NULL,quant=quant, iterations=10000, kn=2, degree=2, robust = TRUE, sparse=FALSE, hyper=NULL,debugging=FALSE)
posterior1=fit1$posterior
coverage1 = coverage(posterior1,coefficient,u.grid,kn,degree)
fit2 = pqrBayes(g, y, u, e=NULL,quant=quant, iterations=10000, kn=2, degree=2, robust = FALSE, sparse=TRUE, hyper=NULL,debugging=FALSE)
posterior2=fit2$posterior
coverage2 = coverage(posterior2,coefficient,u.grid,kn,degree)
fit3 = pqrBayes(g, y, u, e=NULL,quant=quant, iterations=10000, kn=2, degree=2, robust = FALSE, sparse=FALSE, hyper=NULL,debugging=FALSE)
posterior3=fit3$posterior
coverage3 = coverage(posterior3,coefficient,u.grid,kn,degree)
CI_RBGLSS = rbind(CI_RBGLSS,coverage)
CI_RBGL = rbind(CI_RBGL,coverage1)
CI_BGLSS = rbind(CI_BGLSS,coverage2)
CI_BGL = rbind(CI_BGL,coverage3)
cat("iteration = ", h, "\n")
}
# the intercept gamma_0 has not been regularized
cp_RBGLSS = colMeans(CI_RBGLSS) # 95% coverage probabilities for the varying coeffcients under the default setting
cp_BGLSS = colMeans(CI_BGLSS)
cp_RBGL = colMeans(CI_RBGL)
cp_BGL = colMeans(CI_BGL)
This package provides implementation for methods proposed in