Margins after user-written censored ordinal probit

[Celia Vera]

Hi everyone,

I want to estimate marginal effects after a censored ordinal probit. Following some links (http://www.stata.com/statalist/archi.../msg00090.html and http://www.stata.com/statalist/archi.../msg00917.html), I have written a maximum likelihood procedure to get the censored ordinal probit and it works fine. However, I am still unable to get the code for the marginal effects. I would appreciate your advise on how to proceed. Here what I have so far:

gen educ21= educ2==1
gen educ22= educ2==2
gen educ23= educ2==3
                                                     
//define uncensored observation: children not attending school stud==0
gen unc=stud==0   
     
cap program drop lfoprobit5 
program define lfoprobit5
args lnf1 theta1 cut1 cut2 
qui replace `lnf1' = unc*(ln(($ML_y1*normprob(`cut1'-`theta1'))+ ($ML_y2*(normprob(`cut2'-`theta1')-normprob(`cut1'-`theta1')))+ ($ML_y3*(1-normprob(`cut2' -`theta1')))))+ ///
stud*(ln(($ML_y1)+ ($ML_y2*(1-normprob(`cut1'-`theta1')))+ ($ML_y3*(1-normprob(`cut2'-`theta1')))))
ereturn local predict ‘"‘lfoprobit5_p’"’ 
end

cap program drop lfoprobit5_p
program define lfoprobit5_p
                syntax newvarlist [if] [in], [x]
                        tempvar touse xb
                        marksample touse // mark `touse'
                        qui _predict double `xb' 
gen double `varlist'=norm([_cut1]_b[_cons]-`xb') if `touse'
end

ml model lf lfoprobit5 (educ21 educ22 educ23 = status sex edad, nocons)  /cut1 /cut2, svy subpop(if edad>=25)  
ml search, repeat(100) 
ml maximize, difficult

I get :

                            
            Linearized
        Coef.    Std. Err.    t    P>t    [95% Conf.    Interval]
                            
eq1        
    status    -.0615345    .0421942    -1.46    0.145    -.1442564    .0211875
    sex       -.1580939    .0288409    -5.48    0.000    -.2146367    -.1015512
    edad      -.0156747    .001728    -9.07    0.000    -.0190624    -.012287
                            
cut1        
    _cons    -1.687997    .0658898    -25.62    0.000    -1.817174    -1.558819
                            
cut2        
    _cons    -.7937602    .0643738    -12.33    0.000    -.9199653    -.6675551

Then, I use margins to get, for instance, marginal effects on the first outcome (out of three) :

margins, dydx(*) predict(outcome(#1))

I get the same coefficients than the ordered probit:

Average marginal effects                          Number    of obs   =    46705
Model VCE    : Linearized

Expression   : Linear prediction, predict(outcome(#1))
dy/dx w.r.t. : status sex edad

        
Delta-method
dy/dx   Std. Err.      z    P>z    [95% Conf.    Interval]
        
status   -.0615345   .0421942    -1.46   0.145    -.1442336    .0211646
sex      -.1580939   .0288409    -5.48   0.000    -.2146211    -.1015668
edad     -.0156747    .001728    -9.07   0.000    -.0190615    -.0122879

[Richard Williams]

You will need to write your own predict command or else figure out how to use an existing one. In your program you would have something like

ereturn local predict "cprob_p"

The predict routine for oprobit is oprobi_p.ado. Maybe you could modify it.

The other possibility might be to use the expression option on margins. I think you'd have to write out the formula that computes the predicted probability. i.e. convert the xb's into the corresponding predicted probabilities.

[Celia Vera]

Hi Richard,

Thanks a lot for your remark.
I am doing this because my final aim is to estimate an ordinal probit with censored data and the oprobit command does not allow this. Thus, I was first trying to get marginal effects in a regular probit and test it with margins after oprobit. After verifying I get the same, I would only modify the ordinal probit and incorporate censoring observations (I already know how to do this part).

Following your inputs in this thread and in a previous one (http://www.statalist.org/forums/foru...ordinal-probit), I finally made it !

However, my code only runs with mfx. The problem with margins arises when I declare factor variables as such (i.male instead of male in the oprobit equation). Since mfx identifies the nature of the variables from the possible values, there is no problem when I use mfx. I would like to use the advantages that margins have over mix and thus wonder what I am not doing right. I have replicated the problem with a public data. Any input on where should I look at ?

cap program drop lfoprobit1
program lfoprobit1
  args lnf xb a1 a2
  local y "$ML_y1"
  quietly replace `lnf' =  ln(  normal(`a1'-`xb')) if `y'==1
  quietly replace `lnf' =  ln(  normal(`a2'-`xb') - normal(`a1'-`xb'))  if `y'==2
  quietly replace `lnf' =  ln(1-normal(`a2'-`xb')) if `y'==3
end

cap program drop lfoprobit6
program lfoprobit6, eclass
ml model lf lfoprobit1 (rep78 = i.foreign length mpg, nocons)  /cut1 /cut2
ml search, repeat(100) 
ml maximize, difficult
eret local predict="lfoprobit6_p"
eret local cmd="lfoprobit6"
end

cap program drop lfoprobit6_p
program define lfoprobit6_p
version 9, missing
syntax [anything] [if] [in] [, * ]
if index(`"`anything'"',"*") {
        ParseNewVars `0'
        local varspec `s(varspec)'
        local varlist `s(varlist)'
        local typlist `s(typlist)'
        local if `"`s(if)'"'
        local in `"`s(in)'"'
        local options `"`s(options)'"'
}
else {
        local varspec `anything'
        syntax [newvarlist] [if] [in] [, * ]
}
local nvars : word count `varlist'

ParseOptions, `options'
local type `s(type)'
local outcome `s(outcome)'
if "`type'" != "" {
        local `type' `type'
}
else {
        if `"`outcome'"' != "" {
                di in gr "(option pr assumed; predicted probability)"
        }
        else {
                di in gr "(option pr assumed; predicted probabilities)"
        }
}
version 6, missing
tempvar xb2 xb3 xb4 touse
mark `touse' `if' `in'
qui _predict double `xb2' if `touse', xb `offset'
qui _predict double `xb3' if `touse', xb `offset'
qui _predict double `xb4' if `touse', xb `offset'
if `"`outcome'"' == "1"{
gen double `varlist'=normprob(_b[/cut1]-`xb2') if `touse'
}
else if `"`outcome'"' == "2"{
gen double `varlist'=normprob(`xb3' - _b[/cut1]) - normprob(`xb3' - _b[/cut2]) if `touse'
}
else if `"`outcome'"' == "3"{
gen double `varlist'=normprob(`xb4' - _b[/cut2]) if `touse'
}
label var `varlist'    
end                        

cap program drop ParseNewVars
program ParseNewVars, sclass
        version 9, missing
        syntax [anything(name=vlist)] [if] [in] [, * ]

        if missing(e(version)) {
                local old oldologit
        }
        _score_spec `vlist', `old'
        sreturn local varspec `vlist'
        sreturn local if        `"`if'"'
        sreturn local in        `"`in'"'
        sreturn local options   `"`options'"'
