rhierLinearModel {bayesm} | R Documentation |
rhierLinearModel
implements a Gibbs Sampler for hierarchical linear models with a normal prior.
rhierLinearModel(Data, Prior, Mcmc)
Data |
list(regdata,Z) (Z optional). |
Prior |
list(Deltabar,A,nu.e,ssq,nu,V) (optional). |
Mcmc |
list(R,keep) (R required). |
Model: length(regdata) regression equations.
y_i = X_ibeta_i + e_i. e_i ~ N(0,tau_i). nvar X vars in each equation.
Priors:
tau_i ~ nu.e*ssq_i/chi^2_{nu.e}. tau_i is the variance of e_i.
beta_i ~ 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} ~ N(vec(Deltabar),V_{beta} (x) A^{-1}).
V_{beta} ~ 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
regdata[[i]]$X
regdata[[i]]$y
Deltabar
A
nu.e
ssq
nu
V
R
keep
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 |
Peter Rossi, Graduate School of Business, University of Chicago, Peter.Rossi@ChicagoGsb.edu.
For further discussion, see Bayesian Statistics and Marketing
by Rossi, Allenby and McCulloch, Chapter 3.
http://faculty.chicagogsb.edu/peter.rossi/research/bsm.html
## 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) }