Dear Edoardo, Shareen, and others,
I've generalized the .do file for the censored ordered
probit into an .ado file. I tried to include as many
features as possible. At the moment, it is byable,
and like cnsreg takes an argument censored(varname)
where the variable is set to -1 for left censored
observations, 0 for uncensored observations, and 1 for
right censored (so you can use an enrollment binary
for the censoring variable in schooling attainment
regressions). It also takes the ml options. In Stata
8 these include "svy" which invokes the stored survey
settings (may be handy if the attainment data is from
a stratified survey). Also, you may want to specify
"difficult" if it doesn't converge.
I'll see what I can do about writing a full help file.
It might take some time, though, as this is the first
time I've written a command.
In response to your earlier concerns, Edoardo, I think
it's probably okay to allow negative cutpoints. I'm
pretty sure that when you see results with the
cutpoints all above zero, they have been rescaled (I put
up a confused and confusing post about this: "st oprobit
cutpoints"). Also, the .ado file takes initial values
from an ordered probit model (without censoring,
survey settings, etc.).
I hope that people find this new file useful, and
would appreciate any programming advice. In
particular, I would like to know if there's a better
way to pass the name of the censoring variable to the
evaluator (I'm currently using a global macro). Also,
I tried to use vectorized commands rather than loops
in the evaluator, but it is fairly slow, so I'd love
any advice on style and efficiency.
cheers,
Owen
*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
--- edoardo masset <[email protected]> wrote:
> dear Shareen,
>
> I've written this do-file for the estimation of a
> censored ordered probit. Many thanks to William Gould
> for his help. The file estimates school attainment as
> in your case. The likelihood function is a modified
> ordered probit function (as in Maddala pag 48) and is
> the sum of the LF for the censored and the uncensored
> observations. I limited the number of outcomes to 7
> because I found some problems in the estimation:
> - you might get negative cut-off points, I think these
> should be constrained to be larger than zero, though
> I've seen published papers with negative cutoffs. But
> I dont know how to do this.
> - I guess it would work better using initial values
> estimated on an ordinary ord probit or on the first
> equation of the model. But again I have to find out
> how to do this.
>
> hope this can help
> Edoardo
>
>
> ***CENSORED ORDERED PROBIT MODEL*****
>
> /*'att', children attending school are censored
> ovbservations.
> The program maxcens maximises a likelihood function
> which is the sum of the LF of uncensored and
> uncensored observations.
> NOTE: convergence is achieved with difficulty
> cut-off values are out of sensible range*/
>
>
> /*define school attainment dependent variable*/
> gen school1 =school==0
> /*never been to school*/
> gen school2 =school>0 & school<4
> /*first half primary*/
> gen school3 =school>3 & school<7
> /*second half primary*/
> gen school4 =school==7
> /*jss1*/
> gen school5 =school==8
> /*JSS2*/
> gen school6 =school==9
> /*JSS3*/
> gen school7 =school>9
> /*all beyond*/
>
> gen unc=att==0 /*define
> uncensored observation: children not attending
> school*/
>
> cap program drop maxcens
> program define maxcens
> args lnf1 theta1 cut1 cut2 cut3 cut4 cut5 cut6
> qui replace `lnf1' =
> unc*(ln(($ML_y1*normprob(`cut1'-`theta1'))+/*
> */
> ($ML_y2*(normprob(`cut2'-`theta1')-normprob(`cut1'
> -`theta1')))+ /*
> */
> ($ML_y3*(normprob(`cut3'-`theta1')-normprob(`cut2'
> -`theta1')))+ /*
> */
> ($ML_y4*(normprob(`cut4'-`theta1')-normprob(`cut3'
> -`theta1')))+ /*
> */
> ($ML_y5*(normprob(`cut5'-`theta1')-normprob(`cut4'
> -`theta1')))+ /*
> */
> ($ML_y6*(normprob(`cut6'-`theta1')-normprob(`cut5'
> -`theta1')))+ /*
> */
> ($ML_y7*(1-normprob(`cut6' -`theta1')))))+ /*
> */ att*(ln(($ML_y1)+/*
> */
> ($ML_y2*(1-normprob(`cut1'-`theta1')))+ /*
> */
> ($ML_y3*(1-normprob(`cut2'-`theta1')))+ /*
> */
> ($ML_y4*(1-normprob(`cut3'-`theta1')))+ /*
> */
> ($ML_y5*(1-normprob(`cut4'-`theta1')))+ /*
> */
> ($ML_y6*(1-normprob(`cut5'-`theta1')))+ /*
> */
> ($ML_y7*(1-normprob(`cut6' -`theta1')))))
>
> end
>
> ml model lf maxcens (school1 school2 school3 school4
> school5 school6 school7= rural sex age ,nocons) /cut1
> /cut2 /cut3 /cut4 /cut5 /cut6
> ml search, repeat(100)
> ml maximize, difficult
>
> log close
>
>
>
>
>
> __________________________________
> Do you Yahoo!?
> Yahoo! SiteBuilder - Free, easy-to-use web site design software
> http://sitebuilder.yahoo.com
> *
> * For searches and help try:
> * http://www.stata.com/support/faqs/res/findit.html
> * http://www.stata.com/support/statalist/faq
> * http://www.ats.ucla.edu/stat/stata/
*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
