rordprobitGibbs {bayesm} R Documentation

## Gibbs Sampler for Ordered Probit

### Description

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

### Usage

`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β + e with e ~ N(0, I)
y = k if c[k] ≤ z < c[k+1] with k = 1,…,K
cutoffs = {c, , c[k+1]}

β ~ N(betabar, A^{-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, …, 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)`)

### Value

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 = -100 and c[k+1] = 100. c is set to 0 for identification.

The relationship between cut-offs and dstar is:
c = exp(dstar),
c = c + exp(dstar), ...,
c[k] = c[k-1] + exp(dstar[k-2])

### Author(s)

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

### References

Bayesian Statistics and Marketing by Rossi, Allenby, and McCulloch
http://www.perossi.org/home/bsm-1

`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