bayesm (version 2.2-1)

rhierLinearMixture: Gibbs Sampler for Hierarchical Linear Model

Description

rhierLinearMixture implements a Gibbs Sampler for hierarchical linear models with a mixture of normals prior.

Usage

rhierLinearMixture(Data, Prior, Mcmc)

Arguments

Data

list(regdata,Z) (Z optional).

Prior

list(deltabar,Ad,mubar,Amu,nu,V,nu.e,ssq,ncomp) (all but ncomp are optional).

Mcmc

list(R,keep) (R required).

Value

• a list containing
• taudraw
  R/keep x nreg array of error variance draws
• betadraw
  nreg x nvar x R/keep array of individual regression coef draws
• Deltadraw
  R/keep x nz x nvar array of Deltadraws
• nmix
  list of three elements, (probdraw, NULL, compdraw)

concept

• bayes
• MCMC
• Gibbs Sampling
• mixture of normals
• normal mixture
• heterogeneity
• regresssion
• 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$= ZDelta[i,] + $u_i$. Note: here ZDelta refers to Z%*%D, ZDelta[i,] is ith row of this product. Delta is an nz x nvar array. $u_i$ $\sim$ $N(mu_{ind},Sigma_{ind})$. $ind$ $\sim$ multinomial(pvec). $pvec$ $\sim$ dirichlet (a) $delta= vec(Delta)$ $\sim$ $N(deltabar,A_d^{-1})$ $mu_j$ $\sim$ $N(mubar,Sigma_j (x) Amu^{-1})$ $Sigma_j$ $\sim$ IW(nu,V) 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*nvar vector of prior means (def: 0)} Ad{ prior prec matrix for vec(Delta) (def: .01I)} mubar{ nvar x 1 prior mean vector for normal comp mean (def: 0)} Amu{ prior precision for normal comp mean (def: .01I)} nu{ d.f. parm for IW prior on norm comp Sigma (def: nvar+3)} V{ pds location parm for IW prior on norm comp Sigma (def: nuI)} nu.e{ d.f. parm for regression error variance prior (def: 3)} ssq{ scale parm for regression error var prior (def: var($y_i$))} ncomp{ number of components used in normal mixture } 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 Rossi, Allenby and McCulloch, Chapter 3. http://faculty.chicagogsb.edu/peter.rossi/research/bsm.html

See Also

rhierLinearModel

Examples

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

set.seed(66)
nreg=300; nobs=500; nvar=3; nz=2

Z=matrix(runif(nreg*nz),ncol=nz) 
Z=t(t(Z)-apply(Z,2,mean))
Delta=matrix(c(1,-1,2,0,1,0),ncol=nz)
tau0=.1
iota=c(rep(1,nobs))

## create arguments for rmixture

tcomps=NULL
a=matrix(c(1,0,0,0.5773503,1.1547005,0,-0.4082483,0.4082483,1.2247449),ncol=3)
tcomps[[1]]=list(mu=c(0,-1,-2),rooti=a) 
tcomps[[2]]=list(mu=c(0,-1,-2)*2,rooti=a)
tcomps[[3]]=list(mu=c(0,-1,-2)*4,rooti=a)
tpvec=c(.4,.2,.4)                               

regdata=NULL						  # simulated data with Z
betas=matrix(double(nreg*nvar),ncol=nvar)
tind=double(nreg)

for (reg in 1:nreg) {
tempout=rmixture(1,tpvec,tcomps)
betas[reg,]=Delta%*%Z[reg,]+as.vector(tempout$x)
tind[reg]=tempout$z
X=cbind(iota,matrix(runif(nobs*(nvar-1)),ncol=(nvar-1)))
tau=tau0*runif(1,min=0.5,max=1)
y=X%*%betas[reg,]+sqrt(tau)*rnorm(nobs)
regdata[[reg]]=list(y=y,X=X,beta=betas[reg,],tau=tau)
}

## run rhierLinearMixture

Data1=list(regdata=regdata,Z=Z)
Prior1=list(ncomp=3)
Mcmc1=list(R=R,keep=1)

out1=rhierLinearMixture(Data=Data1,Prior=Prior1,Mcmc=Mcmc1)

cat("Summary of Delta draws",fill=TRUE)
summary(out1$Deltadraw,tvalues=as.vector(Delta))
cat("Summary of Normal Mixture Distribution",fill=TRUE)
summary(out1$nmix)

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