rhierLinearMixture: 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)

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)

Arguments

Data: list(regdata, Z)
Prior: list(deltabar, Ad, mubar, Amu, nu, V, nu.e, ssq, ncomp)
Mcmc: list(R, keep, nprint)

Author

Peter Rossi, Anderson School, UCLA, perossichi@gmail.com.

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 $\beta$s 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 `keep`th 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`.

References

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

Examples

Run this code

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

Run the code above in your browser using DataLab

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

`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

`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`.