rhierMnlRwMixture {bayesm} | R Documentation |
MCMC Algorithm for Hierarchical Multinomial Logit with Mixture-of-Normals Heterogeneity
Description
rhierMnlRwMixture
is a MCMC algorithm for a hierarchical multinomial logit with a mixture of normals heterogeneity distribution. This is a hybrid Gibbs Sampler with a RW Metropolis step for the MNL coefficients for each panel unit.
Usage
rhierMnlRwMixture(Data, Prior, Mcmc)
Arguments
Data |
list(lgtdata, Z, p) |
Prior |
list(a, deltabar, Ad, mubar, Amu, nu, V, a, ncomp, SignRes) |
Mcmc |
list(R, keep, nprint, s, w) |
Details
Model and Priors
y_i
\sim
MNL(X_i,\beta_i)
with i = 1, \ldots,
length(lgtdata)
and where \beta_i
is nvar x 1
\beta_i
= Z\Delta
[i,] + u_i
Note: Z\Delta
is the matrix Z * \Delta
and [i,] refers to i
th row of this product
Delta is an nz x nvar
array
u_i
\sim
N(\mu_{ind},\Sigma_{ind})
with 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)
Note: Z
should NOT include an intercept and is centered for ease of interpretation.
The mean of each of the nlgt
\beta
s is the mean of the normal mixture.
Use summary()
to compute this mean from the compdraw
output.
Be careful in assessing prior parameter Amu
: 0.01 is too small for many applications.
See chapter 5 of Rossi et al for full discussion.
Argument Details
Data = list(lgtdata, Z, p)
[Z
optional]
lgtdata: | list of nlgt=length(lgtdata) lists with each cross-section unit MNL data |
lgtdata[[i]]$y: | n_i x 1 vector of multinomial outcomes (1, ..., p) |
lgtdata[[i]]$X: | n_i*p x nvar design matrix for i th unit |
Z: | nreg x nz matrix of unit chars (def: vector of ones) |
p: | number of choice alternatives |
Prior = list(a, deltabar, Ad, mubar, Amu, nu, V, a, ncomp, SignRes)
[all but ncomp
are optional]
a: | ncomp x 1 vector of Dirichlet prior parameters (def: rep(5,ncomp) ) |
deltabar: | nz*nvar x 1 vector of prior means (def: 0) |
Ad: | prior precision matrix for vec(D) (def: 0.01*I) |
mubar: | nvar x 1 prior mean vector for normal component mean (def: 0 if unrestricted; 2 if restricted) |
Amu: | prior precision for normal component mean (def: 0.01 if unrestricted; 0.1 if restricted) |
nu: | d.f. parameter for IW prior on normal component Sigma (def: nvar+3 if unrestricted; nvar+15 if restricted) |
V: | PDS location parameter for IW prior on normal component Sigma (def: nu*I if unrestricted; nu*D if restricted with d_pp = 4 if unrestricted and d_pp = 0.01 if restricted) |
ncomp: | number of components used in normal mixture |
SignRes: | nvar x 1 vector of sign restrictions on the coefficient estimates (def: rep(0,nvar) )
|
Mcmc = list(R, keep, nprint, s, w)
[only R
required]
R: | number of MCMC draws |
keep: | MCMC thinning parameter -- 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) |
s: | scaling parameter for RW Metropolis (def: 2.93/sqrt(nvar) ) |
w: | fractional likelihood weighting parameter (def: 0.1) |
Sign Restrictions
If \beta_ik
has a sign restriction: \beta_ik = SignRes[k] * exp(\beta*_ik)
To use sign restrictions on the coefficients, SignRes
must be an nvar x 1
vector containing values of either 0, -1, or 1. The value 0 means there is no sign restriction, -1 ensures that the coefficient is negative, and 1 ensures that the coefficient is positive. For example, if SignRes = c(0,1,-1)
, the first coefficient is unconstrained, the second will be positive, and the third will be negative.
The sign restriction is implemented such that if the the k
'th \beta
has a non-zero sign restriction (i.e., it is constrained), we have \beta_k = SignRes[k] * exp(\beta*_k)
.
The sign restrictions (if used) will be reflected in the betadraw
output. However, the unconstrained mixture components are available in nmix
. Important: Note that draws from nmix
are distributed according to the mixture of normals but not the coefficients in betadraw
.
Care should be taken when selecting priors on any sign restricted coefficients. See related vignette for additional information.
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:
Deltadraw |
|
betadraw |
|
nmix |
a list containing: |
loglike |
|
SignRes |
|
Note
Note: as of version 2.0-2 of bayesm
, the fractional weight parameter has been changed to a weight between 0 and 1.
w
is the fractional weight on the normalized pooled likelihood. This differs from what is in Rossi et al chapter 5, i.e.
like_i^{(1-w)} x like_pooled^{((n_i/N)*w)}
Large R
values may be required (>20,000).
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=10000} else {R=10}
set.seed(66)
p = 3 # num of choice alterns
ncoef = 3
nlgt = 300 # num of cross sectional units
nz = 2
Z = matrix(runif(nz*nlgt),ncol=nz)
Z = t(t(Z) - apply(Z,2,mean)) # demean Z
ncomp = 3 # num of mixture components
Delta = matrix(c(1,0,1,0,1,2),ncol=2)
comps=NULL
comps[[1]] = list(mu=c(0,-1,-2), rooti=diag(rep(1,3)))
comps[[2]] = list(mu=c(0,-1,-2)*2, rooti=diag(rep(1,3)))
comps[[3]] = list(mu=c(0,-1,-2)*4, rooti=diag(rep(1,3)))
pvec = c(0.4, 0.2, 0.4)
## simulate from MNL model conditional on X matrix
simmnlwX= function(n,X,beta) {
k = length(beta)
Xbeta = X%*%beta
j = nrow(Xbeta) / n
Xbeta = matrix(Xbeta, byrow=TRUE, ncol=j)
Prob = exp(Xbeta)
iota = c(rep(1,j))
denom = Prob%*%iota
Prob = Prob/as.vector(denom)
y = vector("double",n)
ind = 1:j
for (i in 1:n) {
yvec = rmultinom(1, 1, Prob[i,])
y[i] = ind%*%yvec
}
return(list(y=y, X=X, beta=beta, prob=Prob))
}
## simulate data
simlgtdata = NULL
ni = rep(50, 300)
for (i in 1:nlgt) {
betai = Delta%*%Z[i,] + as.vector(rmixture(1,pvec,comps)$x)
Xa = matrix(runif(ni[i]*p,min=-1.5,max=0), ncol=p)
X = createX(p, na=1, nd=NULL, Xa=Xa, Xd=NULL, base=1)
outa = simmnlwX(ni[i], X, betai)
simlgtdata[[i]] = list(y=outa$y, X=X, beta=betai)
}
## plot betas
if(0){
bmat = matrix(0, nlgt, ncoef)
for(i in 1:nlgt) {bmat[i,] = simlgtdata[[i]]$beta}
par(mfrow = c(ncoef,1))
for(i in 1:ncoef) { hist(bmat[,i], breaks=30, col="magenta") }
}
## set parms for priors and Z
Prior1 = list(ncomp=5)
keep = 5
Mcmc1 = list(R=R, keep=keep)
Data1 = list(p=p, lgtdata=simlgtdata, Z=Z)
## fit model without sign constraints
out1 = rhierMnlRwMixture(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)
}
## fit model with constraint that beta_i2 < 0 forall i
Prior2 = list(ncomp=5, SignRes=c(0,-1,0))
out2 = rhierMnlRwMixture(Data=Data1, Prior=Prior2, Mcmc=Mcmc1)
cat("Summary of Delta draws", fill=TRUE)
summary(out2$Deltadraw, tvalues=as.vector(Delta))
cat("Summary of Normal Mixture Distribution", fill=TRUE)
summary(out2$nmix)
## plotting examples
if(0) {
plot(out2$betadraw)
plot(out2$nmix)
}