Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

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


From   Owen Haaga <[email protected]>
To   [email protected]
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 <[email protected]> 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/



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