st: gllamm bivariate probit randome effects

["Frank Muehsam" <f.muehsam@web.de>]

Dear Statalist,

I habe a question on estimating a bivariate probit with random effects and
 would appreciate comments and thoughts. I want to use Gllamm for this model.
I am still trying to understand the "language" of Gllamm to which I'm not used to. 
This question came up a while ago, but was not solved as fas as I can see:
http://www.stata.com/statalist/archive/2010-03/msg00171.html

In an earlier discussion, Joseph Coveney posted an example, how biprobit, 
xtprobit and gllamm can be used to estimate a bivariate probit:
http://www.stata.com/statalist/archive/2005-11/msg00180.html

I repeat the example here (and add  "version 9.0"): 

* Example 1 start*/
version 9.0
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
* Example 1 end

Ok, this works great and I wanted to replicate the example with "real" data but 
it did not work and I would be happy to know why:

* Example 2 start 
clear
use http://www.stata-press.com/data/r11/school, clear
eststo e1:biprobit private vote years logptax loginc
expand 2
bys obs: gen rid = _n-1
gen yvar = private if rid == 0
replace yvar = vote if rid == 1


generate Xyears = years * rid
generate Xlogptax = logptax * rid
generate Xloginc = loginc * rid

xtprobit yvar years logptax loginc rid Xyears Xlogptax Xloginc, ///
  i(obs) intmethod(aghermite) intpoints(30) nolog

scalar scale_factor = sqrt( 1 / (1 + exp(_b[/lnsig2u])))
display _b[years] * scalar(scale_factor)
display _b[logptax] * scalar(scale_factor)
display _b[loginc] * scalar(scale_factor)
display _b[_cons] * scalar(scale_factor)
display (_b[Xyears]+_b[years]) * scalar(scale_factor)
display (_b[Xlogptax]+_b[logptax]) * scalar(scale_factor)
display (_b[Xloginc]+_b[loginc]) * scalar(scale_factor)
display (_b[rid]+_b[_cons]) * scalar(scale_factor)


gllamm yvar years logptax loginc rid Xyears Xlogptax Xloginc, ///
  i(obs) family(binomial) link(probit) nip(30) adapt 

scalar scale_factor = sqrt( 1 / (1 + _b[obs1:_cons]^2))  
display _b[years] * scalar(scale_factor)
display _b[logptax] * scalar(scale_factor)
display _b[loginc] * scalar(scale_factor)
display _b[_cons] * scalar(scale_factor)
display (_b[Xyears]+_b[years]) * scalar(scale_factor)
display (_b[Xlogptax]+_b[logptax]) * scalar(scale_factor)
display (_b[Xloginc]+_b[loginc]) * scalar(scale_factor)
display (_b[rid]+_b[_cons]) * scalar(scale_factor)  
display _b[obs1:_cons]^2 / (1 + _b[obs1:_cons]^2)
* Example 2 end

The results of gllamm and xtprobit look similar but not equal and differ from 
biprobit. Did I construct the reshaped data structure correct or is something 
with the scaling parameter wrong? 

The model I would like to estimate has two additional random effects for each 
equation that are (allowed to be) correlated and bivariate normal distributed. 
So I would like to estimate four components of the error: standard deviaton 
of the random effects, their correlation and the correlation of the error term. 

Then I thought, I should use gllamms "eq" command and specify two random 
intercepts (I hope, I get the syntax correct...). I specified the model similar 
to what is in the great review by Grilli and Ramipichini:
http://www.cmm.bristol.ac.uk/learning-training/multilevel-m-software/reviewgllamm.pdf

*  Example 3 start
clear
use http://www.stata-press.com/data/r11/school, clear
eststo e1:biprobit private vote years logptax loginc
set seed 2728
expand 2
bys obs: gen rid = _n-1
gen yvar = private if rid == 0
replace yvar = vote if rid == 1

generate Xyears = years * rid
generate Xlogptax = logptax * rid
generate Xloginc = loginc * rid

gen rid1 = rid == 0
gen rid2 = rid == 1

eq rid1: rid1
eq rid2: rid2
 
gllamm yvar years logptax loginc rid1 rid2 Xyears Xlogptax Xloginc, ///
  i(obs) nrf(2) eqs(rid1 rid2) lv(rid) fv(rid) family(binomial) ///
  link(probit) nip(7) adapt nocons
* Example 3 end

After many Iterations it converges but is this the model...?  I am still thinking about it and how to recover the
parameters. Unfortunately I do not have a model to compare with - or do you know 
a working example?

Thanks 
Frank







___________________________________________________________
Schon gehört? WEB.DE hat einen genialen Phishing-Filter in die
Toolbar eingebaut! http://produkte.web.de/go/toolbar

*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/