rhierLinearMixture {bayesm} R Documentation

## 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β_i + e_i with e_i ~ N(0, τ_i)

τ_i ~ nu.e*ssq_i/χ^2_{nu.e} where τ_i is the variance of e_i
B = ZΔ + U or β_i = Δ' Z[i,]' + u_i
Δ 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` βs is the mean of the normal mixture. Use `summary()` to compute this mean from the `compdraw` output.

u_i ~ N(μ_{ind}, Σ_{ind})
ind ~ multinomial(pvec)

pvec ~ dirichlet(a)
delta = vec(Δ) ~ N(deltabar, A_d^{-1})
μ_j ~ N(mubar, Σ_j(x) Amu^{-1})
Σ_j ~ 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`.

### 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, perossichi@gmail.com.

### References

For further discussion, see Chapter 5, Bayesian Statistics and Marketing by Rossi, Allenby, and McCulloch.
http://www.perossi.org/home/bsm-1

`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[] = list(mu=c(0,-1,-2),   rooti=a)
tcomps[] = list(mu=c(0,-1,-2)*2, rooti=a)
tcomps[] = 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)