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(p,lgtdata,Z) ( Z is optional)
`Prior` | list(deltabar,Ad,mubar,Amu,nu,V,ncomp) (all but ncomp are optional)
`Mcmc` | list(s,w,R,keep) (R required)

### Details

Model:
y_i ~ MNL(X_i,theta_i). i=1,..., length(lgtdata). theta_i is nvar x 1.

theta_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).

Priors:
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)

Lists contain:

`p`
  p is number of choice alternatives
`lgtdata`
  list of lists with each cross-section unit MNL data
`lgtdata[[i]]$y`
  n_i vector of multinomial outcomes (1,...,m)
`lgtdata[[i]]$X`
  n_i by nvar design matrix for ith unit
`deltabar`
  nz*nvar vector of prior means (def: 0)
`Ad`
  prior prec matrix for vec(D) (def: .01I)
`mubar`
  nvar x 1 prior mean vector for normal comp mean (def: 0)
`Amu`
  prior precision for normal comp mean (def: .01I)
`nu`
  d.f. parm for IW prior on norm comp Sigma (def: nvar+3)
`V`
  pds location parm for IW prior on norm comp Sigma (def: nuI)
`ncomp`
  number of components used in normal mixture
`s`
  scaling parm for RW Metropolis (def: 2.93/sqrt(nvar))
`w`
  fractional likelihood weighting parm (def: .1)
`R`
  number of MCMC draws
`keep`
  MCMC thinning parm: keep every keepth draw (def: 1)

### Value

a list containing:

`Deltadraw` | R/keep x nz*nvar matrix of draws of Delta, first row is initial value
`betadraw` | nlgt x nvar x R/keep array of draws of betas
`probdraw` | R/keep x ncomp matrix of draws of probs of mixture components (pvec)
`compdraw` | list of list of lists (length R/keep)
`loglike` | log-likelihood for each kept draw (length R/keep)

### Note

More on `compdraw` component of return value list:

compdraw[[i]]
  the ith draw of components for mixtures
compdraw[[i]][[j]]
  ith draw of the jth normal mixture comp
compdraw[[i]][[j]][[1]]
  ith draw of jth normal mixture comp mean vector
compdraw[[i]][[j]][[2]]
  ith draw of jth normal mixture cov parm (rooti)

Note: Z does not include an intercept and is centered for ease of interpretation.

Be careful in assessing prior parameter, Amu. .01 is too small for many applications. See Rossi et al, chapter 5 for full discussion.

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)*wi)

Large R values may be required (>20,000).

### Author(s)

Peter Rossi, Graduate School of Business, University of Chicago, Peter.Rossi@ChicagoGsb.edu.

### References

For further discussion, see Bayesian Statistics and Marketing by Rossi, Allenby and McCulloch, Chapter 5.
http://gsbwww.uchicago.edu/fac/peter.rossi/research/bsm.html

### See Also

`rmnlIndepMetrop`

### Examples

    ##
    if(nchar(Sys.getenv("LONG_TEST")) != 0)
    {

    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                                # no 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(.4,.2,.4)

    ## simulate data
    simlgtdata=NULL
    ni=rep(50,300)
    for (i in 1:nlgt)
    {  betai=Delta%*%Z[i,]+as.vector(rmixture(1,pvec,comps)$x)
       X=NULL
       for(j in 1:ni[i])
          { Xone=cbind(rbind(c(rep(0,p-1)),diag(p-1)),runif(p,min=-1.5,max=0))
            X=rbind(X,Xone) }
       outa=simmnlwX(ni[i],X,betai)
       simlgtdata[[i]]=list(y=outa$y,X=X,beta=betai)
    }

    ## plot betas
    if(0){
    ## set if(1) above to produce plots
    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
    nu=ncoef+3
    V=nu*diag(rep(1,ncoef))
    Ad=.01*(diag(rep(1,nz*ncoef)))
    mubar=matrix(rep(0,ncoef),nrow=1)
    deltabar=rep(0,ncoef*nz)
    Amu=matrix(.01,ncol=1)
    a=rep(5,ncoef)

    R=10000
    keep=5
    w=.1
    s=2.93/sqrt(ncoef)
    Data1=list(p=p,lgtdata=simlgtdata,Z=Z)
    Prior1=list(ncomp=ncomp,nu=nu,V=V,Amu=Amu,mubar=mubar,a=a,Ad=Ad,deltabar=deltabar)
    Mcmc1=list(s=s,w=w,R=R,keep=keep)
    out=rhierMnlRwMixture(Data=Data1,Prior=Prior1,Mcmc=Mcmc1)

    if(R < 1000) {begin=1} else {begin=1000/keep}
    end=R/keep
    cat(" Deltadraws ",fill=TRUE)
    mat=apply(out$Deltadraw[begin:end,],2,quantile,probs=c(.01,.05,.5,.95,.99))
    mat=rbind(as.vector(Delta),mat); rownames(mat)[1]="delta"; print(mat)

    tmom=momMix(matrix(pvec,nrow=1),list(comps))
    pmom=momMix(out$probdraw[begin:end,],out$compdraw[begin:end])
    mat=rbind(tmom$mu,pmom$mu)
    rownames(mat)=c("true","post expect")
    cat(" mu and posterior expectation of mu",fill=TRUE)
    print(mat)
    mat=rbind(tmom$sd,pmom$sd)
    rownames(mat)=c("true","post expect")
    cat(" std dev and posterior expectation of sd",fill=TRUE)
    print(mat)
    mat=rbind(as.vector(tmom$corr),as.vector(pmom$corr))
    rownames(mat)=c("true","post expect")
    cat(" corr and posterior expectation of corr",fill=TRUE)
    print(t(mat))

    if(0) {
    ## set if(1) to produce plots
    par(mfrow=c(4,1))
    plot(out$betadraw[1,1,])
    abline(h=simlgtdata[[1]]$beta[1])
    plot(out$betadraw[2,1,])
    abline(h=simlgtdata[[2]]$beta[1])
    plot(out$betadraw[100,3,])
    abline(h=simlgtdata[[100]]$beta[3])
    plot(out$betadraw[101,3,])
    abline(h=simlgtdata[[101]]$beta[3])
    par(mfrow=c(4,1))
    plot(out$Deltadraw[,1])
    abline(h=Delta[1,1])
    plot(out$Deltadraw[,2])
    abline(h=Delta[2,1])
    plot(out$Deltadraw[,3])
    abline(h=Delta[3,1])
    plot(out$Deltadraw[,6])
    abline(h=Delta[3,2])
    begin=1000/keep
    end=R/keep
    ngrid=50
    grid=matrix(0,ngrid,ncoef)
    rgm=matrix(c(-3,-7,-10,3,1,0),ncol=2)
    for(i in 1:ncoef) {rg=rgm[i,]; grid[,i]=rg[1] + ((1:ngrid)/ngrid)*(rg[2]-rg[1])}
    mden=eMixMargDen(grid,out$probdraw[begin:end,],out$compdraw[begin:end])
    par(mfrow=c(2,ncoef))
    for(i in 1:ncoef)
    {plot(grid[,i],mden[,i],type="l")}
    for(i in 1:ncoef)
    tden=mixDen(grid,pvec,comps)
    for(i in 1:ncoef)
    {plot(grid[,i],tden[,i],type="l")}
    }
    }

[Package bayesm version 2.0-8 Index]