bayesm (version 3.0-2)

rhierLinearModel: Gibbs Sampler for Hierarchical Linear Model

Description

rhierLinearModel implements a Gibbs Sampler for hierarchical linear models with a normal prior.

Usage

rhierLinearModel(Data, Prior, Mcmc)

Arguments

Data

list(regdata,Z) (Z optional).

Prior

list(Deltabar,A,nu.e,ssq,nu,V) (optional).

Mcmc

list(R,keep,nprint) (R required).

Value

a list containing

betadraw

nreg x nvar x R/keep array of individual regression coef draws

taudraw

R/keep x nreg array of error variance draws

Deltadraw

R/keep x nz x nvar array of Deltadraws

Vbetadraw

R/keep x nvar*nvar array of Vbeta draws

Details

Model: length(regdata) regression equations. $y_i=X_i{\beta }_i+e_i$. $e_i$ $\sim$ $N(0,{\tau }_i)$. nvar X vars in each equation.

Priors: ${\tau }_i$ $\sim$ nu.e*$ss{q}_i/{\chi }_{nu.e}^2$. ${\tau }_i$ is the variance of $e_i$. ${\beta }_i$ $\sim$ N(Z$\mathrm{\Delta }$[i,],$V_{\beta }$). Note: Z$\mathrm{\Delta }$ is the matrix Z * $\mathrm{\Delta }$; [i,] refers to ith row of this product.

$vec(\mathrm{\Delta })$ given $V_{\beta }$ $\sim$ $N(vec(Deltabar),V_{\beta }(x)A^{-1})$. $V_{\beta }$ $\sim$ $IW(nu,V)$. $Delta,Deltabar$ are nz x nvar. $A$ is nz x nz. $V_{\beta }$ is nvar x nvar.

Note: if you don't have any Z vars, omit Z in the Data argument and a vector of ones will be inserted for you. In this case (of no Z vars), the matrix $\mathrm{\Delta }$ will be 1 x nvar and should be interpreted as the mean of all unit $\beta$ s.

List arguments contain:

• regdata list of lists with X,y matrices for each of length(regdata) regressions

• regdata[[i]]$X X matrix for equation i

• regdata[[i]]$y y vector for equation i

• Deltabar nz x nvar matrix of prior means (def: 0)

• A nz x nz matrix for prior precision (def: .01I)

• nu.e d.f. parm for regression error variance prior (def: 3)

• ssq scale parm for regression error var prior (def: var($y_i$))

• nu d.f. parm for Vbeta prior (def: nvar+3)

• V Scale location matrix for Vbeta prior (def: nu*I)

• R number of MCMC draws

• keep MCMC thinning parm: keep every keepth draw (def: 1)

• nprint print the estimated time remaining for every nprint'th draw (def: 100)

References

For further discussion, see Bayesian Statistics and Marketing by Rossi, Allenby and McCulloch, Chapter 3. http://www.perossi.org/home/bsm-1

See Also

rhierLinearMixture

Examples

##
if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10}

nreg=100; nobs=100; nvar=3
Vbeta=matrix(c(1,.5,0,.5,2,.7,0,.7,1),ncol=3)
Z=cbind(c(rep(1,nreg)),3*runif(nreg)); Z[,2]=Z[,2]-mean(Z[,2])
nz=ncol(Z)
Delta=matrix(c(1,-1,2,0,1,0),ncol=2)
Delta=t(Delta) # first row of Delta is means of betas
Beta=matrix(rnorm(nreg*nvar),nrow=nreg)%*%chol(Vbeta)+Z%*%Delta
tau=.1
iota=c(rep(1,nobs))
regdata=NULL
for (reg in 1:nreg) { X=cbind(iota,matrix(runif(nobs*(nvar-1)),ncol=(nvar-1)))
	y=X%*%Beta[reg,]+sqrt(tau)*rnorm(nobs); regdata[[reg]]=list(y=y,X=X) }

Data1=list(regdata=regdata,Z=Z)
Mcmc1=list(R=R,keep=1)
out=rhierLinearModel(Data=Data1,Mcmc=Mcmc1)

cat("Summary of Delta draws",fill=TRUE)
summary(out$Deltadraw,tvalues=as.vector(Delta))
cat("Summary of Vbeta draws",fill=TRUE)
summary(out$Vbetadraw,tvalues=as.vector(Vbeta[upper.tri(Vbeta,diag=TRUE)]))

if(0){
## plotting examples
plot(out$betadraw)
plot(out$Deltadraw)
}