end

cap program drop ParseOptions
program ParseOptions, sclass
        version 9, missing
        syntax [,                       ///
                Outcome(string)         ///
                EQuation(string)        ///
                Index                   ///
                XB                      ///
                STDP                    ///
                Pr                      ///
                noOFFset                ///
                SCores                  ///
                SCore(string)           ///
        ]

        // check options that take arguments
        if `"`equation'"' != "" & `"`score'"' != "" {
                di as err ///
                "options score() and equation() may not be combined"
                exit 198
        }
        if `"`score'"' != "" & `"`outcome'"' != "" {
                di as err ///
                "options score() and outcome() may not be combined"
                exit 198
        }
        if `"`equation'"' != "" & `"`outcome'"' != "" {
                di as err ///
                "options equation() and outcome() may not be combined"
                exit 198
        }
        local eq `"`score'`equation'`outcome'"'

        // check switch options
        local type `index' `xb' `stdp' `pr' `scores'
        if `:word count `type'' > 1 {
                local type : list retok type
                di as err "the following options may not be combined: `type'"
                exit 198
        }
        if !inlist("`type'","","scores") & `"`score'"' != "" {
                di as err "options `type' and score() may not be combined"
                exit 198
        }
        if `"`score'"' != "" {
                local scores scores
        }

        // save results
        sreturn clear
        sreturn local type      `type'
        sreturn local outcome   `"`eq'"'
end

clear
sysuse auto
replace rep=1 if rep==2 | rep==5| rep==.
recode rep (3=2)
recode rep (4=3)

lfoprobit6
margins, dydx(*) predict(outcome(1)) atmeans
I get
Expression   : predict(outcome(1))
dy/dx w.r.t. : 1.foreign length mpg
at           : 0.foreign       =    .7027027 (mean)
1.foreign       =    .2972973 (mean)
length          =    187.9324 (mean)
mpg             =     21.2973 (mean)

        
Delta-method
dy/dx   Std. Err.      z    P>z    [95% Conf. Interval]
        
1.foreign   (not estimable)
length     (not estimable)
mpg        (not estimable)
        
Note: dy/dx for factor levels is the discrete change    from    the base level.

Best.
Celia P.

[Alok Bohara]

Hi Celia:

I am curious to see your code of long-hand oprobit with a margins / expression command. Mine returns ... Alok

>
. ml model lf oprobit_prog ///
(xc:num_cigs = $m1, noconstant if current_ind == 1) /cut1 /cut2 /cut3 [pweight ///
= FinalWgt], vce(cluster sch_class) technique(nr bfgs)

margins, dydx(sp_pca_1) expression(predict(eq(xc)))

gives me:

Average marginal effects Number of obs = 249
Model VCE : Robust

Expression : Linear prediction, predict(eq(xc))
dy/dx w.r.t. : sp_pca_1

------------------------------------------------------------------------------
| Delta-method
| dy/dx Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
sp_pca_1 | . (not estimable)
------------------------------------------------------------------------------