rhierNegbinRw {bayesm} | R Documentation

## MCMC Algorithm for Negative Binomial Regression

### Description

`rhierNegbinRw` implements an MCMC strategy for the hierarchical Negative Binomial (NBD) regression model. Metropolis steps for each unit level set of regression parameters are automatically tuned by optimization. Over-dispersion parameter (alpha) is common across units.

### Usage

    rhierNegbinRw(Data, Prior, Mcmc)

### Arguments

`Data` | list(regdata,Z)
`Prior` | list(Deltabar,Adelta,nu,V,a,b)
`Mcmc` | list(R,keep,s_beta,s_alpha,c,Vbeta0,Delta0)

### Details

Model: y_i ~ NBD(mean=lambda, over-dispersion=alpha).
lambda=exp(X_ibeta_i)

Prior: beta_i ~ N(Delta'z_i,Vbeta).

vec(Delta|Vbeta) ~ N(vec(Deltabar),Vbeta (x) Adelta).
Vbeta ~ IW(nu,V).
alpha ~ Gamma(a,b).
note: prior mean of alpha = a/b, variance = a/(b^2)

list arguments contain:

`regdata`
  list of lists with data on each of nreg units
`regdata[[i]]$X`
  nobs_i x nvar matrix of X variables
`regdata[[i]]$y`
  nobs_i x 1 vector of count responses
`Z`
  nreg x nz mat of unit chars (def: vector of ones)
`Deltabar`
  nz x nvar prior mean matrix (def: 0)
`Adelta`
  nz x nz pds prior prec matrix (def: .01I)
`nu`
  d.f. parm for IWishart (def: nvar+3)
`V`
  location matrix of IWishart prior (def: nuI)
`a`
  Gamma prior parm (def: .5)
`b`
  Gamma prior parm (def: .1)
`R`
  number of MCMC draws
`keep`
  MCMC thinning parm: keep every keepth draw (def: 1)
`s\_beta`
  scaling for beta| alpha RW inc cov (def: 2.93/sqrt(nvar))
`s\_alpha`
  scaling for alpha | beta RW inc cov (def: 2.93)
`c`
  fractional likelihood weighting parm (def:2)
`Vbeta0`
  starting value for Vbeta (def: I)
`Delta0`
  starting value for Delta (def: 0)

### Value

a list containing:

`llike` | R/keep vector of values of log-likelihood
`betadraw` | nreg x nvar x R/keep array of beta draws
`alphadraw` | R/keep vector of alpha draws
`acceptrbeta` | acceptance rate of the beta draws
`acceptralpha` | acceptance rate of the alpha draws

### Note

The NBD regression encompasses Poisson regression in the sense that as alpha goes to infinity the NBD distribution tends to the Poisson.
For "small" values of alpha, the dependent variable can be extremely variable so that a large number of observations may be required to obtain precise inferences.

For ease of interpretation, we recommend demeaning Z variables.

### Author(s)

Sridhar Narayanam & 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

`rnegbinRw`

### Examples

    ##
    if(nchar(Sys.getenv("LONG_TEST")) != 0)
    {
    ##
    set.seed(66)
    simnegbin =
    function(X, beta, alpha) {
    #   Simulate from the Negative Binomial Regression
    lambda = exp(X %*% beta)
    y=NULL
    for (j in 1:length(lambda))
        y = c(y,rnbinom(1,mu = lambda[j],size = alpha))
    return(y)
    }

    nreg = 100        # Number of cross sectional units
    T = 50            # Number of observations per unit
    nobs = nreg*T
    nvar=2            # Number of X variables
    nz=2              # Number of Z variables

    # Construct the Z matrix
    Z = cbind(rep(1,nreg),rnorm(nreg,mean=1,sd=0.125))

    Delta = cbind(c(0.4,0.2), c(0.1,0.05))
    alpha = 5
    Vbeta = rbind(c(0.1,0),c(0,0.1))

    # Construct the regdata (containing X)
    simnegbindata = NULL
    for (i in 1:nreg) {
        betai = as.vector(Z[i,]%*%Delta) + chol(Vbeta)%*%rnorm(nvar)
        X = cbind(rep(1,T),rnorm(T,mean=2,sd=0.25))
        simnegbindata[[i]] = list(y=simnegbin(X,betai,alpha), X=X,beta=betai)
    }

    Beta = NULL
    for (i in 1:nreg) {Beta=rbind(Beta,matrix(simnegbindata[[i]]$beta,nrow=1))}

    Data = list(regdata=simnegbindata, Z=Z)
    Deltabar = matrix(rep(0,nvar*nz),nrow=nz)
    Vdelta = 0.01 * diag(nvar)
    nu = nvar+3
    V = 0.01*diag(nvar)
    a = 0.5
    b = 0.1
    Prior = list(Deltabar=Deltabar, Vdelta=Vdelta, nu=nu, V=V, a=a, b=b)

    R=10000
    keep =1
    s_beta=2.93/sqrt(nvar)
    s_alpha=2.93
    c=2
    Mcmc = list(R=R, keep = keep, s_beta=s_beta, s_alpha=s_alpha, c=c)
    out = rhierNegbinRw(Data, Prior, Mcmc)

    # Unit level mean beta parameters
    Mbeta = matrix(rep(0,nreg*nvar),nrow=nreg)
    ndraws = length(out$alphadraw)
    for (i in 1:nreg) { Mbeta[i,] = rowSums(out$Betadraw[i, , ])/ndraws }

    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)
    cat(" alphadraws ",fill=TRUE)
    mat=apply(matrix(out$alphadraw),2,quantile,probs=c(.01,.05,.5,.95,.99))
    mat=rbind(as.vector(alpha),mat); rownames(mat)[1]="alpha"; print(mat)
    }

[Package bayesm version 2.0-8 Index]