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