Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

# st: gllamm bivariate probit randome effects

 From "Frank Muehsam" <[email protected]> To [email protected] Subject st: gllamm bivariate probit randome effects Date Fri, 18 Feb 2011 10:46:20 +0100 (CET)

```

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/
```

© Copyright 1996–2018 StataCorp LLC   |   Terms of use   |   Privacy   |   Contact us   |   Site index