# qrjoint: Joint Estimation of Linear Quantile Planes In qrjoint: Joint Estimation in Linear Quantile Regression

## Description

Estimate intercept and slope functions within a joint linear regression model of the quantiles, with possible right censoring of the response.

## Usage

 ``` 1 2 3 4 5 6 7 8 9 10``` ``` qrjoint(x, y, nsamp = 1e3, thin = 10, cens = rep(0,length(y)), incr = 0.01, par = "prior", nknots = 6, hyper = list(sig = c(.1,.1), lam = c(6,4), kap = c(.1,.1,1)), shrink = FALSE, prox.range = c(.2,.95), acpt.target = 0.15, ref.size = 3, blocking = "std5", temp = 1, expo = 2, blocks.mu, blocks.S, fix.nu=FALSE) ## S3 method for class 'qrjoint' update(object, nadd, append = TRUE, ...) ```

## Arguments

 `x` numeric design matrix of explanatory variables (do not include a column of ones). Could be given as a data frame. Rows are observations and columns are variables. Could be given as a dimensionless vector when there is a single predictor. `y` numeric vector of response data. `nsamp` number of posterior samples to be saved; defaults to 1000. `thin` thinning rate for the Markov chain sampler – one posterior sample is saved per `thin` iterations. Defaults to 10. The Markov chain sampler runs for a total of `nsamp * thin` many iterations. `cens` censoring status of response. Must be a vector of 0s and 1s of length same as length(y), with 1 indicating right censoring, and, 0 indicating no censoring. Defaults to all zeros. `incr` tau grid increment. Defaults to 0.01. `par` character string indicating how the sampler is to be initialized. Only two options are currently supported: "prior" to initialize at a random draw from the prior; "RQ" to initialize at a model space approximation of the estimates from `rq`. `nknots` number of knots to be used for low rank approximation of the Gaussian process priors. Defaults to 6. `hyper` hyperparameters of the prior distribution. Must be a list with some of all of the following fields: `sig`: a two vector giving the parameters of the inverse-gamma distribution on sigma-square that is used when `shrink`=TRUE, `lam`: a two vector giving the parameters of the beta distribution on proximity = exp(-0.01 * lambda^2), and `kap`: a vector to be coerced into a `3 * nkap` matrix, with `nkap` being the number of components in the mixture of gamma prior on `kappa`, and each column of the matrix gives the shape, rate and mixing weight of a component. `shrink` for applying shrinkage to gamma[0] and gamma. Defaults to FALSE, in which case a right Haar prior is used on (gamma[0], gamma, sigma2). If TRUE then a horseshoe type prior is used. `prox.range` for specifying the range of length-scale parameter of the Gaussian process prior. `acpt.target` target acceptance rate of the adaptive Metropolis sampler; defaults to 0.15 `ref.size` adaptation rate of the adaptive Metropolis sampler. The proposal density is updated once every `ref.size` iterations. Could be a single number or a vector of length same as the number of blocks. `blocking` type of blocking to be applied. Either a character string specifying one to be chosen from the supplied menu (see Details), or a list giving user specified blocks. In the latter case, each element of the list is a logical vector of length equal to the total number of model parameters, which equals `(nknots+1)*(ncol(X)+1) + 2` indicating which model parameters belong to the block. `temp` temperature of the log-likelihood function. The log-likelihood function is raised to the power of `temp`. Defaults to 1. `expo` the exponent to be used in the covariance kernel of the Gaussian process priors. Defaults to 2, giving the standard squared-exponential covariance kernel. `blocks.mu` initial block specific means in the form of a list. If left unspecified then will be automatically generated as a list of vectors of zeros of appropriate lengths matching the corresponding block sizes. `blocks.S` initial block specific covariance matrices in the form of a list. If left unspecified then will be automatically generated as a list of identity matrices of appropriate dimensions matching the corresponding block sizes. When `blocking` is chosen as one of the menu items of the form "std*", known prior covariance information and estimated variance matrices from `rq` are used. `fix.nu` either the logical FALSE indicating that nu should be learned, or a positive real number giving the fixed value of nu, which is then excluded from MCMC updates `object` a fitted model of the class 'qrjoint'. `nadd` number of additional MCMC samples. `append` logical indicating whether new samples should be appended to old ones. If FALSE then old samples are discarded. `...` no additional arguments are allowed

## Details

The following menu choices are available for blocking the parameter vector. Below, `p = ncol(X)`.

`"single"`: a single block containing all parameters

`"single2"`: one block containing all parameters and an additional block containing only (gamma[0], gamma, sigma, nu)

`"single3"`: like `"single2"`, but the second block is split into two further blocks, one with (gamma[0], gamma), the other with (sigma, nu)

`"std0"`: `p+1` blocks, (j+1)-th contains (W*[j], gamma[j], sigma, nu), j = 0, ..., p.

`"std1"`: total `p+2` blocks. First `p+1` blocks same as `"std0"` and one additional block of (gamma[0], gamma, sigma, nu).

`"std2"`: total `p+3` blocks. First `p+1` blocks same as `"std0"` and two additional blocks of (gamma[0], gamma) and (sigma, nu)

`"std3"`: total `p+3` blocks. First `p+1` blocks are W*[j], j=0,...,p, last two are (gamma[0], gamma) and (sigma, nu)

`"std4"`: total `p+3` blocks. First `p+1` blocks are (W*[j], gamma[j]), j=0,...,p, last two are (gamma[0], gamma) and (sigma, nu)

