rhierMnlDP              package:bayesm              R Documentation

_M_C_M_C _A_l_g_o_r_i_t_h_m _f_o_r _H_i_e_r_a_r_c_h_i_c_a_l _M_u_l_t_i_n_o_m_i_a_l _L_o_g_i_t _w_i_t_h _D_i_r_i_c_h_l_e_t _P_r_o_c_e_s_s _P_r_i_o_r _H_e_t_e_r_o_g_e_n_e_i_t_y

_D_e_s_c_r_i_p_t_i_o_n:

     'rhierMnlDP' is a MCMC algorithm for a hierarchical multinomial
     logit with a Dirichlet Process Prior for the distribution of
     heteorogeneity.  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 there are panel units.  This is a hybrid Gibbs
     Sampler with a RW Metropolis step for the MNL  coefficients for
     each panel unit.  This procedure can be interpreted as a Bayesian
     semi-parameteric method in the sense that the DP prior can
     accomodate heterogeniety of an unknown form.

_U_s_a_g_e:

     rhierMnlDP(Data, Prior, Mcmc)

_A_r_g_u_m_e_n_t_s:

    Data: list(p,lgtdata,Z) ( Z is optional) 

   Prior: list(deltabar,Ad,Prioralpha,lambda_hyper) (all are optional)

    Mcmc: list(s,w,R,keep) (R required)

_D_e_t_a_i_l_s:

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

     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. 

     beta_i ~ N(mu_i,Sigma_i). 

     Priors: 
      theta_i=(mu_i,Sigma_i) ~ DP(G_0(lambda),alpha)
      G_0(lambda):
      mu_i | Sigma_i ~ N(0,Sigma_i (x) a^{-1})
      Sigma_i ~ IW(nu,nu*v*I)

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

     alpha ~ (1-(alpha-alphamin)/(alphamax-alphamin))^power 
      alpha= alphamin then expected number of components = Istarmin 
      alpha= alphamax then expected number of components = Istarmax 

     Lists contain: 

     Data:

   '_p' p is number of choice alternatives

   '_l_g_t_d_a_t_a' list of lists with each cross-section unit MNL data

   '_l_g_t_d_a_t_a[[_i]]$_y' n_i vector of multinomial outcomes (1,...,m)

   '_l_g_t_d_a_t_a[[_i]]$_X' n_i by nvar design matrix for ith unit

     Prior: 

   '_d_e_l_t_a_b_a_r' nz*nvar vector of prior means (def: 0)

   '_A_d' prior prec matrix for vec(D) (def: .01I)

     Prioralpha:

   '_I_s_t_a_r_m_i_n' expected number of components at lower bound of support
        of alpha def(1)

   '_I_s_t_a_r_m_a_x' expected number of components at upper bound of support
        of alpha (def: min(50,.1*nlgt))

   '_p_o_w_e_r' power parameter for alpha prior (def: .8)

     lambda_hyper:

   '_a_l_i_m' defines support of a distribution,def:c(.01,2) 

   '_n_u_l_i_m' defines support of nu distribution, def:c(.01,3) 

   '_v_l_i_m' defines support of v distribution, def:c(.1,4) 

     Mcmc:

   '_R' number of mcmc draws

   '_k_e_e_p' thinning parm, keep every keepth draw

   '_m_a_x_u_n_i_q' storage constraint on the number of unique components

   '_g_r_i_d_s_i_z_e' number of discrete points for hyperparameter priors,def:
        20

_V_a_l_u_e:

     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

    nmix: list of 3 components, probdraw, NULL, compdraw 

   adraw: R/keep draws of hyperparm a

   vdraw: R/keep draws of hyperparm v

  nudraw: R/keep draws of hyperparm nu

Istardraw: R/keep draws of number of unique components

alphadraw: R/keep draws of number of DP tightness parameter

 loglike: R/keep draws of log-likelihood

_N_o_t_e:

     As is well known, Bayesian density estimation involves computing
     the predictive distribution of a "new" unit parameter, theta_{n+1}
     (here "n"=nlgt). This is done by averaging the normal base
     distribution over draws from the distribution of theta_{n+1} given
     theta_1, ..., theta_n,alpha,lambda,Data. To facilitate this, we
     store those draws from the predictive distribution of theta_{n+1}
     in a list structure compatible  with other 'bayesm' routines that
     implement a finite mixture of normals.

     More on nmix list:
       contains the draws from the predictive distribution of a "new"
     observations parameters.  These are simply the parameters of one
     normal distribution.  We enforce compatibility with a mixture of k
     components in order to utilize generic summary  plotting
     functions.  

     Therefore,'probdraw' is a vector of ones.  'zdraw' (indicator
     draws) is omitted as it is not necessary for density estimation.
     'compdraw' contains the draws of the theta_{n+1} as a list of list
     of lists.

     More on 'compdraw' component of return value list:

   _c_o_m_p_d_r_a_w[[_i]] ith draw of components for mixtures

   _c_o_m_p_d_r_a_w[[_i]][[_1]] ith draw of the thetanp1

   _c_o_m_p_d_r_a_w[[_i]][[_1]][[_1]] ith draw of mean vector

   _c_o_m_p_d_r_a_w[[_i]][[_1]][[_2]] ith draw of parm (rooti)

     We parameterize the prior on Sigma_i 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.

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

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

_A_u_t_h_o_r(_s):

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

_R_e_f_e_r_e_n_c_e_s:

     For further discussion, see _Bayesian Statistics and Marketing_ by
     Rossi, Allenby and McCulloch, Chapter 5. 
      <URL:
     http://faculty.chicagogsb.edu/peter.rossi/research/bsm.html>

_S_e_e _A_l_s_o:

     'rhierMnlRwMixture'

_E_x_a_m_p_l_e_s:

     ##
     if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=20000} 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                                # 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(2,3)))
     comps[[2]]=list(mu=c(0,-1,-2)*2,rooti=diag(rep(2,3)))
     comps[[3]]=list(mu=c(0,-1,-2)*4,rooti=diag(rep(2,3)))
     pvec=c(.4,.2,.4)

     simmnlwX= function(n,X,beta) {
       ##  simulate from MNL model conditional on X matrix
       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 with a mixture of 3 normals
     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(1){
     ## 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 Data and Mcmc lists
     keep=5
     Mcmc1=list(R=R,keep=keep)
     Data1=list(p=p,lgtdata=simlgtdata,Z=Z)

     out=rhierMnlDP(Data=Data1,Mcmc=Mcmc1)

     cat("Summary of Delta draws",fill=TRUE)
     summary(out$Deltadraw,tvalues=as.vector(Delta))

     if(0) {
     ## plotting examples
     plot(out$betadraw)
     plot(out$nmix)
     }

