st: bivariate ordered probit / ordinal probit - ctd.

[Thomas Corneli�en <cornelissen@mbox.iqw.uni-hannover.de>]

Dear Statalist,
I am referring to last year's threads "bivariate ordinal probit model" and 
"[multivariate ordered probit models in Stata]".

Joseph Coveney suggested to use -gllamm- to estimate a bivariate ordered 
probit model and gave an example do-file.

Before applying -gllamm- to a bivariate ordered probit I wanted to check 
whether I was using it correctly and therefore wanted to estimate a 
bivariate probit by -gllamm- and compared it to Stata's results 
form -biprobit-.

Unfortunately, I do not get similar estimates, although the log likelihood 
at convergence is similar.

I feel that the difference might be due to the fact that -gllamm- estimates 
individual random effects while biprobit does not, but I am not sure whether 
I understand that correctly.
Moreover, -gllamm- does not seem to give me an estimate of the parameter 
'rho' estimated by the bivariate probit model?

The results and the code are below.

Are there any suggestions about how to have the same estimates from -gllamm- 
as from -biprobit-, so that I could then transfer it to a bivariate ordered 
probit ?
Thank for any comments!
Thomas
-------------------------------------------------------------------------
Thomas Cornelissen
Institute of Empirical Economic Research
University of Hannover, Germany
cornelissen@mbox.iqw.uni-hannover.de


+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-biprobit- delivers:

Iteration 4:   log pseudolikelihood = -213.50427

Seemingly unrelated bivariate probit              Number of obs   = 
250
                                                 Wald chi2(2)    = 
120.85
Log pseudolikelihood = -213.50427                 Prob > chi2     = 
0.0000

                                 (Std. Err. adjusted for 250 clusters in 
pid)
------------------------------------------------------------------------------
            |               Robust
            |      Coef.   Std. Err.      z    P>|z|     [95% Conf. 
Interval]
-------------+----------------------------------------------------------------
y0           |
         x1 |   1.010417     .11046     9.15   0.000     .7939194 
1.226915
      _cons |    .083602   .0920831     0.91   0.364    -.0968775 
.2640815
-------------+----------------------------------------------------------------
y1           |
         x1 |   1.189647   .1300007     9.15   0.000     .9348505 
1.444444
      _cons |  -.0102275   .0950683    -0.11   0.914    -.1965579 
.1761028
-------------+----------------------------------------------------------------
    /athrho |   .9527857   .1551993     6.14   0.000     .6486006 
1.256971
-------------+----------------------------------------------------------------
        rho |   .7410416   .0699728                      .5707272 
.8502268
------------------------------------------------------------------------------
Wald test of rho=0:                 chi2(1) =  37.6887    Prob > chi2 = 
0.0000


++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-gllamm- delivers:
(x2 is the dummy to identify the equations and iac is the interaction of x1 
and x2)

gllamm model

log likelihood = -213.50427


Variances and covariances of random effects
------------------------------------------------------------------------------

***level 2 (pid)

   var(1): 2.8616262 (1.0396932)
------------------------------------------------------------------------------

------------------------------------------------------------------------------
          y |      Coef.   Std. Err.      z    P>|z|     [95% Conf. 
Interval]
-------------+----------------------------------------------------------------
y            |
         x1 |   1.985573   .3561421     5.58   0.000     1.287547 
2.683598
         x2 |  -.1843845   .1846338    -1.00   0.318      -.54626 
.177491
        iac |   .3522056    .288057     1.22   0.221    -.2123757 
.916787
      _cons |   .1642864   .1813658     0.91   0.365    -.1911839 
.5197568
-------------+----------------------------------------------------------------
pid1         |
      _cons |   1.691634   .3073044     5.50   0.000     1.089329 
2.29394
------------------------------------------------------------------------------

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
And here is the code, a modification of Joseph's suggestion from the earlier 
threads::

set more off
drawnorm res2 res3, ///
corr(1 0.7 \ 0.7 1) ///
n(250) seed(20040318) clear
generate float x1 = invnorm(uniform())
forvalues dep = 2/3 {
generate float y`dep' = 0.7 * x1 + ///
sqrt(1 - 0.7^2) * res`dep'
}
generate byte y0 = 0
replace y0 = y0 + (norm(y2) > 0.5)
generate byte y1 = 0
replace y1 = y1 + (norm(y3) > 0.5)

drop y2 y3 res*
generate int pid = _n
*
biprobit (y0 x1) (y1 x1), cluster(pid)
gen i=1
keep x1 y* pid
reshape long y, i(pid) j(x2)
generate float iac = x1 * x2
gllamm y x1 x2 iac, i(pid) family(binomial) link(probit) nip(30) adapt trace 
allc
exit 
*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/