`"std5"`: total `p+4` blocks. First `p+3` are same as `"std4"` and one additional block containing all parameters.

## Value

`qrjoint(x, y, ...)` returns a ‘qrjoint’ class object to be used by `update`, `coef` and `summary`.

`update(object,...)` runs additional Markov chain iterations and appends thinned draws to an existing ‘qrjoint’ object `object`. All relevant details are inherited from `object`.

Returned object is a list containing the following variables.

 `par` latest draw of the parameter vector `x` scaled and centered design matrix `y` response vector `cens` censoring status vector `shrink` shrinkage indicator `hyper` completed list of hyper-parameters `dim` model dimension vector of the form c(n, p, length of tau grid, position of tau0 on the grid, nknots, length of lambda grid, nkap, total number of MCMC iterations, thin, nsamp) `gridmats` details of covariance matrix factors etc, intended for internal use. `tau.g` the tau grid `muV` list of means for parameter blocks `SV` list of covariance matrices for parameter blocks `blocks` list of blocks `blocks.size` vector of block lengths `dmcmcpar` numeric vector containing details of adaptive MCMC runs, equals c(temp, decay rate of adaptation, vector of target acceptance rates for the blocks, vector of increment scales used in adaptation). Intended strictly for internal use. `imcmcpar` numeric vector containing details of adaptive MCMC runs, equals c(number of parameter blocks, ref.size, indicator on whether details are to be printed during MCMC progress, rate of details printing, a vector of counters needed for printing). Intended strictly for internal use. `parsamp` a long vector containing the parameter draws. Could be coerced into a matrix of dim `npar * nsamp`. Intended primarily for use by `summary` and `coef`. `acptsamp` a long vector containing rates of acceptance statistics for parameter blocks. Could be coerced into a matrix of dim `nblocks * nsamp`. Not very informative, because thinning times and adaptation times may not be exactly synced. `lpsamp` vector of log posterior values for the saved MCMC draws. `prox` vector of proximity (exp(-0.01*lambda^2)) grid values `reg.ix` positions of the regular tau grid on the expanded tail-appended grid `runtime` run time of the MCMC

## References

Yang, Y., and Tokdar, S. T. (2015). Joint Estimation of Quantile Planes over Arbitrary Predictor Spaces.

`summary.qrjoint` and `coef.qrjoint`.
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66``` ``` ## Plasma data analysis # recoding variables data(plasma) plasma\$Sex <- as.factor(plasma\$Sex) plasma\$SmokStat <- as.factor(plasma\$SmokStat) plasma\$VitUse <- 3 - plasma\$VitUse plasma\$VitUse <- as.factor(plasma\$VitUse) # creating predictors and response (beta carotene concentration in the plasma) X <- model.matrix(BetaPlasma ~ Age + Sex + SmokStat + Quetelet + VitUse + Calories + Fat + Fiber + Alcohol + Cholesterol + BetaDiet, data = plasma)[,-1] Y <- plasma\$BetaPlasma # model fitting with 40 posterior samples from 80 iterations (thin = 2) # this is of course for illustration, for practical model fitting you # ought to try at least something like nsamp = 500, thin = 20 fit.qrj <- qrjoint(X, Y, nsamp = 40, thin = 2) summary(fit.qrj, more = TRUE) ## Not run: # additional MCMC runs to get 10 more samples (20 additional iterations) fit.qrj <- update(fit.qrj, 10) summary(fit.qrj, more = TRUE) ## End(Not run) ## Not run: ### UIS data analysis (with right censoring) data(uis) y <- uis\$Y cens <- 1 - uis\$CENSOR X <- model.matrix(~ TREAT + NDT + IV3 + BECK + FRAC + RACE + AGE + SITE, data = uis)[,-1] uis.qrj <- qrjoint(X, y, cens = cens, nsamp = 50, thin = 2, fix.nu = 1e5) summary(uis.qrj, more = TRUE) betas <- coef(uis.qrj, plot = TRUE, col = "darkgreen") tau.grid <- uis.qrj\$tau.g[uis.qrj\$reg.ix] L <- length(tau.grid) beta.samp <- betas\$beta.samp # survival curve estimation for 9 randomly chosen subjects n <- nrow(X) ix.sel <- sort(sample(n, 9)) Xsel <- X[ix.sel,] Qsel.gp <- apply(beta.samp, 2, function(b) return(tcrossprod(matrix(b, nrow = L), cbind(1, Xsel)))) k <- length(ix.sel) colRGB <- col2rgb("darkgreen")/255 colTrans <- rgb(colRGB[1], colRGB[2], colRGB[3], alpha = 0.05) par(mfrow = c(3,3), mar = c(4,3,2,1) + .1) for(i in 1:k){ plot(exp(Qsel.gp[(i-1)*L + 1:L,1]), 1 - tau.grid, ty = "n", ann = FALSE, bty = "n", xlim = exp(c(2, 8)), ylim = c(0,1), lty = 2, log = "x") for(j in 1:ncol(beta.samp)) lines(exp(Qsel.gp[(i-1)*L + 1:L,j]), 1 - tau.grid, col = colTrans, lwd = 1) title(xlab = "Return time (days)", ylab = "Survival function", line = 2) title(main = bquote(Obs.Id == .(ix.sel[i]))) grid() } ## End(Not run) ```