Re: st: Right Censored ordered probit with Stata 11

[Nick Cox <njcoxstata@gmail.com>]

It works for me as long as you unwrap wrapped lines.

. set trace on
. set traced 1

before running and the program will show you each error.

Nick

On Sat, Sep 24, 2011 at 5:28 AM, Urmi Bhattacharya <ub3@indiana.edu> wrote:
> 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)
>>>>

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