BLR: Bayesian Linear Regression

Description

The BLR (‘Bayesian Linear Regression’) function
was designed to fit parametric regression models using different
types of shrinkage methods. An earlier version of this program was presented in de los Campos et al. (2009).

Usage

1
2

Arguments

y

(numeric, n) the data-vector (NAs allowed).

XF

(numeric, n x pF) incidence matrix for bF, may be NULL.

XR

(numeric, n x pR) incidence matrix for bR, may be NULL.

XL

(numeric, n x pL) incidence matrix for bL, may be NULL.

GF

(list) providing an $ID (integer, n) linking observations to groups
(e.g., lines or sires) and a (co)variance structure ($A, numeric, pU x pU) between effects of the grouping factor
(e.g., line or sire effects). Note: ID must be an integer taking values from 1 to pU; ID[i]=q indicates that
the ith observation in y belongs to cluster q whose (co)variance function is in the qth row (column) of A.
GF may be NULL.

weights

(numeric, n) a vector of weights, may be NULL.

nIter,burnIn, thin

(integer) the number of iterations, burn-in and thinning.

saveAt

(string) this may include a path and a pre-fix that will be added to the name of the files that are saved as the program runs.

prior

(list) containing the following elements,

prior$varE, prior$varBR, prior$varU: (list) each providing degree of freedom ($df) and scale ($S). These are the parameters of the scaled inverse-χ^2
distributions assigned to variance components, see Eq. (2) below. In the parameterization used by BLR() the prior expectation of variance parameters is S/(df-2).

prior$lambda: (list) providing $value (initial value for lambda); $type (‘random’ or ‘fixed’) this argument specifies
whether lambda should be kept fixed at the value provided by $value or updated with samples from the posterior
distribution; and, either $shape and $rate (this when a Gamma prior is desired on lambda^2) or $shape1, $shape2 and
$max, in this case p(lambda|max,alpha1,alpha2) = Beta(lamda/max |alpha1,alpha2). For detailed description of these priors see de los Campos et al. (2009).

thin2

This value controls wether the running means are saved to disk or not. If thin2 is greater than nIter the running
means are not saved (default, thin2=1e10).

minAbsBeta

The minimum absolute value of the components of bL to avoid numeric problems when sampling from \boldsymbol τ^2, default 1e-9

Details

The program runs a Gibbs sampler for the Bayesian regression model described below.

Likelihood. The equation for the data is:

y=1*mu + XF*bF + XR*bR + XL*bL + Z*u + e ...(1)

where y, the response is a n x 1 vector (NAs allowed); mu is
an intercept; XF, XR, XL and Z are incidence matrices
used to accommodate different
types of effects (see below), and; e is a vector of model residuals assumed to be
distributed as e ~ MVN(0,Diag(varE/wi^2)),
here varE is an (unknown)
variance parameter and w_i are (known) weights that allow for heterogeneous-residual variances.

Any of the elements in the right-hand side of the linear predictor, except mu and e , can be omitted;
by default the program runs an intercept model.

Prior. The residual variance is assigned a scaled inverse-χ^2 prior with degree of freedom and scale parameter
provided by the user, that is, varE ~ Inv_Scaled_chisq(varE | dfe,Se). The regression coefficients {mu,bF,bR,bL,u} are assigned priors
that yield different type of shrinkage. The intercept and the vector of regression coefficients bF are assigned flat priors
(i.e., estimates are not shrunk). The vector of regression coefficients bR is assigned a
Gaussian prior with variance common to all effects, that is,
bRj ~ NIID(0,varBr). This prior is
the Bayesian counterpart of Ridge Regression. The variance parameter varBr,
is treated as unknown and it is assigned a scaled inverse-χ^2 prior, that is,
varBr ~ Inv_Scaled_chisq(varBr | dfBr,SBr) with degrees of freedom
dfBr, and scale SBr provided by the user.

The vector of regression coefficients bL is treated as in
the Bayesian LASSO of Park and Casella (2008). Specifically,

where, Exp(.|.) is an exponential prior and p(lambda) can either be: (a)
a mass-point at some value (i.e., fixed lambda); (b) p(lambda^2)~Gamma(r,delta) this
is the prior suggested by Park and Casella (2008); or, (c) p(lambda|max,alpha1,alpha2) = Beta(lamda/max |alpha1,alpha2), see de los Campos et al. (2009) for details. It can be shown that the marginal prior of regression coefficients bL_k, Integrate(N(bL_k|0,varE * tau_k^2) * Exp(tau_k^2 | lambda^2) d tau_k^2), is Double-Exponential. This prior has thicker tails and higher peak of mass at zero than the Gaussian prior used for bR, inducing a different type of shrinkage.

The vector u is used to model the so called ‘infinitesimal effects’, and is assigned a prior u~MVN(0,varU),
where, A is a positive-definite matrix (usually a relationship matrix computed from a pedigree) and varU is an unknow variance, whose prior is
varU ~ Inv_Scaled_chisq(varU | dfu,Su).

$nIter

the number of iterations made in the Gibbs sampler.

$burnIn

the nuber of iteratios used as burn-in.

$thin

the thin used.

$y

original data-vector.

The posterior means returned by BLR are calculated after burnIn is
passed and at a thin as specified by the user.

Save. The routine will save samples of mu, variance components and lambda and running means
(rm*.dat). Running means are computed using the thinning specified by
the user (see argument thin above); however these running means are
saved at a thinning specified by argument thin2 (by default, thin2=1e10 so that running means are computed
as the sampler runs but not saved to the disc).