bayesm (version 3.0-2)

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

Details

Model: length(regdata) regression equations. \(y_i = X_i\beta_i + e_i\). \(e_i\) \(\sim\) \(N(0,\tau_i)\). nvar is the number of 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\). \(B = Z\Delta + U\) or \(\beta_i = \Delta' Z[i,]' + u_i\). \(\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

• deltabarnz*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\)))

• a Dirichlet prior parameter (def: 5)

• ncomp number of components used in normal mixture

• 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

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