Bayesian Bivariate Ordered Probit Model Stata: using "llevaluator" feature

[Behram Wali]

I'm trying to estimate a Bayesian Bivariate Ordered Probit Model using Stata's "llevaluator" feature for the first time in Stata/SE 14.1 .
I am using Zurab Sajaia's bioprobit routine in the "llevaluator" in my "bayesmh" command.

The reference to Stata's bioprobit is:
Sajaia, Zurab. "Maximum likelihood estimation of a bivariate ordered probit model: implementation and Monte Carlo simulations." The Stata Journal 4.2 (2008): 1-18.

The code i am trying is:

bayesmh (x50 x22_5 x52_2) (x51 x22_5 x53_2), llevaluator(bioprobit) 
prior ({x50:} {x51:}, normal(0,1000))

Stata returns with giving me the following error:
. bayesmh (x50 x22_5 x52_2) (x51 x22_5 x53_2), llevaluator(bioprobit)
no prior specified for x50:x22_5

I also tried:

bayesmh (x50 x22_5 x52_2) (x51 x22_5 x53_2), llevaluator(bioprobit) 
prior ({x50:} {x51:}, normal(0,1000)) 
prior({x50:_cons} {x51:_cons}, normal(0,1000)) 
prior({x50:x22_5} {x50:x52_5} {x51:x22_5} {x51:x53_2}, normal(0,1000))

Stata returns with the following error:
. bayesmh (x50 x22_5 x52_2) (x51 x22_5 x53_2), llevaluator(bioprobit)
no prior specified for x50:x22_5

I suspect that i am making a huge mistake in my code. Can anyone please guide me on this?

Thank you in advance for your help!
Behram

[Behram Wali]

Following upon my above query, i tried this code too. Can anyone help please?

bayesmh (x50) (x51), llevaluator(bioprobit)
prior ({x50:}, normal(0,1))
prior ({x51:}, normal(0,1))
prior ({x50: cut11 cut12 cut13}, exponential({lambda=30}))
prior ({x51: cut21 cut22 cut23}, exponential({lambda=30}))
prior({rho}, uniform(10,40)) block(lambda) dots

Error:no prior specified for x50:_cons

[Nikolay Balov (StataCorp)]

Behram,

I would expect your bayesmh specification to be something like this

. bayesmh (x50) (x51), llevaluator(bioprobit,                   ///
    param({rho}{x50:cut11 cut12 cut13}{x51:cut21 cut22 cut23})) ///
    prior ({x50:_cons}, normal(0,1))                            ///
    prior ({x51:_cons}, normal(0,1))                            ///
    prior ({x50: cut11 cut12 cut13}, exponential({lambda=30}))  ///
    prior ({x51: cut21 cut22 cut23}, exponential({lambda=30}))  ///
    prior({rho}, uniform(10,40)) block(lambda) dots

I cannot be sure though unless you show me the bioprobit evaluator you are using.

Nikolay

[Behram Wali]

Thanks a lot for your valuable reply. I am using Zurab Sajaia's bioprobit routine in the "llevaluator" in my "bayesmh" command. I am copying the script from bioprobit.ado and is as follow:

*! version 1.11  17Aug2007
#delimit ;

