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
Urmi Bhattacharya <[email protected]>
To
[email protected]
Subject
Re: st: Right Censored ordered probit with Stata 11
Date
Mon, 26 Sep 2011 08:20:32 +0530
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 <[email protected]> wrote:
> 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/
>
*
* 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/