rordprobitGibbs function

Gibbs Sampler for Ordered Probit

rordprobitGibbs implements a Gibbs Sampler for the ordered probit model with a RW Metropolis step for the cut-offs.

rordprobitGibbs(Data, Prior, Mcmc)

Arguments

• Data: list(y, X, k)
• Prior: list(betabar, A, dstarbar, Ad)
• Mcmc: list(R, keep, nprint, s)

Details

Model and Priors

$z = X\beta + e$ with $e$ $~$ $N(0, I)$

$y = k$ if c[k] $\le z <$ c[k+1] with $k = 1,\ldots,K$

cutoffs = {c[1], $\ldots$, c[k+1]}

$\beta$ $~$ $N(betabar, A^{-1})$

$dstar$ $~$ $N(dstarbar, Ad^{-1})$

Be careful in assessing prior parameter Ad: 0.1 is too small for many applications.

Argument Details

Data = list(y, X, k)

y:	$n x 1$ vector of observations, ( $1, \ldots, k$ )
X:	$n x p$ Design Matrix
k:	the largest possible value of y

Prior = list(betabar, A, dstarbar, Ad) [optional]

betabar:	$p x 1$ prior mean (def: 0)
A:	$p x p$ prior precision matrix (def: 0.01*I)
dstarbar:	$ndstar x 1$ prior mean, where $ndstar=k-2$ (def: 0)
Ad:	$ndstar x ndstar$ prior precision matrix (def: I)

Mcmc = list(R, keep, nprint, s) [only R required]

R:	number of MCMC draws
keep:	MCMC thinning parameter -- keep every keep th draw (def: 1)
nprint:	print the estimated time remaining for every nprint 'th draw (def: 100, set to 0 for no print)
s:	scaling parameter for RW Metropolis (def: 2.93/ sqrt(p) )

Returns

A list containing: - betadraw: $R/keep x p$ matrix of betadraws

• cutdraw: $R/keep x (k-1)$ matrix of cutdraws

• dstardraw: $R/keep x (k-2)$ matrix of dstardraws

• accept: acceptance rate of Metropolis draws for cut-offs

Note

set c[1] = -100 and c[k+1] = 100. c[2] is set to 0 for identification.

The relationship between cut-offs and dstar is:

c[3] = exp(dstar[1]),

c[4] = c[3] + exp(dstar[2]), ...,

c[k] = c[k-1] + exp(dstar[k-2])

References

Bayesian Statistics and Marketing by Rossi, Allenby, and McCulloch.

Author(s)

Peter Rossi, Anderson School, UCLA, perossichi@gmail.com .

See Also

rbprobitGibbs

Examples

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

## simulate data for ordered probit model

simordprobit=function(X, betas, cutoff){
  z = X%*%betas + rnorm(nobs)   
  y = cut(z, br = cutoff, right=TRUE, include.lowest = TRUE, labels = FALSE)  
  return(list(y = y, X = X, k=(length(cutoff)-1), betas= betas, cutoff=cutoff ))
}

nobs = 300 
X = cbind(rep(1,nobs),runif(nobs, min=0, max=5),runif(nobs,min=0, max=5))
k = 5
betas = c(0.5, 1, -0.5)       
cutoff = c(-100, 0, 1.0, 1.8, 3.2, 100)
simout = simordprobit(X, betas, cutoff) 
  
Data=list(X=simout$X, y=simout$y, k=k)

## set Mcmc for ordered probit model
   
Mcmc = list(R=R)   
out = rordprobitGibbs(Data=Data, Mcmc=Mcmc)

cat(" ", fill=TRUE)
cat("acceptance rate= ", accept=out$accept, fill=TRUE)
 
## outputs of betadraw and cut-off draws
  
cat(" Summary of betadraws", fill=TRUE)
summary(out$betadraw, tvalues=betas)
cat(" Summary of cut-off draws", fill=TRUE) 
summary(out$cutdraw, tvalues=cutoff[2:k])

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