Bookmark and Share

Notice: On March 31, it was announced that Statalist is moving from an email list to a forum. The old list will shut down on April 23, and its replacement, statalist.org is already up and running.


[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: st: Right Censored ordered probit with Stata 11


From   Austin Nichols <austinnichols@gmail.com>
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: Right Censored ordered probit with Stata 11
Date   Sat, 24 Sep 2011 10:44:44 -0400

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/


© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   Site index