rhierLinearMixture {bayesm} | R Documentation |
rhierLinearMixture
implements a Gibbs Sampler for hierarchical linear models with a mixture of normals prior.
rhierLinearMixture(Data, Prior, Mcmc)
Data |
list(regdata,Z) (Z optional). |
Prior |
list(deltabar,Ad,mubar,Amu,nu,V,nu.e,ssq,ncomp) (all but ncomp are optional). |
Mcmc |
list(R,keep) (R required). |
Model: length(regdata) regression equations.
y_i = X_ibeta_i + e_i. e_i ~ N(0,tau_i). nvar X vars in each equation.
Priors:
tau_i ~ nu.e*ssq_i/chi^2_{nu.e}. tau_i is the variance of e_i.
beta_i= ZDelta[i,] + u_i.
Note: here ZDelta refers to Z%*%D, ZDelta[i,] is ith row of this product.
Delta is an nz x nvar array.
u_i ~ N(mu_{ind},Sigma_{ind}). ind ~ multinomial(pvec).
pvec ~ dirichlet (a)
delta= vec(Delta) ~ N(deltabar,A_d^{-1})
mu_j ~ N(mubar,Sigma_j (x) Amu^{-1})
Sigma_j ~ IW(nu,V)
List arguments contain:
regdata
regdata[[i]]$X
regdata[[i]]$y
deltabar
Ad
mubar
Amu
nu
V
nu.e
ssq
ncomp
R
keep
a list containing
taudraw |
R/keep x nreg array of error variance draws |
betadraw |
nreg x nvar x R/keep array of individual regression coef draws |
Deltadraw |
R/keep x nz x nvar array of Deltadraws |
nmix |
list of three elements, (probdraw, NULL, compdraw) |
More on probdraw
component of nmix return value list:
this is an R/keep by ncomp array of draws of mixture component probs (pvec)
More on compdraw
component of nmix return value list:
Note: Z should not include an intercept and should be centered for ease of interpretation.
Be careful in assessing prior parameter, Amu. .01 can be too small for some applications. See
Rossi et al, chapter 5 for full discussion.
Peter Rossi, Graduate School of Business, University of Chicago, Peter.Rossi@ChicagoGsb.edu.
For further discussion, see Bayesian Statistics and Marketing
by Rossi, Allenby and McCulloch, Chapter 3.
http://faculty.chicagogsb.edu/peter.rossi/research/bsm.html
## 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=.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(.4,.2,.4) regdata=NULL # simulated data with Z 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) if(0){ ## plotting examples plot(out1$betadraw) plot(out1$nmix) plot(out1$Deltadraw) }