rhierLinearDPParallel: MCMC Algorithm for Hierarchical Linear Model with Dirichlet...

In scalablebayesm: Distributed Markov Chain Monte Carlo for Bayesian Inference in Marketing

MCMC Algorithm for Hierarchical Linear Model with Dirichlet Process Prior Heterogeneity

Description

rhierLinearDPParallel is an MCMC algorithm for a hierarchical linear model with a Dirichlet Process prior for the distribution of heterogeneity. A base normal model is used so that the DP can be interpreted as allowing for a mixture of normals with as many components as panel units. The function implements a Gibbs sampler on the coefficients for each panel unit. This procedure can be interpreted as a Bayesian semi-parametric method since the DP prior accommodates heterogeneity of unknown form.

Usage

rhierLinearDPParallel(Data, Prior, Mcmc, verbose = FALSE)

Arguments

Data	A list of: 'regdata' - A nreg size list of regdata, and optional 'Z'- nreg \times nz matrix of unit characteristics.
Prior	A list with one required parameter: 'ncomp', and optional parameters: 'deltabar', 'Ad', 'mubar', 'Amu', 'nu', 'V', 'nu.e', and 'ssq'.
Mcmc	A list with one required parameter: 'R'-number of iterations, and optional parameters: 'keep' and 'nprint'.
verbose	If TRUE, enumerates model parameters and timing information.

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 \times nvar matrix

u_i \sim N(\mu_{ind}, \Sigma_{ind})

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)

\theta_i = (\mu_i, \Sigma_i) \sim DP(G_0(\lambda), alpha)

G_0(\lambda):
\mu_i | \Sigma_i \sim N(0, \Sigma_i (x) a^{-1})
\Sigma_i \sim IW(nu, nu*v*I)
delta = vec(\Delta) \sim N(deltabar, A_d^{-1})

\lambda(a, nu, v):
a \sim uniform[alim[1], alimb[2]]
nu \sim dim(data)-1 + exp(z)
z \sim uniform[dim(data)-1+nulim[1], nulim[2]]
v \sim uniform[vlim[1], vlim[2]]

Z should not include an intercept and should be centered for ease of interpretation. The mean of each of the nreg \betas is the mean of the normal mixture. Use summary() to compute this mean from the compdraw output.

The prior on \Sigma_i is parameterized such that mode(\Sigma) = nu/(nu+2) vI. The support of nu enforces a non-degenerate IW density; nulim[1] > 0

The default choices of alim, nulim, and vlim determine the location and approximate size of candidate "atoms" or possible normal components. The defaults are sensible given a reasonable scaling of the X variables. You want to insure that alim is set for a wide enough range of values (remember a is a precision parameter) and the v is big enough to propose Sigma matrices wide enough to cover the data range.

A careful analyst should look at the posterior distribution of a, nu, v to make sure that the support is set correctly in alim, nulim, vlim. In other words, if we see the posterior bunched up at one end of these support ranges, we should widen the range and rerun.

If you want to force the procedure to use many small atoms, then set nulim to consider only large values and set vlim to consider only small scaling constants. Set alphamax to a large number. This will create a very "lumpy" density estimate somewhat like the classical Kernel density estimates. Of course, this is not advised if you have a prior belief that densities are relatively smooth.

Argument Details

Data = list(regdata, Z) [Z optional]

regdata:	A nreg/s_shard size list of regdata
regdata[[i]]$X:	n_i \times nvar design matrix for equation i
regdata[[i]]$y:	n_i \times 1 vector of observations for equation i
Z:	A list of s partitions where each partition include (nreg/s_shard) \times nz matrix of unit characteristics

Prior = list(deltabar, Ad, Prioralphalist, lambda_hyper, nu, V, nu_e, mubar, Amu, ssq, ncomp) [all but ncomp are optional]

deltabar:	(nz \times nvar) \times 1 vector of prior means (def: 0)
Ad:	prior precision matrix for vec(D) (def: 0.01*I)
mubar:	nvar \times 1 prior mean vector for normal component mean (def: 0)
Amu:	prior precision for normal component mean (def: 0.01)
nu.e:	d.f. parameter for regression error variance prior (def: 3)
V:	PDS location parameter for IW prior on normal component Sigma (def: nu*I)
ssq:	scale parameter for regression error variance prior (def: var(y_i))
ncomp:	number of components used in normal mixture

Mcmc = list(R, keep, nprint) [only R required]

R:	number of MCMC draws
keep:	MCMC thinning parameter -- 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)

Value

A list containing:

compdraw	A list (length: R/keep) where each list contains 'mu' (vector, length: 'ncomp') and 'rooti' (matrix, shape: ncomp \times ncomp)
probdraw	A (R/keep) \times (ncomp) matrix that reports the probability that each draw came from a particular component
Deltadraw	A (R/keep) \times (nz \times nvar) matrix of draws of Delta, first row is initial value. Deltadraw is NULL if Z is NULL

Author(s)

Federico Bumbaca, Leeds School of Business, University of Colorado Boulder, [e-mail]

References

Bumbaca, Federico (Rico), Sanjog Misra, and Peter E. Rossi (2020), "Scalable Target Marketing: Distributed Markov Chain Monte Carlo for Bayesian Hierarchical Models", Journal of Marketing Research, 57(6), 999-1018.

Chapter 5, Bayesian Statistics and Marketing by Rossi, Allenby, and McCulloch.

See Also

partition_data, drawPosteriorParallel, combine_draws, rheteroLinearIndepMetrop

Examples

######### Single Component with rhierLinearDPParallel########
R = 1000

nreg = 2000
nobs = 5 #number of observations
nvar = 3 #columns
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)) 

#Default
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) 
tpvec=c(1)                               

regdata=NULL						  
betas=matrix(double(nreg*nvar),ncol=nvar) 
tind=double(nreg) 

for (reg in 1:nreg) {
 tempout=bayesm::rmixture(1,tpvec,tcomps) 
 if (is.null(Z)){
   betas[reg,]= as.vector(tempout$x)  
 }else{
   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)
}


Data1=list(list(regdata=regdata,Z=Z))
s = 1
Data2=scalablebayesm::partition_data(Data1,s=s)

Prior1=list(ncomp=1)
Mcmc1=list(R=R,keep=1)

out2 = parallel::mclapply(Data2, FUN = rhierLinearDPParallel, Prior = Prior1,
Mcmc = Mcmc1, mc.cores = s, mc.set.seed = FALSE)

######### Multiple Components with rhierLinearDPParallel########
R = 1000

set.seed(66)
nreg=2000
nobs=5 #number of observations
nvar=3 #columns
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)) 

#Default
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						  
betas=matrix(double(nreg*nvar),ncol=nvar) 
tind=double(nreg) 

for (reg in 1:nreg) {
 tempout=bayesm::rmixture(1,tpvec,tcomps) 
 if (is.null(Z)){
   betas[reg,]= as.vector(tempout$x)  
 }else{
   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)
}


Data1=list(list(regdata=regdata,Z=Z))
s = 1
Data2=scalablebayesm::partition_data(Data1, s=s)

Prior1=list(ncomp=3)
Mcmc1=list(R=R,keep=1)

out2 = parallel::mclapply(Data2, FUN = rhierLinearDPParallel, Prior = Prior1,
Mcmc = Mcmc1, mc.cores = s, mc.set.seed = FALSE)

scalablebayesm documentation built on April 3, 2025, 7:55 p.m.