R/qreg.R

In carter-allen/lbayes:

#' Bayesian quantile regression
#'
#' This function allows you to fit Bayesian QR.
#' @param formula Model formula object
#' @param data Data frame to use
#' @param quantile Quantile to model
#' @keywords Quantile regression
#' @export qreg
#' @examples


#### Helper function ####
# rho function as defined in Reed and Yu
rho <- function(u,p)
{
    return(ifelse(u >= 0,p*abs(u),(1-p)*abs(u)))
}

qreg <- function(formula,data,quantile = 0.5)
{
    f <- as.formula(formula)
    y <- model.frame(formula, data)[,1]
    X <- model.matrix(formula, data)
    n <- length(y)
    p <- quantile
    
    #### Priors ####
    c <- d <- 0.001 # gamma priors for tau
    k <- ncol(X) # mvnorm dimension
    b0 <- rep(0,k) # prior mean of beta
    T0 <- diag(0.01,k) # prior precision of beta
    
    #### Initialize ####
    beta <- rep(0,k)
    
    #### Sample Storage ####
    nsim <- 1000 # Number of MCMC Iterations
    thin <- 1	# Thinning interval
    burn <- 0	# Burnin
    lastit <- (nsim-burn)/thin # Last stored value
    Beta <- matrix(0,lastit,k) # storage for each beta
    Sigma2<-rep(0,lastit) # storage for each 1/tau
    
    #### Gibbs ####
    start.time<-proc.time()
    print(paste("Started MCMC of",nsim,"iterations with a burn in of",burn))
    pb <- txtProgressBar(min = 0, max = nsim, style = 3)
    for(i in 1:lastit)
    {
        # Update tau
        c.i <- c + n
        d.i <- d + sum(rho(y - X %*% beta,p))
        tau <- rgamma(1,c.i,d.i)
        
        # Update v
        v.vec <- rinvgauss(n,1/(p*(1-p)*abs(y - X %*% beta)),tau/(2*p*(1-p)))
        v <- diag(v.vec)
        
        # Compute u
        u = y - (1-2*p)/(p*(1-p)*v.vec)
        
        # Update beta
        V <- solve((tau*p*(1-p)/2)*t(X)%*%v%*%X + T0)
        M <- V %*% ((tau*p*(1-p)/2) * t(X)%*%v%*%u+T0%*%b0)
        beta <- c(rmvnorm(1,M,V))
        
        if (i> burn & i%%thin==0) 
        {
            j<-(i-burn)/thin
            Beta[j,]<-beta 
            Sigma2[j]<-1/tau
        }
        setTxtProgressBar(pb, i)
    }
    close(pb)
    run.time<-proc.time()-start.time
    print(paste("Finished MCMC after",run.time[1],"seconds"))
    return(Beta)
}