BSquare (version 1.1)

qreg: Quantile regression with parametric basis functions.)

Description

MCMC code for the quantile regression model of Reich and Smith, 2013.

Usage

qreg(X,Y=NULL,Y_low=NULL,Y_high=NULL,status=NULL, L=4,base="Gaussian",varying_effect=NULL, tau=seq(0.05,0.95,0.05), burn=10000,iters=50000)

Arguments

X

Matrix of predictors with the first column consisting of all ones and all other values between -1 and 1.

Y

A vector of responses.

Y_low,Y_high

Vectors of endpoints for interval-censored values.

status

Censoring status taking values 0 if uncensored 1 if left-censored on the interval (-Inf,Y) 2 if right-censored on the interval (Y,Inf) 3 if censored on the interval (Y_low,Y_high).

L

The number of basis functions in quantile function

base

The centering distribution which can take values "Gaussian", "t", "logistic", "gamma", "weibull", or "ALAP."

varying_effect

If varying_effect = j, then only the covariates in the first j columns of X have different effects on different quantile levels.

tau

Vector of quantile levels for output.

burn

Number of MCMC samples to discard as burn-in.

iters

Number of MCMC samples to generate after the burn-in.

Value

q

Posterior samples of the quantile function.

LPML

Log pseudo-maximum likelihood statistic for model comparisons.

Details

See http://www4.stat.ncsu.edu/~reich/QR/ for more detailed descriptions and examples.

References

Reich BJ, Smith LB (2013). Bayesian quantile regression for censored data. In press, Biometrics. Smith LB, Fuentes M, Herring AH, Reich BJ (2013) Bayesian dependent quantile regression processes for birth outcomes. Submitted. Reich BJ (2012) Spatiotemporal quantile regression for detecting distributional changes in environmental processes. JRSS-C, 64, 535-553. Reich BJ, Fuentes M, Dunson DB (2011) Bayesian spatial quantile regression. JASA, 106, 620.

See Also

dqreg qr_plot

Examples

  #Continuous data example
  #Load the air quality data
  data(airquality)
  ozone<-airquality[,1]
  solar<-airquality[,2]
  
  #Remove missing observations
  missing<-is.na(ozone+solar)
  ozone<-ozone[!missing]
  solar<-solar[!missing]
  solar_std<-1.8*(solar - min(solar))/(max(solar)-min(solar)) - 0.9
  
  #Fit the model and plot results
  X<-cbind(1,solar_std)
  #use longer chains in practice
  fit<-qreg(X,ozone,L=4,base="gamma", iters = 1000, burn = 1000)
  qr_plot(fit,2, main = "Solar Effect")

  
  #Right-censored data example
  
  library(survival)
  data(veteran)
  
  trt<-ifelse(veteran[,1]==2,-1,1)
  logtime<-log(veteran[,3])
  event<-veteran[,4]
  status<-ifelse(event==1,0,2)
  X<-cbind(1,trt)
  #use longer chains in practice
  fit<-qreg(X,Y=logtime,status=status,iters =1000, burn = 1000)
  qr_plot(fit,index=2,main="Treatment effect")