blr: Bayesian Linear Regression

Description

This function fits a Bayesian linear regression model using scale invariant prior.

Usage

blr(formula, data = NULL, mcmc = list(), prior = list(), marginal.likelihood = TRUE)

Value

An object of class blm representing the Bayesian spectral analysis model fit. Generic functions such as print and fitted have methods to show the results of the fit.

The MCMC samples of the parameters in the model are stored in the list mcmc.draws and the posterior samples of the fitted values are stored in the list fit.draws. The output list also includes the following objects:

post.est: posterior estimates for all parameters in the model.
lmarg: log marginal likelihood.
rsquarey: correlation between $y$ and $\hat{y}$.
call: the matched call.
mcmctime: running time of Markov chain from system.time().

Arguments

formula: an object of class “formula”
data: an optional data frame.
mcmc: a list giving the MCMC parameters. The list includes the following integers (with default values in parentheses): nblow (1000) giving the number of MCMC in transition period, nskip (1) giving the thinning interval, smcmc (1000) giving the number of MCMC for analysis.
prior: a list giving the prior information. The list includes the following parameters (default values specify the non-informative prior): beta_m0 and beta_v0 giving the hyperparameters of the multivariate normal distribution for parametric part including intercept, sigma2_m0 and sigma2_v0 giving the prior mean and variance of the inverse gamma prior for the scale parameter of response.
marginal.likelihood: a logical variable indicating whether the log marginal likelihood is calculated.

Details

This generic function fits a Bayesian linear regression model using scale invariant prior.

Let $y_i$ and $w_i$ be the response and the vector of parametric predictors, respectively. The model for regression function is as follows.

$$y_i = w_i^T\beta + \epsilon_i, ~ i=1,\ldots,n, $$ where the error terms $\{\epsilon_i\}$ are a random sample from a normal distribution, $N(0,\sigma^2)$.

The conjugate priors are assumed for $\beta$ and $\sigma$: $$\beta | \sigma \sim N(m_{0,\beta}, \sigma^2V_{0,\beta}), \quad \sigma^2 \sim IG\Big(\frac{r_{0,\sigma}}{2}, \frac{s_{0,\sigma}}{2}\Big)$$

Examples

Run this code

#####################
# Simulated example #
#####################

# Simulate data
set.seed(1)

n <- 100
w <- runif(n)
y <- 3 + 2*w + rnorm(n, sd = 0.8)

# Fit the model with default priors and mcmc parameters
fout <- blr(y ~ w)

# Summary
print(fout); summary(fout)

# Fitted values
fit <- fitted(fout)

# Plots
plot(fout)