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

Re: st: censored ordered probit model


From   Owen Haaga <[email protected]>
To   [email protected]
Subject   Re: st: censored ordered probit model
Date   Fri, 22 Aug 2003 12:20:18 -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


*
*   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