program define bioprobit;
    version 9;
    if replay() {;
        if ("`e(cmd)'" ~= "bioprobit") error 301;
        bioprobit_replay `0';
    }; // end if
    else bioprobit_mx `0';
end; // end program bioprobit

program define bioprobit_replay;
    syntax [,Level(cilevel)];
    local rho diparm(athrho  , tanh label("rho"));
    _coef_table_header;
    display;
    _coef_table, level(`level') `rho' notest;
       display    in green e(chi2_ct) " test of indep. eqns. :"
        _col(38) "chi2(" in yellow "1" in green ") = " in yellow %8.2f e(chi2_c)
        _col(59) in green  "Prob > chi2 = " in yellow %6.4f e(p_c);
        if strpos("`e(title)'", "Simultaneous")    display    in smcl in green "{hline 78}";

end; // end program bioprobit_reply

program define bioprobit_mx, eclass;

    gettoken first : 0, match(paren);
    if missing("`paren'") {;                     // syntax 1, bivariate ordered probit models
           syntax varlist [if] [in] [pweight iweight fweight]
                                        [, offset1(varname) offset2(varname) COLlinear
                                   Robust CLuster(varname) Level(cilevel) end *];
        gettoken y1 varlist : varlist;
        gettoken y2 varlist : varlist;
        local x1 `varlist';
        local x2 `varlist';
        local pref "B";
    };
    else {;                                      // syntax 2, seemingly unrelated bivariate ordered probit model
           syntax anything(id="equation id" equalok) [if] [in] [pweight iweight fweight]
                                       [, offset1(varname) offset2(varname) COLlinear
                                   Robust CLuster(varname) Level(cilevel) end *];
        gettoken eq1 eq2_ :  anything,  match(parns) bind;
        tokenize "`eq1'", parse("=");
        if ("`2'"!="=") {;
            tokenize "`eq1'";
            local y1 `1';
            macro shift;
            local x1 `*';
        };
        else {;
            local y1 `1';
            local x1 `3';
        };
        gettoken eq2 : eq2_ , match(parns) bind;
        if ("`parns'"!="(") local eq2 "`eq2_'";
        tokenize "`eq2'", parse("=");
        if ("`2'"!="=") {;
            tokenize "`eq2'";
            local y2 `1';
            macro shift;
            local x2 `*';
        };
        else {;
            local y2 `1';
            local x2 `3';
        };
        local pref "Seemingly unrelated b";
    }; // if
    global END =~missing("`end'");

    quietly {;

    // define standard dep vars
    tempvar sy1 sy2;
    egen `sy1' = group(`y1');
    egen `sy2' = group(`y2');

    marksample touse;
    markout `touse' `sy1' `sy2' `x1' `x2' `cluster' `offset1' `offset2', strok;

    mlopts mlopts, `options';

    if (~missing("`offset1'")) local offo1 "offset(`offset1')";
    if (~missing("`offset2'")) local offo2 "offset(`offset2')";
       if (~missing("`weight'"))  local weight "[`weight'`exp']";
    if (~missing("`cluster'")) local clopt cluster(`cluster');
    if "`weight'" == "pweight" | (~missing("`cluster'")) local robust "robust";

    // Remove collinear variables
    noisily _rmdcoll `sy1' `x1' `weight' if `touse', `collinear';
    local x1 "`r(varlist)'";
noi tabulate `sy2' if `touse';

    noisily _rmdcoll `sy2' `x2' `weight' if `touse', `collinear';
    local x2 "`r(varlist)'";

    summarize `sy1'; global NC1  = r(max);
     summarize `sy2'; global NC2  = r(max);    // number of categiries
    if ($NC1==1) {;
        noisily display as error "`y1' does not vary";
        exit 2000;
    };
    if ($NC2==1) {;
        noisily display as error "`y2' does not vary";
        exit 2000;
    };
    local NC1_1   = $NC1-1;
    local NC2_1   = $NC2-1;

    // define initial values
    tempname Ib_op1 Ib_op2 Ic_op Ic_op1 Ic_op2 I_rho TMP ll_0;
       scalar `ll_0' = 0;

       oprobit `sy1' `x1' `weight' if `touse', `offo1';
        matrix             `TMP' = e(b);
        matrix       `Ib_op1' = `TMP'[1,"`sy1':"];
        matrix coleq `Ib_op1' = `y1';
        matrix          `Ic_op1' = `TMP'[1,"cut1:_cons".."cut`NC1_1':_cons"];
    scalar `ll_0'=`ll_0'+e(ll);
    if $END {;
        tempvar xb1;
        matrix score `xb1' = `Ib_op1';
    };
    oprobit `sy2' `xb1' `x2' `weight' if `touse', `offo2';
        matrix             `TMP' = e(b);
        matrix          `Ib_op2' = `TMP'[1,"`sy2':"];
        matrix coleq `Ib_op2' = `y2';
        matrix          `Ic_op2' = `TMP'[1,"cut1:_cons".."cut`NC2_1':_cons"];
    scalar `ll_0'=`ll_0'+e(ll);


    matrix `Ic_op'  = `Ic_op1'[1,1];
    forvalues i = 2/`NC1_1' {; matrix `Ic_op' =`Ic_op',sqrt(`Ic_op1'[1,`i']-`Ic_op1'[1,`i'-1]); };
    matrix `Ic_op'  = `Ic_op', `Ic_op2'[1,1];
    forvalues i = 2/`NC2_1' {; matrix `Ic_op' =`Ic_op',sqrt(`Ic_op2'[1,`i']-`Ic_op2'[1,`i'-1]);    };

    correlate `sy1' `sy2' if `touse';
    matrix          `I_rho' = .5;//atanh(r(rho));
    matrix colnames `I_rho' = athrho:_cons;

    local cuts;
    forvalues k = 1/`NC1_1' {; local cuts "`cuts' /cut1`k'"; local ceqs "`ceqs' cut1`k'"; };
    forvalues k = 1/`NC2_1' {; local cuts "`cuts' /cut2`k'"; local ceqs "`ceqs' cut2`k'"; };
    matrix coleq    `Ic_op' = `ceqs';
    matrix colnames `Ic_op' = _cons;
    }; // end quietly

     // maximization
    if $END {;
        matrix `I_rho'  =`I_rho',`Ib_op2'[1,1];
        matrix colnames `I_rho' = athrho:_cons gamma:_cons;

        matrix `Ib_op2' =`Ib_op2'[1,2...];
        local gamma "/gamma";
        local title "title(Simultaneous bivariate ordered probit regression)";
    };
    else local title "title(`pref'ivariate ordered probit regression)";

    ml model d2 bioprobit_d2
         (`y1' : `sy1' = `x1', noconstant)
         (`y2' : `sy2' = `x2', noconstant)
         /athrho
         `gamma'
         `cuts'
         `weight'
         if `touse'
           ,
         `title'
         difficult
         collinear
         missing
         search(on)
         init(`Ib_op1' `Ib_op2' `I_rho' `Ic_op')
         maximize
         `clopt'
         `robust'
           `mlopts'
         ;

    tempname b D D1 V;
    matrix `b'  =e(b);
    matrix `V'  =e(V);
    local f = colnumb(`b',"cut11:_cons") -1;
    matrix `D' = I(`f'), J(`f', `NC1_1'+`NC2_1',0);                                   // matrix of derivatives
    matrix `D1' = (J(1,`f',0), 1, J(1,`NC1_1'+`NC2_1'-1, 0));
    matrix `D'  = `D' \ `D1';
    forvalues k = 2/`NC1_1' {;
        matrix `D1'[1,`f'+`k'] = 2*`b'[1,`f'+`k'];
        matrix `D'  = `D' \ `D1';
           matrix `b'[1,`f'+`k'] =`b'[1,`f'+`k'-1]+`b'[1,`f'+`k']^2;
    };
    local f = colnumb(`b',"cut21:_cons") -1;
    if (`NC2_1' > 1) matrix `D1' = (J(1,`f', 0), 1, J(1,`NC2_1'-1, 0));
    else             matrix `D1' = (J(1,`f', 0), 1);
       matrix `D'  = `D' \ `D1';
       forvalues k = 2/`NC2_1' {;
        matrix `D1'[1,`f'+`k']  = 2*`b'[1,`f'+`k'];
        matrix `D'  = `D' \ `D1';
        matrix `b'[1,`f'+`k'] =`b'[1,`f'+`k'-1]+`b'[1,`f'+`k']^2;
    };
    matrix `V' = `D'*`V'*`D'';
    ereturn repost b=`b' V=`V';

    ereturn scalar ll_0 = `ll_0';               // loglikelihood for non-correlated case
    ereturn scalar k_aux= `NC1_1'+`NC2_1';      // identify ancillary parameters

    ereturn local cmd    "bioprobit";
    ereturn local predict "bioprobit_p";
    ereturn local depvar "`y1' `y2'";
    ereturn local offset1 `offset1';
    ereturn local offset2 `offset2';

    if missing("`robust'") {;
        ereturn local chi2_ct "LR";
        ereturn scalar chi2_c = 2 * (e(ll) - `ll_0');
    };
    else {;
        ereturn local chi2_ct "Wald";
        quietly test [athrho]_b[_cons];
        ereturn scalar chi2_c = r(chi2);
        if $END {;
            ereturn local
            quietly test [athrho]_b[_cons]+[gamma]_b[_cons]==0;
            ereturn scalar chi2_c2 = r(chi2);
            ereturn scalar p_c2    = chiprob(2, e(chi2_c2));
        };
    };
    ereturn scalar p_c = chiprob(1, e(chi2_c));
    global END;
    bioprobit_replay, level(`level');
end; // end program bioprobit_mx

I also tried the specification you suggested but Stata gives the following error:

. bayesmh (x50) (x51), llevaluator(bioprobit,                   ///
>     param({rho}{x50:cut11 cut12 cut13}{x51:cut21 cut22 cut23})) ///
>     prior ({x50:_cons}, normal(0,1))                            ///
>     prior ({x51:_cons}, normal(0,1))                            ///
>     prior ({x50: cut11 cut12 cut13}, exponential({lambda=30}))  ///
>     prior ({x51: cut21 cut22 cut23}, exponential({lambda=30}))  ///
>     prior({rho}, uniform(10,40)) block(lambda) dots
     _check_duplicates():  3301  subscript invalid
                 <istmt>:     -  function returned error
r(3301);

end of do-file

r(3301);

Please note that there are two other ado files too (under) bioprobit, named as bioprobit_d2.ado and bioprobit_p.ado. I am not sure which ado file should i use for my case.

I thank you for your time and attention.
-Behram

[Nikolay Balov (StataCorp)]

Behram,

I cannot reproduce the error because I don't have the data, but regardless, the bioprobit program you are trying to use is a maximum likelihood evaluator and it is not compatible with bayesmh. Log-likelihood evaluators in bayesmh have specific syntax and structure, see http://www.stata.com/manuals14/bayesbayesmhevaluators.pdf.

If you understand the bioprobit model you may try to extract the relevant code from the Zurab Sajaia's bioprobit routine which computes the log-likelihood for a given set of parameters and write a bayesmh log-likelihood evaluator. Unfortunately, I am not familiar with the model and presently don't have such an evaluator ready.

Nikolay

[Behram Wali]

Thank you Nikolay for your input. I have just started studying how to program user specified models, by studying and practicing examples in Maximum Likelihood Estimation With Stata by Gould and Pitblado. Unfortunately, bioprobit's evaluator is too advanced for me to program at this stage.

I would appreciate if someone can guide me on log-likelihood evaluator for bioprobit in bayesmh. Thanks much

[Martin Andresen]

I have a related problem that appears when I try to use prior(...,density(...)) for two or more parameters. Here is the simplest possible example I could cook up, based on the examples from bayesmh.

sysuse auto, clear
bayesmh mpg weight trunk, likelihood(normal({var})) prior({mpg:weight}, density(normalden({mpg:weight})))  prior({mpg:trunk}, density(normalden({mpg:trunk}))) prior({var}, jeffreys) initial({mpg:weight} -0.6 {mpg:trunk} 2 {var} 11.83) nocons

I am aware of the difference betwen linear combinations (option xb in paramspec) and a single coefficient, but just can't get this to work any way I try. As far as I understand {mpg:trunk} above should refer to the parameter, not the linear combination.... Does anyone have any advice?

[Martin Andresen]

I would appreciate if someone can guide me on log-likelihood evaluator for bioprobit in bayesmh. Thanks much

Wouldn't this mostly involve summarizing the individual log likelihood contributions and adding them to the scalar `lnden', which you need to define as the first argument to the evaluator?