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\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 .
|
Value
A list containing:
taudraw |
|
betadraw |
|
Deltadraw |
|
nmix |
a list containing: |
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.
See Also
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)
}