bayesm (version 1.1-2)

rhierLinearModel: Gibbs Sampler for Hierarchical Linear Model

Description

rhierLinearModel implements a Gibbs Sampler for hierarchical linear models.

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) (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

concept

• bayes
• MCMC
• Gibbs Sampling
• hierarchical models
• linear model

Details

Model: length(regdata) regression equations. $y_i = X_ibeta_i + e_i$. $e_i$ $\sim$ $N(0,tau_i)$. nvar X vars in each equation. Priors: $tau_i$ $\sim$ nu.e*$ssq_i/\chi^2_{nu.e}$. $tau_i$ is the variance of $e_i$. $beta_i$ $\sim$ N(ZDelta[i,],$V_{beta}$). Note: ZDelta is the matrix Z * Delta; [i,] refers to ith row of this product. $vec(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, set Z=iota (nreg x 1). 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)}

References

For further discussion, see Bayesian Statistics and Marketing by Allenby, McCulloch, and Rossi, Chapter 3. http://gsbwww.uchicago.edu/fac/peter.rossi/research/bsm.html

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) }

nu.e=3
ssq=NULL
for(reg in 1:nreg) {ssq[reg]=var(regdata[[reg]]$y)}
nu=nvar+3
V=nu*diag(c(rep(1,nvar)))
A=diag(c(rep(.01,nz)),ncol=nz)
Deltabar=matrix(c(rep(0,nz*nvar)),nrow=nz)


Data=list(regdata=regdata,Z=Z)
Prior=list(Deltabar=Deltabar,A=A,nu.e=nu.e,ssq=ssq,nu=nu,V=V)
Mcmc=list(R=R,keep=1)
out=rhierLinearModel(Data=Data,Mcmc=Mcmc)

cat("Deltadraws ",fill=TRUE)
mat=apply(out$Deltadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
mat=rbind(as.vector(Delta),mat); rownames(mat)[1]="delta"; print(mat)
cat("Vbetadraws ",fill=TRUE)
mat=apply(out$Vbetadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
mat=rbind(as.vector(Vbeta),mat); rownames(mat)[1]="Vbeta"; print(mat)