sur_sample_powerprior {surbayes}

Sample from SUR posterior via power prior

Description

This function uses Gibbs sampling to sample from the posterior density of a SUR model using the power prior.

Usage

sur_sample_powerprior(
  formula.list,
  data,
  histdata,
  M,
  Sigma0 = NULL,
  a0 = 1,
  burnin = 0,
  thin = 1
)

Arguments

formula.list

A list of formulas, each element giving the formula for the corresponding endpoint.

data

A data.frame containing all the variables contained in formula.list]

histdata

A data.frame of historical data to apply power prior on

M

Number of samples to be drawn

Sigma0

A J \times J matrix giving the initial covariance matrix. Default is the MLE.

a0

A scalar between 0 and 1 giving the power prior parameter

burnin

A non-negative integer giving the burn-in parameter

thin

A positive integer giving the thin parameter

Value

A list. First element is posterior draws. Second element is list of JxJ covariance matrices.

Examples

## Taken from bayesm package
if(nchar(Sys.getenv("LONG_TEST")) != 0) {M=1000} else {M=10}
set.seed(66)
## simulate data from SUR
beta1 = c(1,2)
beta2 = c(1,-1,-2)
nobs = 100
nreg = 2
iota = c(rep(1, nobs))
X1 = cbind(iota, runif(nobs))
X2 = cbind(iota, runif(nobs), runif(nobs))
Sigma = matrix(c(0.5, 0.2, 0.2, 0.5), ncol = 2)
U = chol(Sigma)
E = matrix( rnorm( 2 * nobs ), ncol = 2) %*% U
y1 = X1 %*% beta1 + E[,1]
y2 = X2 %*% beta2 + E[,2]
X1 = X1[, -1]
X2 = X2[, -1]
data = data.frame(y1, y2, X1, X2)
names(data) = c( paste0( 'y', 1:2 ), paste0('x', 1:(ncol(data) - 2) ))
## run DMC sampler
formula.list = list(y1 ~ x1, y2 ~ x2 + x3)

## fit using historical data as current data set--never done in practice
out = sur_sample_powerprior( formula.list, data, histdata = data, M = M )