[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

From |
Owen Haaga <ohaaga@yahoo.com> |

To |
statalist@hsphsun2.harvard.edu |

Subject |
Re: st: censored ordered probit model (includes .ado) |

Date |
Fri, 22 Aug 2003 10:10:21 -0700 (PDT) |

Dear Edoardo, Shareen, and others, I've generalized the .do file for the censored ordered probit into an .ado file. I tried to include as many features as possible. At the moment, it is byable, and like cnsreg takes an argument censored(varname) where the variable is set to -1 for left censored observations, 0 for uncensored observations, and 1 for right censored (so you can use an enrollment binary for the censoring variable in schooling attainment regressions). It also takes the ml options. In Stata 8 these include "svy" which invokes the stored survey settings (may be handy if the attainment data is from a stratified survey). Also, you may want to specify "difficult" if it doesn't converge. I'll see what I can do about writing a full help file. It might take some time, though, as this is the first time I've written a command. In response to your earlier concerns, Edoardo, I think it's probably okay to allow negative cutpoints. I'm pretty sure that when you see results with the cutpoints all above zero, they have been rescaled (I put up a confused and confusing post about this: "st oprobit cutpoints"). Also, the .ado file takes initial values from an ordered probit model (without censoring, survey settings, etc.). I hope that people find this new file useful, and would appreciate any programming advice. In particular, I would like to know if there's a better way to pass the name of the censoring variable to the evaluator (I'm currently using a global macro). Also, I tried to use vectorized commands rather than loops in the evaluator, but it is fairly slow, so I'd love any advice on style and efficiency. cheers, Owen *Start cprbt.ado --- capture program drop coprbt program define coprbt , eclass byable(recall) version 8.0 *Replay support. if replay() { if "`e(cmd)'" != "coprbt" { error 301 } syntax [, level(passthru)] } *Main program else { syntax varlist [if] [in] , censored(varname) [level(passthru) *] marksample touse preserve qui keep if `touse' gettoken depvar indepvars : varlist *Rescale the independent variable into categories one through maxvalue tempvar newdepvar egen `newdepvar'=group(`depvar') qui sum `newdepvar' local maxvalue =r(max) *Check to make sure there are uncensored obs for each value forvalues i=1(1)`maxvalue' { tempvar ucsample sumucs qui gen `ucsample'=0 qui replace `ucsample'=(`censored'!=1) if (`newdepvar'==`i') qui gen `sumucs'=sum(`ucsample') local errorcount =0 capture assert `sumucs'[_N]>0&`sumucs'[_N]<. if _rc==9 { di as error "WARNING: value number `i' has no uncensored obs" } drop `sumucs' drop `ucsample' } *Store in local k the number of cutpoints that will be necessary local k =`maxvalue'-1 *Pass the name of the censoring variable to the evaluator (as a global macro) global OGLFXH "`censored'" /*Passes the name of the censoring variable*/ *Estimate a regular ordered probit to get initial value estimates qui oprobit `newdepvar' `indepvars' *Take the coefficients and cutpoints from the ordered probit and store them. tempname starting matrix `starting'=e(b) tempname cutvector matrix `cutvector'=`starting'[1,"_cut1".."_cut`k'"] *Put together local macros with the starting values forvalues i=1(1)`k' { local cut`i'=`cutvector'[1,`i'] local initlist "`initlist' /_cut`i'=`cut`i''" local cutlist "`cutlist' /_cut`i'" } foreach var of local indepvars { local `var'_b =`starting'[1,colnumb(`starting',"`var'")] } foreach var of local indepvars { local initlist "`initlist' `var'=``var'_b'" } *And finally, estimate the model. ml model lf coprbt_ll (`newdepvar'=`indepvars' , nocons ) `cutlist' /// if `touse' , maximize title("Homemade Censored Ordered Probit") /// init(`initlist') `options' ereturn local cmd "coprbt" } ml display , `level' end *END coprbt.ado *START coprbt_ll.ado *! v. 1.0 for use with coprbt.ado -Owen Haaga 22 August 2003 capture program drop coprbt_ll program define coprbt_ll local i=0 while "``i''"!="" { local ++i } local --i /*i should contain the number of arguments*/ local maxvalue=`i'-1 /*2 other arguments plus k cutpoints is i, i-1=k+1=maxvalue*/ local k=`maxvalue'-1 /*k is the total number of cutpoints*/ forvalues i=1(1)`k' { local cuts "`cuts' cut`i'" } args lnf theta `cuts' tempname cutvector matrix `cutvector' =J(1,`k',0) forvalues i=1(1)`k' { local cutvalue=`cut`i''[1] /*WEAK. What if the first obs is not in the sample?*/ matrix `cutvector'[1,`i']=`cutvalue' } local censored "$OGLFXH" /*Any other way to pass this?*/ di "$OGLFXH is in OGLFXH" tempvar y1 ymax yother qui gen `y1'=($ML_y1==1) qui gen `ymax'=($ML_y1==`maxvalue') qui gen `yother'=($ML_y1!=1)&($ML_y1!=`maxvalue') qui replace `lnf'=ln(norm(`cutvector'[1,$ML_y1]-`theta')) if `y1'==1&`censored'==0 qui replace `lnf'=ln(norm(`cutvector'[1,$ML_y1]-`theta')-norm(`cutvector'[1,$ML_y1-1]-`theta')) if `yother'==1&`censored'==0 qui replace `lnf'=ln(1-norm(`cutvector'[1,$ML_y1-1]-`theta')) if `ymax'==1&`censored'==0 qui replace `lnf'=ln(1) if `y1'==1&`censored'==1 qui replace `lnf'=ln(1-norm(`cutvector'[1,$ML_y1-1]-`theta')) if `yother'==1&`censored'==1 qui replace `lnf'=ln(1-norm(`cutvector'[1,$ML_y1-1]-`theta')) if `ymax'==1&`censored'==1 qui replace `lnf'=ln(norm(`cutvector'[1,$ML_y1]-`theta')) if `y1'==1&`censored'==-1 qui replace `lnf'=ln(norm(`cutvector'[1,$ML_y1]-`theta')) if `yother'==1&`censored'==-1 qui replace `lnf'=ln(1) if `ymax'==1&`censored'==-1 end *END coprbt_ll.ado --- edoardo masset <edomasset@yahoo.com> wrote: > dear Shareen, > > I've written this do-file for the estimation of a > censored ordered probit. Many thanks to William Gould > for his help. The file estimates school attainment as > in your case. The likelihood function is a modified > ordered probit function (as in Maddala pag 48) and is > the sum of the LF for the censored and the uncensored > observations. I limited the number of outcomes to 7 > because I found some problems in the estimation: > - you might get negative cut-off points, I think these > should be constrained to be larger than zero, though > I've seen published papers with negative cutoffs. But > I dont know how to do this. > - I guess it would work better using initial values > estimated on an ordinary ord probit or on the first > equation of the model. But again I have to find out > how to do this. > > hope this can help > Edoardo > > > ***CENSORED ORDERED PROBIT MODEL***** > > /*'att', children attending school are censored > ovbservations. > The program maxcens maximises a likelihood function > which is the sum of the LF of uncensored and > uncensored observations. > NOTE: convergence is achieved with difficulty > cut-off values are out of sensible range*/ > > > /*define school attainment dependent variable*/ > gen school1 =school==0 > /*never been to school*/ > gen school2 =school>0 & school<4 > /*first half primary*/ > gen school3 =school>3 & school<7 > /*second half primary*/ > gen school4 =school==7 > /*jss1*/ > gen school5 =school==8 > /*JSS2*/ > gen school6 =school==9 > /*JSS3*/ > gen school7 =school>9 > /*all beyond*/ > > gen unc=att==0 /*define > uncensored observation: children not attending > school*/ > > cap program drop maxcens > program define maxcens > args lnf1 theta1 cut1 cut2 cut3 cut4 cut5 cut6 > qui replace `lnf1' = > unc*(ln(($ML_y1*normprob(`cut1'-`theta1'))+/* > */ > ($ML_y2*(normprob(`cut2'-`theta1')-normprob(`cut1' > -`theta1')))+ /* > */ > ($ML_y3*(normprob(`cut3'-`theta1')-normprob(`cut2' > -`theta1')))+ /* > */ > ($ML_y4*(normprob(`cut4'-`theta1')-normprob(`cut3' > -`theta1')))+ /* > */ > ($ML_y5*(normprob(`cut5'-`theta1')-normprob(`cut4' > -`theta1')))+ /* > */ > ($ML_y6*(normprob(`cut6'-`theta1')-normprob(`cut5' > -`theta1')))+ /* > */ > ($ML_y7*(1-normprob(`cut6' -`theta1')))))+ /* > */ att*(ln(($ML_y1)+/* > */ > ($ML_y2*(1-normprob(`cut1'-`theta1')))+ /* > */ > ($ML_y3*(1-normprob(`cut2'-`theta1')))+ /* > */ > ($ML_y4*(1-normprob(`cut3'-`theta1')))+ /* > */ > ($ML_y5*(1-normprob(`cut4'-`theta1')))+ /* > */ > ($ML_y6*(1-normprob(`cut5'-`theta1')))+ /* > */ > ($ML_y7*(1-normprob(`cut6' -`theta1'))))) > > end > > ml model lf maxcens (school1 school2 school3 school4 > school5 school6 school7= rural sex age ,nocons) /cut1 > /cut2 /cut3 /cut4 /cut5 /cut6 > ml search, repeat(100) > ml maximize, difficult > > log close > > > > > > __________________________________ > Do you Yahoo!? > Yahoo! SiteBuilder - Free, easy-to-use web site design software > http://sitebuilder.yahoo.com > * > * 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/ * * 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/

- Prev by Date:
**st: RE: RE: Towards publication quality output** - Next by Date:
**Re: st: SP/RM[1,2] EMS** - Previous by thread:
**st: RE: RE: Towards publication quality output** - Next by thread:
**st: Using egen or similar** - Index(es):

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