rmvpGibbs: Gibbs Sampler for Multivariate Probit

Description Usage Arguments Details Value Author(s) References See Also Examples

View source: R/rmvpgibbs_rcpp.r

Description

rmvpGibbs implements the Edwards/Allenby Gibbs Sampler for the multivariate probit model.

Usage

1
rmvpGibbs(Data, Prior, Mcmc)

Arguments

Data

list(y, X, p)

Prior

list(betabar, A, nu, V)

Mcmc

list(R, keep, nprint, beta0 ,sigma0)

Details

Model and Priors

w_i = X_iβ + e with e ~ N(0,Σ). Note: w_i is p x 1.
y_{ij} = 1 if w_{ij} > 0, else y_i = 0. j = 1, …, p

beta and Sigma are not identifed. Correlation matrix and the betas divided by the appropriate standard deviation are. See reference or example below for details.

β ~ N(betabar, A^{-1})
Σ ~ IW(nu, V)

To make X matrix use createX

Argument Details

Data = list(y, X, p)

X: n*p x k Design Matrix
y: n*p x 1 vector of 0/1 outcomes
p: dimension of multivariate probit

Prior = list(betabar, A, nu, V) [optional]

betabar: k x 1 prior mean (def: 0)
A: k x k prior precision matrix (def: 0.01*I)
nu: d.f. parameter for Inverted Wishart prior (def: (p-1)+3)
V: PDS location parameter for Inverted Wishart prior (def: nu*I)

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

R: number of MCMC draws
keep: MCMC thinning parameter -- keep every keepth draw (def: 1)
nprint: print the estimated time remaining for every nprint'th draw (def: 100, set to 0 for no print)
beta0: initial value for beta
sigma0: initial value for sigma

Value

A list containing:

betadraw

R/keep x k matrix of betadraws

sigmadraw

R/keep x p*p matrix of sigma draws – each row is the vector form of sigma

Author(s)

Peter Rossi, Anderson School, UCLA, [email protected].

References

For further discussion, see Chapter 4, Bayesian Statistics and Marketing by Rossi, Allenby, and McCulloch.
http://www.perossi.org/home/bsm-1

See Also

rmnpGibbs

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10}
set.seed(66)

simmvp = function(X, p, n, beta, sigma) {
  w = as.vector(crossprod(chol(sigma),matrix(rnorm(p*n),ncol=n))) + X%*%beta
  y = ifelse(w<0, 0, 1)
  return(list(y=y, X=X, beta=beta, sigma=sigma))
}

p = 3
n = 500
beta = c(-2,0,2)
Sigma = matrix(c(1, 0.5, 0.5, 0.5, 1, 0.5, 0.5, 0.5, 1), ncol=3)
k = length(beta)
I2 = diag(rep(1,p))
xadd = rbind(I2)
for(i in 2:n) { xadd=rbind(xadd,I2) }
X = xadd

simout = simmvp(X,p,500,beta,Sigma)

Data1 = list(p=p, y=simout$y, X=simout$X)
Mcmc1 = list(R=R, keep=1)

out = rmvpGibbs(Data=Data1, Mcmc=Mcmc1)

ind = seq(from=0, by=p, length=k)
inda = 1:3
ind = ind + inda
cat(" Betadraws ", fill=TRUE)
betatilde = out$betadraw / sqrt(out$sigmadraw[,ind])
attributes(betatilde)$class = "bayesm.mat"
summary(betatilde, tvalues=beta/sqrt(diag(Sigma)))

rdraw = matrix(double((R)*p*p), ncol=p*p)
rdraw = t(apply(out$sigmadraw, 1, nmat))
attributes(rdraw)$class = "bayesm.var"
tvalue = nmat(as.vector(Sigma))
dim(tvalue) = c(p,p)
tvalue = as.vector(tvalue[upper.tri(tvalue,diag=TRUE)])
cat(" Draws of Correlation Matrix ", fill=TRUE)
summary(rdraw, tvalues=tvalue)

## plotting examples
if(0){plot(betatilde, tvalues=beta/sqrt(diag(Sigma)))}

Example output

Table of y values
y
  0   1 
752 748 
 
Starting Gibbs Sampler for MVP
   500  obs of  3  binary indicators;  3  indep vars (including intercepts)
 
Prior Parms:
betabar
[1] 0 0 0
A
     [,1] [,2] [,3]
[1,] 0.01 0.00 0.00
[2,] 0.00 0.01 0.00
[3,] 0.00 0.00 0.01
nu
[1] 6
V
     [,1] [,2] [,3]
[1,]    6    0    0
[2,]    0    6    0
[3,]    0    0    6
 
MCMC Parms:
   10  reps; keeping every  1 th draw  nprint=  100
initial beta=  0 0 0
initial sigma= 
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1
 
 MCMC Iteration (est time to end - min) 
 Total Time Elapsed: 0.00 
 Betadraws 
fewer than 100 draws submitted 
 Draws of Correlation Matrix 
fewer than 100 draws submitted 

bayesm documentation built on July 21, 2017, 7:18 p.m.