rhierLinearMixture          package:bayesm          R Documentation

_G_i_b_b_s _S_a_m_p_l_e_r _f_o_r _H_i_e_r_a_r_c_h_i_c_a_l _L_i_n_e_a_r _M_o_d_e_l

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

     'rhierLinearMixture' implements a Gibbs Sampler for hierarchical
     linear models with a mixture of normals prior.

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

     rhierLinearMixture(Data, Prior, Mcmc)

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

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

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

     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:

   '_r_e_g_d_a_t_a' list of lists with X,y matrices for each of
        length(regdata) regressions

   '_r_e_g_d_a_t_a[[_i]]$_X' X matrix for equation i 

   '_r_e_g_d_a_t_a[[_i]]$_y' y vector for equation i 

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

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

   '_m_u_b_a_r' nvar x 1 prior mean vector for normal comp mean (def: 0)

   '_A_m_u' prior precision for normal comp mean (def: .01I)

   '_n_u' 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)

   '_n_u._e' d.f. parm for regression error variance prior (def: 3)

   '_s_s_q' scale parm for regression error var prior (def: var(y_i))

   '_n_c_o_m_p' number of components used in normal mixture 

   '_R' number of MCMC draws

   '_k_e_e_p' MCMC thinning parm: keep every keepth draw (def: 1)

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

     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)

_N_o_t_e:

     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: 

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

_c_o_m_p_d_r_a_w[[_i]][[_j]] ith draw of the jth normal mixture comp

_c_o_m_p_d_r_a_w[[_i]][[_j]][[_1]] ith draw of jth normal mixture comp mean vector

_c_o_m_p_d_r_a_w[[_i]][[_j]][[_2]] ith draw of jth normal mixture cov parm (rooti) 

     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.

_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 3. 
      <URL:
     http://faculty.chicagogsb.edu/peter.rossi/research/bsm.html>

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

     'rhierLinearModel'

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

     ##
     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)
     }

