bayesQR (version 2.0)

QRc: Bayesian quantile regression with continuous dependent variable

Description

QRc implements a Bayesian method for estimating quantile regression parameters. This method uses the asymmetric Laplace distribution for this goal (see references). To improve the speed of the routine, the Markov Chain Monte Carlo (MCMC) part of the algorithm is programmed in Fortran and is called from within the R function QRc.

Usage

QRc(Data, Prior, Mcmc)

Arguments

Data

list(y, X, p)

Prior

list(beta0, V0, shape0, scale0)

Mcmc

list(R, keep)

Value

• A list containing:
• method
  a string containing 'QRc'
• p
  the quantile that was estimated
• betadraw
  R/keep x nvar(X) matrix of beta draws
• sigmadraw
  R/keep vector of sigma draws

Details

ll{ Model: y = Xbeta + e e ~ ALD(location=0, scale=sigma, quantile=p) Priors: beta ~ N(beta0, V0) sigma ~ invGamma(shape0,scale0) } List arguments contain:

• X: n x nvar(X) matrix of predictors (first column should be a vector of ones if intercept is desired)
• y: n x 1 vector of observations (dependent variable)
• p: quantile of interest (should be between 0 and 1)
• beta0: nvar(X) x 1 vector of prior means (default: 0)
• V0: nvar(X) x nvar(X) prior covariance matrix (default: 100*diag(ncol(X)))
• shape0: shape parameter for inverse Gamma prior for sigma (default: 0.01)
• scale0: scale parameter for inverse Gamma prior for sigma (default: 0.01)
• R: number of MCMC draws
• keep: thinning parameter, i.e. keep every keepth draw (default: 1)

References

The algorithm is an improved (Gibbs instead of Metropolis-Hastings) implementation of: Yu, K. and Moyeed, R. (2001). Bayesian quantile regression, Statistics and Probability Letters, 54, 437-447. Also see: Kozumi, H. and Kobayashi, G. (2011). Gibbs sampling methods for Bayesian quantile regression, Journal of Statistical Computation and Simulation, 81(11), 1565-1578.

Examples

# Simulate data from heteroskedastic regression
n <- 200
X <- runif(n=n,min=0,max=10)
X <- cbind(1,X)
y <- 1 + 2*X[,2] + rnorm(n=n, mean=0, sd=.6*X[,2])

# Initiate plot
## Plot datapoints
plot(X[,2], y, main="", cex=.6, xlab="X")

# Initialize the inputs for QRc
Data = list(y=y, X=X, p=.05)
Prior = list(beta0=rep(0,ncol(X)),V0=100*diag(ncol(X)), shape0=0.01, scale0=0.01)
Mcmc = list(R=5000, keep=1)

# Write loop to analyze 5 quantiles
for (i in 1:5) {
    if (i==1) p = .05
    if (i==2) p = .25
    if (i==3) p = .50
    if (i==4) p = .75
    if (i==5) p = .95

    Data$p = p
    out = QRc(Data=Data, Prior=Prior, Mcmc=Mcmc)

    ## Add quantile regression lines to the plot (exclude first 500 burn-in draws)
    abline(a=mean(out$betadraw[500:Mcmc$R,1]),b=mean(out$betadraw[500:Mcmc$R,2]),lty=i,col=i)
}

# Estimate and plot OLS model
outOLS = lm(y~X[,2])
abline(outOLS,lty=1,lwd=2,col=6)

# Add legend to plot
legend(x=0,y=max(y),legend=c(.05,.25,.50,.75,.95,"OLS"),lty=c(1,2,3,4,5,1),lwd=c(1,1,1,1,1,2),col=c(1:6),title="Quantile")