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

From |
Charlotte Gary <charlottegary83@yahoo.com> |

To |
statalist@hsphsun2.harvard.edu |

Subject |
Re: st: RE: RE: AGAIN on mvnp code Cappellari & Jenkins 2006 |

Date |
Wed, 8 Jul 2009 08:08:05 -0700 (PDT) |

Dear Nick, really thank you now it works perfectly! May I ask you once again another suggestion. I'm trying to compare the results obtained using ml (paragraph 3.4 cappellari and jenkins ) with those obtained using mvprobit instead of ml. I just would like to be sure that I’m estimating the predicted probabilities of multiple outcome variables after mvprobit correctly. hence my question is: Does anyone see something wrong in the following of proceeding? I followed what Cappellari and Jenkins SJ 2006 suggest at page 166, that is: • Step 1→ estimate the multivariate probit model using mvprobit • Step 2→ save the upper integration point variables – i.e. the 3 (M=3) linear index variables Im =bM’xM which can be derived by using mvppred xb after mvprobit – In this way I obtained xb1, xb2, xb3. • Step 3→ save the V matrix. I used mvppred stdp after mvprobit. I obtained stdp1, stdp2, stdp3. Then I used correlate stdp1 stdp2 stdp3 and I named this matrix V. • Step 4→ create the Cholesky matrix using: matrix c=cholesky(V) • Step 5→ create the signs variable k1=2*y1-1, k2=2*y2-1, k3=2*y3-1 • Step 6→ use mdraws to generate the draws variables – e.g.:mdraws, dr(250) neq(3) prefix(z) random seed(123456789) antithetics replace • Step 7→ calculate the probabilities by using the egen command with the Im variables as arguments and referring to a matrix equal to cholesky(V ) in the chol() option. i.e.: egen pr_joint=mvnp( xb1 xb2 xb3), dr(100) chol(V) prefix(ff) signs(k1 k2 k3) --- On Wed, 7/8/09, Nick Cox <n.j.cox@durham.ac.uk> wrote: > From: Nick Cox <n.j.cox@durham.ac.uk> > Subject: st: RE: RE: AGAIN on mvnp code Cappellari & Jenkins 2006 > To: statalist@hsphsun2.harvard.edu > Date: Wednesday, July 8, 2009, 1:00 PM > You definitely should not change the > code > > egen `sp' = > > to > > egen sp = > > as that produces the error you report. > > From my understanding, what I would suggest is > > 0. Before you call this code, > > gen safecopy = . > > 1. Within your -myll- code add as last line > > replace safecopy = `sp' > > 2. Access the variable -safecopy- afterwards. > > Nick > n.j.cox@durham.ac.uk > > > Charlotte Gary > > My purpose is to get the probabilities of every possible > multiple outcome, therefore sp in the code of Cappellari and > Jenkins 2006 (p.168) is the variable of my interest. > However, in the code (see again below) sp is defined > as a tempvar hence it is not saved as a new variable in the > dataset. I tried to modify the program omitting sp from the > tempvar list and changing this line: > egen `sp'= mvnp(`xb1' `xb2' `xb3'), chol(`C') draws($dr) > prefix(z) signs(`k1' `k2' `k3' ) > > into this: > egen sp= mvnp(`xb1' `xb2' `xb3'), chol(`C') draws($dr) > prefix(z) signs(`k1' `k2' `k3' ) > > However, this was of proceeding is not correct as after ml > maximize an error message appears: “sp is already > defined”. > Hence my question is: how is it possible to obtain the > probabilities of every possible multiple outcome at unit of > analysis level? > Really thank you in advance! > > > program define myll > args lnf xb1 xb2 xb3 c21 c31 c32 > tempvar sp k1 k2 k3 > quietly { > gen double `k1' = 2*$ML_dummy_C - 1 > gen double `k2' = 2*$ML_dummy_I - 1 > gen double `k3' = 2*$ML_dummy_S - 1 > tempname cf21 cf22 cf31 cf32 cf33 C > su `c21', meanonly > scalar `cf21' = r(mean) > su `c31', meanonly > scalar `cf31' = r(mean) > su `c32', meanonly > scalar `cf32' = r(mean) > scalar `cf22' = sqrt( 1 - `cf21'^2 ) > scalar `cf33' = sqrt( 1 - `cf31'^2 - `cf32'^2 ) > mat `C' = (1, 0, 0 \ `cf21', `cf22', 0 \ `cf31',`cf32', > `cf33') > egen `sp'= mvnp(`xb1' `xb2' `xb3'), chol(`C') draws($dr) > prefix(z) signs(`k1' `k2' `k3' ) > replace `lnf'= ln(`sp') > > } > end > > quietly { > probit y1 x1 > mat b1 = e(b) > mat coleq b1 = y1 > probit y2 x1 x2 > mat b2 = e(b) > mat coleq b2 = y2 > probit y3 x1 x2 x3 > mat b3 = e(b) > mat coleq b3 = y3 > mat b0 = b1, b2, b3 > } > > mdraws, dr(250) neq(3) prefix(z) random seed(123456789) > antithetics replace > > > global dr = r(n_draws) > ml model lf myll (y1: y1=x1) (y2: y2=x1 x2) (y3: y3 = x1 x2 > x3)/c21 /c31 /c32, title("MV Probit by MSL, $dr pseudorandom > draws") > ml init b0 > ml maximize > > * > * 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/ > * * 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/

- Prev by Date:
**Re: st: Is anyone aware of the more efficient way to use reshape command** - Next by Date:
**Re: st: data points within radius?** - Previous by thread:
**st: Is anyone aware of the more efficient way to use reshape command** - Next by thread:
**Re: st: RE: RE: AGAIN on mvnp code Cappellari & Jenkins 2006** - Index(es):

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