Gibbs Sampler for Hierarchical Linear Model with Mixture-of-Normals Heterogeneity
rhierLinearMixture implements a Gibbs Sampler for hierarchical linear models with a mixture-of-normals prior.
rhierLinearMixture(Data, Prior, Mcmc)
Data |
list(regdata, Z) |
Prior |
list(deltabar, Ad, mubar, Amu, nu, V, nu.e, ssq, ncomp) |
Mcmc |
list(R, keep, nprint) |
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.
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 Detailsnmix 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.
|
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: |
Peter Rossi, Anderson School, UCLA, perossichi@gmail.com.
For further discussion, see Chapter 5, Bayesian Statistics and Marketing by Rossi, Allenby, and McCulloch.
http://www.perossi.org/home/bsm-1
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)
}Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.