Re: st: Right Censored ordered probit with Stata 11

[Urmi Bhattacharya <ub3@indiana.edu>]

Hi Austin,

Thanks for your example. But when i run the program and repeat the
same example you provide  I get an error message saying
 "invalid syntax
r(198); "
I am using Stata version 11.2 on a mac.

Is there something I am doing wrong here?

I am providing the whole code i typed, which is essentially running
Owen Haaga's code, followed by the example you provided:


*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



drawnorm x1 x2 x3 e,n(10000) clear
(obs 10000)
g xb=(x1+x2+x3)/2
g y=(xb+e<-2) +2*inrange(xb+e,-2,0)+3*(xb+e>0)

tabulate y
y Freq. Percent Cum.
1 674 6.74 6.74
2 4,380 43.80 50.54
3 4,946 49.46 100.00
Total 10,000 100.00

g c=uniform()>.9
 replace y=1 if c==1
(947 real changes made)
 coprbt y x1 x2 x3, censored(c)
invalid syntax
r(198);

Any help would be greatly appreciated.

Many Thanks

Urmi

On Fri, Sep 23, 2011 at 6:57 PM, Austin Nichols <austinnichols@gmail.com> wrote:
> Urmi Bhattacharya <ub3@indiana.edu> :
> Owen's program runs fine for me, and is far superior to -oprobit- in
> this simple example:
>
> drawnorm x1 x2 x3 e, n(10000) clear
> g xb=(x1+x2+x3)/2
> g y=(xb+e<-2)+2*inrange(xb+e,-2,0)+3*(xb+e>0)
> g c=uniform()>.9
> replace y=1 if c==1
> oprobit y x1 x2 x3
> coprbt y x1 x2 x3, censored(c)
>
>
> On Fri, Sep 23, 2011 at 9:11 AM, Urmi Bhattacharya <ub3@indiana.edu> wrote:
>> Hi Austin,
>>
>> Thanks for the putting the original code.
>>
>> This was the code I referred to ( though I should have put this up
>> myself). It seems to me that the global macros like $ML_Y1, theta1,etc
>> are not defined here which led me to think a part of the code is
>> missing maybe. Or could be that I am going horribly wrong somewhere.
>>
>> Thanks
>>
>> Urmi
>> On Fri, Sep 23, 2011 at 6:17 PM, Austin Nichols <austinnichols@gmail.com> wrote:
>>> Nick--
>>> Looks like the original code is here:
>>> http://www.stata.com/statalist/archive/2003-07/msg00917.html
>>> but I have never tried to estimate a censored ordered probit, and I
>>> strongly suspect convergence will be difficult to be confident of.
>>> It looks like Owen just adapted that original code into an ado for his
>>> own use, and I don't see any further development,
>>> but I am guessing Maarten Buis will have some input on this
>>> problem--wasn't his dissertation research related to this problem?
>>> (see what I did there?  we call that passing the buck in this country)
>>>
>>> On Fri, Sep 23, 2011 at 8:23 AM, Nick Cox <njcoxstata@gmail.com> wrote:
>>>> I don't know what "the code provided by Edoardo" is.
>>>>
>>>> All I can see is the code in the Statalist posting your referred to in
>>>> your first post. I see no obvious problems in the use of global macros
>>>> there.
>>>>
>>>> Someone else may be able to help (much) more, but you'd increase your
>>>> chance of better help by giving precise details of what you mean.
>>>>
>>>> Nick
>>>>
>>>> On Fri, Sep 23, 2011 at 11:03 AM, Urmi Bhattacharya <ub3@indiana.edu> wrote:
>>>>> Thanks Nick.
>>>>> I have been trying to use the code provided by Edoardo, but it seems
>>>>> the global macros are undefined, suggesting a bit of the log file has
>>>>> not been pasted. As for Owen Haaga's code, I can't figure out how to
>>>>> frame the command after running the ado files defining the program.
>>>>> I am hoping someone on the list has run  into a similar problem and
>>>>> would be able to provide me with suggestions on how to work it.
>>>>> Thanks
>>>>>
>>>>> Urmi
>>>>> On Fri, Sep 23, 2011 at 12:47 PM, Nick Cox <njcoxstata@gmail.com> wrote:
>>>>>>
>>>>>> That posting included code, which you can copy. I have no idea how good it is.
>>>>>>
>>>>>> Owen Haaga hasn't been back to Statalist since, but the person with
>>>>>> the same name who works at the Urban Institute may be the same person.
>>>>>> Odds are that Austin Nichols will know.
>>>>>>
>>>>>> Nick
>>>>>>
>>>>>> On Fri, Sep 23, 2011 at 8:02 AM, Urmi Bhattacharya <ub3@indiana.edu> wrote:
>>>>>>
>>>>>> > I am interested in running  a right censored ordered probit for grade
>>>>>> > attainment. I saw in the statalist archive that Owen Haaga had written
>>>>>> > a generalized ado file for doing the same for version 8.
>>>>>> > http://www.stata.com/statalist/archive/2003-08/msg00426.html
>>>>>> >
>>>>>> >
>>>>>> > I was wondering how could I access the ado file  and run it ?
>>>>>> >
>>>>>> > Does someone know how if there is a different program for doing the same?
>>>
>
> *
> *   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/