Re: st: Right Censored ordered probit with Stata 11

[Urmi Bhattacharya <ub3@indiana.edu>]

Hi Nick and Austin,

Thank you so much for your suggestions. I put the
set trace on
set traced 1

followed by Owen's code and the simulation. It identified some errors
I could correct by  wrapping up the unwrapped lines . But I am getting
an error now and I can't understand what is causing it. I am providing
the result that I am getting

- version 11.0
  - if replay() {
    if "`e(cmd)'" != "coprbt" {
    error 301
    }	
    syntax [, level(passthru)]
    }
  - else {
  - syntax varlist [if] [in] , censored(varname) [level(passthru) *]
  - marksample touse
  - preserve
  - qui keep if `touse'
  = qui keep if __000000
  - gettoken depvar indepvars : varlist
  - tempvar newdepvar
  - egen `newdepvar'=group(`depvar')
  = egen __000001=group(y)
  - qui sum `newdepvar'
  = qui sum __000001
  - local maxvalue =r(max)
  - forvalues i=1(1)`maxvalue' {
  = forvalues i=1(1)3 {
  - tempvar ucsample sumucs
  - qui gen `ucsample'=0
  = qui gen __000002=0
  - qui replace `ucsample'=(`censored'!=1) if (`newdepvar'==`i')
  = qui replace __000002=(c!=1) if (__000001==1)
  - qui gen `sumucs'=sum(`ucsample')
  = qui gen __000003=sum(__000002)
  - local errorcount =0
  - capture assert `sumucs'[_N]>0&`sumucs'[_N]<.
  = capture assert __000003[_N]>0&__000003[_N]<.
  - if _rc==9 {
    di as error "WARNING: value number `i' has no uncensored obs"
    }
  - drop `sumucs'
  = drop __000003
  - drop `ucsample'
  = drop __000002
  - }
  - tempvar ucsample sumucs
  - qui gen `ucsample'=0
  = qui gen __000004=0
  - qui replace `ucsample'=(`censored'!=1) if (`newdepvar'==`i')
  = qui replace __000004=(c!=1) if (__000001==2)
  - qui gen `sumucs'=sum(`ucsample')
  = qui gen __000005=sum(__000004)
  - local errorcount =0
  - capture assert `sumucs'[_N]>0&`sumucs'[_N]<.
  = capture assert __000005[_N]>0&__000005[_N]<.
  - if _rc==9 {
    di as error "WARNING: value number `i' has no uncensored obs"
    }
  - drop `sumucs'
  = drop __000005
  - drop `ucsample'
  = drop __000004
  - }
  - tempvar ucsample sumucs
  - qui gen `ucsample'=0
  = qui gen __000006=0
  - qui replace `ucsample'=(`censored'!=1) if (`newdepvar'==`i')
  = qui replace __000006=(c!=1) if (__000001==3)
  - qui gen `sumucs'=sum(`ucsample')
  = qui gen __000007=sum(__000006)
  - local errorcount =0
  - capture assert `sumucs'[_N]>0&`sumucs'[_N]<.
  = capture assert __000007[_N]>0&__000007[_N]<.
  - if _rc==9 {
    di as error "WARNING: value number `i' has no uncensored obs"
    }
  - drop `sumucs'
  = drop __000007
  - drop `ucsample'
  = drop __000006
  - }
  - local k =`maxvalue'-1
  = local k =3-1
  - global OGLFXH "`censored'"
  = global OGLFXH "c"
  - qui oprobit `newdepvar' `indepvars'
  = qui oprobit __000001  x1 x2 x3
  - tempname starting
  - matrix `starting'=e(b)
  = matrix __000008=e(b)
  - tempname cutvector
  - matrix `cutvector'=`starting'[1,"_cut1".."_cut`k'"]
  = matrix __000009=__000008[1,"_cut1".."_cut2"]

_cut2 not found

    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'"
    }
    ml model lf coprbt_ll (`newdepvar'=`indepvars' , nocons ) `cutlist'
    if `touse' , maximize title("Homemade Censored Ordered Probit")
init(`initlist') `options'
    ereturn local cmd "coprbt"
    }

r(111);

end of do-file

r(111);

The "_cut2 not found" is in red suggesting that is where the problem
is, but I cannot figure out what is causing it.

Any suggestions?

Thanks a lot.

Urmi
On Sat, Sep 24, 2011 at 8:14 PM, Austin Nichols <austinnichols@gmail.com> wrote:
> Urmi Bhattacharya <ub3@indiana.edu>:
> I also guess that you have not got the line structure right in your
> files, but I will also note that a quick simulation I ran indicated
> that -oprobit- run only on uncensored obs outperforms unrestricted
> -oprobit- and -coprbit-.  So start with
>
> .  oprobit y x1 x2 x3 if censored==0
>
> I am also guessing that other models would outperform all 3 options
> mentioned here, but that will have to wait for another day.
>
> On Sat, Sep 24, 2011 at 4:08 AM, Nick Cox <njcoxstata@gmail.com> wrote:
>> 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/
>

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