Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

Re: st: bivariate ordered probit / ordinal probit - ctd.


From   Joseph Coveney <jcoveney@bigplanet.com>
To   Thomas Cornelißen <cornelissen@mbox.iqw.uni-hannover.de>, Statalist <statalist@hsphsun2.harvard.edu>
Subject   Re: st: bivariate ordered probit / ordinal probit - ctd.
Date   Mon, 07 Nov 2005 15:24:23 +0900

This fleshes out the first response a little.  To answer Thomas Cornelißen's
question regarding the discrepancy between coefficients fitted by -biprobit-
and those by -gllamm-:  the coefficients will differ by a scale factor of
sqrt(1 / (1 + random effect variance).  If Thomas multiplies the
coefficients obtained from -gllamm- by this scale factor, then results
from -biprobit- will match those from -gllamm-.  This is explained by
Leonardo Grilli and Carla Rampichini whose working paper's URL I gave in the
first response.

I've illustrated this in the do-file below.  The do-file shows the use of
both -xtprobit- (which is faster) and -gllamm-.  Thomas was interested in
rho, and -xtprobit- also gives rho directly, whereas with -gllamm- you need
to calculate rho from the outputted parameter estimate as described in the
first reply.  The illustration should generalize to ordered probit
regression with -reoprob- and -gllamm , link(oprobit)- using minus the first
cut point for the intercept (if I recall Stata's parameterization
correctly).

Note that -xtprobit- constrains rho to be positive.  If you have a negative
tetrachoric correlation coefficient for the two measures, then you can
use -gllamm- with an equation for the random effect as described in the
working paper by Leonardo Grilli & Carla Rampichini.  If you use -xtprobit-,
you'll need to flip the responses for one of them and note the change of
sign in the coeffients, including rho.

Joseph Coveney

set more off
drawnorm predictor latent0 latent1, ///
  corr(1 0.5 0.5 \ 0.5 1 0.5 \ 0.5 0.5 1) ///
  n(250) seed(`=date("2005-11-07", "ymd")') clear
forvalues i = 0/1 {
    generate byte manifest`i' = 0
    quietly replace manifest`i' = manifest`i' + (latent`i' > 0)
    drop latent`i'
}
generate int row = _n
biprobit (manifest0 predictor) (manifest1 predictor), nolog
quietly reshape long manifest, i(row) j(measure)
generate float interaction = measure * predictor
xtprobit manifest predictor measure interaction, ///
  i(row) intmethod(aghermite) intpoints(30) nolog
scalar scale_factor = sqrt( 1 / (1 + exp(_b[/lnsig2u])))
display _b[predictor] * scalar(scale_factor)
display _b[_cons] * scalar(scale_factor)
display (_b[predictor] + _b[interaction]) * scalar(scale_factor)
display (_b[_cons] + _b[measure]) * scalar(scale_factor)
gllamm manifest predictor measure interaction, ///
  i(row) family(binomial) link(probit) nip(30) adapt nolog
scalar scale_factor = sqrt( 1 / (1 + _b[row1:_cons]^2))
display _b[predictor] * scalar(scale_factor)
display _b[_cons] * scalar(scale_factor)
display (_b[predictor] + _b[interaction]) * scalar(scale_factor)
display (_b[_cons] + _b[measure]) * scalar(scale_factor)
display _b[row1:_cons]^2 / (1 + _b[row1:_cons]^2)
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/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index