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/
