rhierLinearMixture {bayesm}

Gibbs Sampler for Hierarchical Linear Model with Mixture-of-Normals Heterogeneity

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)

Prior

list(deltabar, Ad, mubar, Amu, nu, V, nu.e, ssq, ncomp)

Mcmc

list(R, keep, nprint)

Details

Model and Priors

nreg regression equations with nvar as the number of X vars in each equation
y_i = X_i\beta_i + e_i with e_i \sim N(0, \tau_i)

\tau_i \sim nu.e*ssq_i/\chi^2_{nu.e} where \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 matrix

Z should not include an intercept and should be centered for ease of interpretation. The mean of each of the nreg \betas is the mean of the normal mixture. Use summary() to compute this mean from the compdraw output.

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)

Be careful in assessing the prior parameter Amu: 0.01 can be too small for some applications. See chapter 5 of Rossi et al for full discussion.

Argument Details

Data = list(regdata, Z) [Z optional]

regdata:	list of lists with X and y matrices for each of nreg=length(regdata) regressions
regdata[[i]]$X:	n_i x nvar design matrix for equation i
regdata[[i]]$y:	n_i x 1 vector of observations for equation i
Z:	nreg x nz matrix of unit characteristics (def: vector of ones)

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

deltabar:	nz x nvar vector of prior means (def: 0)
Ad:	prior precision matrix for vec(Delta) (def: 0.01*I)
mubar:	nvar x 1 prior mean vector for normal component mean (def: 0)
Amu:	prior precision for normal component mean (def: 0.01)
nu:	d.f. parameter for IW prior on normal component Sigma (def: nvar+3)
V:	PDS location parameter for IW prior on normal component Sigma (def: nu*I)
nu.e:	d.f. parameter for regression error variance prior (def: 3)
ssq:	scale parameter for regression error variance prior (def: var(y_i))
a:	Dirichlet prior parameter (def: 5)
ncomp:	number of components used in normal mixture

Mcmc = list(R, keep, nprint) [only R required]

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, set to 0 for no print)

nmix Details

nmix is a list with 3 components. Several functions in the bayesm package that involve a Dirichlet Process or mixture-of-normals return nmix. Across these functions, a common structure is used for nmix in order to utilize generic summary and plotting functions.

probdraw:	ncomp x R/keep matrix that reports the probability that each draw came from a particular component
zdraw:	R/keep x nobs matrix that indicates which component each draw is assigned to (here, null)
compdraw:	A list of R/keep lists of ncomp lists. Each of the inner-most lists has 2 elemens: a vector of draws for mu and a matrix of draws for the Cholesky root of Sigma.

Value

A list containing:

taudraw

R/keep x nreg matrix of error variance draws

betadraw

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

Deltadraw

R/keep x nz*nvar matrix of Deltadraws

nmix

a list containing: probdraw, zdraw, compdraw (see “nmix Details” section)

Author(s)

Peter Rossi, Anderson School, UCLA, [e-mail].

References

For further discussion, see Chapter 5, Bayesian Statistics and Marketing by Rossi, Allenby, and McCulloch.

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 = 0.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(0.4, 0.2, 0.4)

## simulated data with Z
regdata = NULL
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)

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