MCMC Algorithm for Hierarchical Binary Logit
This function has been deprecated. Please use rhierMnlRwMixture instead.
rhierBinLogit implements an MCMC algorithm for hierarchical binary logits with a normal heterogeneity distribution. This is a hybrid sampler with a RW Metropolis step for unit-level logit parameters.
rhierBinLogit is designed for use on choice-based conjoint data with partial profiles. The Design matrix is based on differences of characteristics between two alternatives. See Appendix A of Bayesian Statistics and Marketing for details.
rhierBinLogit(Data, Prior, Mcmc)
Data |
list(lgtdata, Z) |
Prior |
list(Deltabar, ADelta, nu, V) |
Mcmc |
list(R, keep, sbeta) |
y_{hi} = 1 with Pr = exp(x_{hi}'β_h) / (1+exp(x_{hi}'β_h) and β_h is nvar x 1
h = 1, …, length(lgtdata) units (or "respondents" for survey data)
β_h = ZDelta[h,] + u_h
Note: here ZDelta refers to Z%*%Delta with ZDelta[h,] the hth row of this product
Delta is an nz x nvar array
u_h ~ N(0, V_{beta}).
delta = vec(Delta) ~ N(vec(Deltabar), V_{beta}(x) ADelta^{-1})
V_{beta} ~ IW(nu, V)
Data = list(lgtdata, Z) [Z optional]
lgtdata: |
list of lists with each cross-section unit MNL data |
lgtdata[[h]]$y: |
n_h x 1 vector of binary outcomes (0,1) |
lgtdata[[h]]$X: |
n_h x nvar design matrix for h'th unit |
Z: |
nreg x nz mat of unit chars (def: vector of ones) |
Prior = list(Deltabar, ADelta, nu, V) [optional]
Deltabar: |
nz x nvar matrix of prior means (def: 0) |
ADelta: |
prior precision matrix (def: 0.01I) |
nu: |
d.f. parameter for IW prior on normal component Sigma (def: nvar+3) |
V: |
pds location parm for IW prior on normal component Sigma (def: nuI) |
Mcmc = list(R, keep, sbeta) [only R required]
R: |
number of MCMC draws |
keep: |
MCMC thinning parm -- keep every keepth draw (def: 1) |
sbeta: |
scaling parm for RW Metropolis (def: 0.2) |
A list containing:
Deltadraw |
R/keep x nz*nvar matrix of draws of Delta |
betadraw |
nlgt x nvar x R/keep array of draws of betas |
Vbetadraw |
R/keep x nvar*nvar matrix of draws of Vbeta |
llike |
R/keep x 1 vector of log-like values |
reject |
R/keep x 1 vector of reject rates over nlgt units |
Some experimentation with the Metropolis scaling paramter (sbeta) may be required.
Peter Rossi, Anderson School, UCLA, perossichi@gmail.com.
For further discussion, see Chapter 5, Bayesian Statistics and Marketing by Rossi, Allenby, and McCulloch.
http://www.perossi.org/home/bsm-1
if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=10000} else {R=10}
set.seed(66)
nvar = 5 ## number of coefficients
nlgt = 1000 ## number of cross-sectional units
nobs = 10 ## number of observations per unit
nz = 2 ## number of regressors in mixing distribution
Z = matrix(c(rep(1,nlgt),runif(nlgt,min=-1,max=1)), nrow=nlgt, ncol=nz)
Delta = matrix(c(-2, -1, 0, 1, 2, -1, 1, -0.5, 0.5, 0), nrow=nz, ncol=nvar)
iota = matrix(1, nrow=nvar, ncol=1)
Vbeta = diag(nvar) + 0.5*iota%*%t(iota)
lgtdata=NULL
for (i in 1:nlgt) {
beta = t(Delta)%*%Z[i,] + as.vector(t(chol(Vbeta))%*%rnorm(nvar))
X = matrix(runif(nobs*nvar), nrow=nobs, ncol=nvar)
prob = exp(X%*%beta) / (1+exp(X%*%beta))
unif = runif(nobs, 0, 1)
y = ifelse(unif<prob, 1, 0)
lgtdata[[i]] = list(y=y, X=X, beta=beta)
}
Data1 = list(lgtdata=lgtdata, Z=Z)
Mcmc1 = list(R=R)
out = rhierBinLogit(Data=Data1, Mcmc=Mcmc1)
cat("Summary of Delta draws", fill=TRUE)
summary(out$Deltadraw, tvalues=as.vector(Delta))
cat("Summary of Vbeta draws", fill=TRUE)
summary(out$Vbetadraw, tvalues=as.vector(Vbeta[upper.tri(Vbeta,diag=TRUE)]))
if(0){
## plotting examples
plot(out$Deltadraw,tvalues=as.vector(Delta))
plot(out$betadraw)
plot(out$Vbetadraw,tvalues=as.vector(Vbeta[upper.tri(Vbeta,diag=TRUE)]))
}Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.