Dear WinBUGS users:
Last week, I asked for advice on how to implement a multivariate probit
model (with 3 equations).
I received many helpful answers. I pasted them below.
I truly thank all those who offered advice for their generous support and
their insightful suggestions. I am working my way through them.
Best regards,
giacomo
##############################################################
1.
You can place a Wishart prior on the precision matrix. This will give you
an estimate of the sigma matrix. You can then use the sigmas to scale the
regression coefficients to obtain probit coeffcients. This is consistent
with Chib and Greenberg who write: "Suppose that (gamma, Q) is an
alternative
parameterisation, where gamma is the regression parameter vector and Q is
the covariance matrix. Then it is easy to show that
pr(y|gamma, Q) = pr(y|beta, Sigma), where beta[j] = (q[jj])^-.5*gamma[j],
Sigma = CQC and C = diag{q[11]^-.5, ..., q[JJ]^-.5}.
From your code:
df <- 3
omega[1:3,1:3] ~ dwish(R[ , ], df)
sigma[1:3,1:3] <- inverse(omega[,])
for (j in 1:3) {
for (k in 1:3) {
bstar[j,k] <- b[j,k]/sqrt(sigma[j,j])
}
}
################################################################
2.
Isn't the matrix guaranteed to be positive definite? I guess the
problem is that you are getting correlations that are very close to 1
or -1, and then the numerical routine used to invert the matrix fails.
Have you tried to use a more informative prior for your
correlations? Maybe start with a Beta(20,20) and see if it runs.
And then start decreasing the parameters of the beta. Another option
might be to try to hard code the inversion of the correlation matrix,
if you are only interested in the trivariate probit model.
#################################################################
3.
This might help....
van den Berg SM, Setiawan A, Bartels M, Polderman TJ, van der Vaart AW,
Boomsma
DI. 2006 Individual differences in puberty onset in girls: Bayesian
estimation
of heritabilities and genetic correlations.
Behav Genet. 2006 Mar;36(2):261-70.
#####################################################################
4.
I may have missed something but can you orthogonalize (Cholesky?) to
three univariate normals and then transform back? That's typically how
you generate samples from a MVN in any case, I think, so BUGs should like
it.
It shouldn't be too bad for a 3x3.
#######################################################################
5.
Practically, probit model with full variance-covariance matrix is
normalised by
setting one of the diagonal elements in the covariance matrix of error
differences to 1. To my understanding, your code seemed to set all the
diagonal
elements of the covariance matrix to 1.
Try Robert E. McCulloch, Nicholas G. Polson and Peter E. Rossi (2000) A
Bayesian analysis of the multinomial probit model with fully identified
parameters, Journal of Econometrics, Vol.99 (1), 173-193.
and the follow-up comment and response on Vol.99 (2)
########################################################################
6.
The problem is that you have allowed the covariance to not have row
diagonal
dominance. An easy fix is to let the covariances range form -.499999 to
.499999. If the sum of the absolute values of the off diagonal elements of
a
row are less than the absolute value of the diagonal element, then the
matrix is
non-singular. Nonsingular covariance matirces should have this property. I
think it is a similar criterion to the Cauchy Swartz rule that forces the
correlation between -1 and 1.
#########################################################################
7.
One idea,
Let x1, x2, x3 be three i.i.d. N(0,1) random variables,
Now let a1, a2, a3 be three vectors with three elements. Also let a1, a2
and
a3 lie on a unit sphere, then
x1*a1,x2*a2,x3*a3 are three MVN random vectors of length three with unit
variances. These vectors span the space of all possible MVN random vectors
with unit variances.
Now, without loss of generality the covariance matrix is orthogonal, so
that
any rigid rotation applied equally to all a1, a2 and a3 would give the
same
transformations. So we transform a1 = a11,0,0, a2= a21,a22,0 and
a3=a31,a32,a33 as the transformation has at least two dimensions.
Since the rotation is rigid, the new vectors also lie on the unit sphere.
The rotation is as follows. First pick a vector and rotate the frame so
that
this vector is just 1,0,0. Next rotate the frame around this vector so the
second vector lies in the x,y plane, the we can express the three vectors
(at least one way), using say spherical coordinates.
Using spherical coordinates, let
a11=1,
a21=cos(c1), a22= sin(c1),
a31=cos(c2)*cos(c3), a32=sin(c2)*cos(c3) and a33=sin(c3)
e.g.
a1=c(1,0,0)
a2=c(cos(c1),sin(c1),0)
a3=c(cos(c2)*cos(c3), sin(c2)*cos(c3),sin(c3))
The correlations are:
cor(ai,ai)=1
cor(a1,a2)=cos(c1)
cor(a1,a3)=cos(c2)*cos(c3)
cor(a2,a3)=cos(c1)*cos(c2)*cos(c3)+sin(c1)*sin(c2)*cos(c3)
So I would try
pi<-3.1415926
c1~dunif(-pi,pi)
c2~dunif(-pi,pi)
c30~dflat(-1,1)
c3<-asin(c30) ###This gives uniform distribution of a3 on sphere.
rho[1]<-cos(c1)
rho[2]<-cos(c2)*cos(c3)
rho[3]<-cos(c1)*cos(c2)*cos(c3)+sin(c1)*sin(c2)*cos(c3)
I hope this helps. I have not had time to try it out.
Since the rotations are not unique, you can have other rho values.
This might do for a start.
-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask]
To mail the BUGS list, mail to [log in to unmask]
Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html
Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list
|