Email discussion lists for the UK Education and Research communities

## BUGS@JISCMAIL.AC.UK

#### View:

 Message: [ First | Previous | Next | Last ] By Topic: [ First | Previous | Next | Last ] By Author: [ First | Previous | Next | Last ] Font: Proportional Font

#### Options

Subject:

Summary: multivariate probit model

From:

Date:

Wed, 28 Feb 2007 01:37:44 -0500

Content-Type:

TEXT/PLAIN

Parts/Attachments:

 TEXT/PLAIN (185 lines)
 ```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 ```