Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: st: Right Censored ordered probit with Stata 11
From
Austin Nichols <[email protected]>
To
[email protected]
Subject
Re: st: Right Censored ordered probit with Stata 11
Date
Sat, 24 Sep 2011 10:44:44 -0400
Urmi Bhattacharya <[email protected]>:
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 <[email protected]> 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 <[email protected]> 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 <[email protected]> wrote:
>>> Urmi Bhattacharya <[email protected]> :
>>> 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 <[email protected]> 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 <[email protected]> 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/