rhierLinearModel {bayesm} | R Documentation

## Gibbs Sampler for Hierarchical Linear Model

### Description

`rhierLinearModel` implements a Gibbs Sampler for hierarchical linear models.

### Usage

    rhierLinearModel(Data, Prior, Mcmc)

### Arguments

`Data` | list(regdata,Z) (Z optional).
`Prior` | list(Deltabar,A,nu.e,ssq,nu,V) (optional).
`Mcmc` | list(R,keep) (R required).

### Details

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 ~ N(ZDelta[i,],V_{beta}).
Note: ZDelta is the matrix Z * Delta; [i,] refers to ith row of this product.

vec(Delta) given V_{beta} ~ N(vec(Deltabar),V_{beta} (x) A^{-1}).
V_{beta} ~ IW(nu,V).
Delta, Deltabar are nz x nvar. A is nz x nz. V_{beta} is nvar x nvar.

Note: if you don't have any z vars, set Z=iota (nreg x 1).

List arguments contain:

`regdata`
  list of lists with X,y matrices for each of length(regdata) regressions
`regdata[[i]]$X`
  X matrix for equation i
`regdata[[i]]$y`
  y vector for equation i
`Deltabar`
  nz x nvar matrix of prior means (def: 0)
`A`
  nz x nz matrix for prior precision (def: .01I)
`nu.e`
  d.f. parm for regression error variance prior (def: 3)
`ssq`
  scale parm for regression error var prior (def: var(y_i))
`nu`
  d.f. parm for Vbeta prior (def: nvar+3)
`V`
  Scale location matrix for Vbeta prior (def: nu*I)
`R`
  number of MCMC draws
`keep`
  MCMC thinning parm: keep every keepth draw (def: 1)

### Value

a list containing

`betadraw` | nreg x nvar x R/keep array of individual regression coef draws
`taudraw` | R/keep x nreg array of error variance draws
`Deltadraw` | R/keep x nz x nvar array of Deltadraws
`Vbetadraw` | R/keep x nvar*nvar array of Vbeta draws

### 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 3.
http://gsbwww.uchicago.edu/fac/peter.rossi/research/bsm.html

### Examples

    ##
    if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10}

    nreg=100; nobs=100; nvar=3
    Vbeta=matrix(c(1,.5,0,.5,2,.7,0,.7,1),ncol=3)
    Z=cbind(c(rep(1,nreg)),3*runif(nreg)); Z[,2]=Z[,2]-mean(Z[,2])
    nz=ncol(Z)
    Delta=matrix(c(1,-1,2,0,1,0),ncol=2)
    Delta=t(Delta) # first row of Delta is means of betas
    Beta=matrix(rnorm(nreg*nvar),nrow=nreg)%*%chol(Vbeta)+Z%*%Delta
    tau=.1
    iota=c(rep(1,nobs))
    regdata=NULL
    for (reg in 1:nreg) { X=cbind(iota,matrix(runif(nobs*(nvar-1)),ncol=(nvar-1)))
            y=X%*%Beta[reg,]+sqrt(tau)*rnorm(nobs); regdata[[reg]]=list(y=y,X=X) }

    nu.e=3
    ssq=NULL
    for(reg in 1:nreg) {ssq[reg]=var(regdata[[reg]]$y)}
    nu=nvar+3
    V=nu*diag(c(rep(1,nvar)))
    A=diag(c(rep(.01,nz)),ncol=nz)
    Deltabar=matrix(c(rep(0,nz*nvar)),nrow=nz)

    Data=list(regdata=regdata,Z=Z)
    Prior=list(Deltabar=Deltabar,A=A,nu.e=nu.e,ssq=ssq,nu=nu,V=V)
    Mcmc=list(R=R,keep=1)
    out=rhierLinearModel(Data=Data,Mcmc=Mcmc)

    cat(" Deltadraws ",fill=TRUE)
    mat=apply(out$Deltadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
    mat=rbind(as.vector(Delta),mat); rownames(mat)[1]="delta"; print(mat)
    cat(" Vbetadraws ",fill=TRUE)
    mat=apply(out$Vbetadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
    mat=rbind(as.vector(Vbeta),mat); rownames(mat)[1]="Vbeta"; print(mat)

[Package bayesm version 2.0-8 Index]