Maximum Likelihood Estimation Using ml command.

[Siddiq Niaz]

Dear all,

I am trying to estimate a skewed-logistic (or Type 1 logistic) binary choice model. This is essentially similar to the one Stata's scobit command does. However, the reason i want to code this model using Stata's ml command will soon be evident.

Based on Stata's manual rscobit.pdf , the log-likelihood function for skewed logit can be constructed simply by using a "skewed-logit" density for error terms instead of using "normal" or "logit" density in case of simple probit or logit model respectively.

The skewed logit CDF is (Wikipedia source):

Stata fits this model as (from rscobit.pdf):

Having said all this, I am using the below code, but i get an error of "unvalid syntax" each time i run ml check on the program.

capture program drop abc
program define abc
version 14
   args lnf xb alpha
    local f  "(1 - 1/(1+exp(`xb'))^(-`alpha'))"
    quietly replace `lnf' = ln(   `f'(`xb')) if $ML_y1==1    /// f is the CDF of skewed #logit i assume???
    quietly replace `lnf' = ln(1 - `f'(`xb')) if $ML_y1==0
  end
  ml check
ml model lf abc (foreign = mpg rep78 turn)
ml maximize, difficult tolerance(1e-6) ltolerance(1e-7)

When i run ml check, i get this:

model not defined

And if i omit running ml check, and run ml model and ml maximize command, i get:

unknown function ()

Can anyone please help me correcting any silly mistake I may be doing?

Ultimately, I am want to estimate a ordered skewed logit model where simple ordered logit model can be derived by forcing alpha == 1, and ordered skewed logit can be estimated by calculating the alpha parameter in CDF. Of course, a binary choice model is the first step to doing this, so am trying to replicate results of scobit.


Thank you!

[Maarten Buis]

Here are a couple of fixes. It still does not find feasible values, which suggest there are more problems left in your code.

clear all
sysuse auto
program define abc
version 14
    args lnf xb alpha
    local f  "(1 - 1/(1+exp(`xb'))^(-`alpha'))"
    quietly replace `lnf' = ln(   `f') if $ML_y1==1    // f is the CDF of skewed #logit i assume???
    quietly replace `lnf' = ln(1 - `f') if $ML_y1==0
end
ml model lf abc (foreign = mpg rep78 turn) /alpha
ml check
ml maximize, difficult tolerance(1e-6) ltolerance(1e-7)

You could also look at viewsource scobit.ado and viewsource scob_lf.ado to see the exact code used by StataCorp.

[Siddiq Niaz]

Dear Maarten,
Thank you for your valuable guidance and critical corrections to my preliminary code. Yes, I studied and referred to viewsource scob_lf.ado throughout my attempts of trying to replicate scobit. The scob_lf.ado codes the skewed logistic distribution as:

local f "((1+exp(`xb'))^(-`alpha'))"

However, in the manual, it states that F(z) should be:

1-1/{1+exp(xb)}^gamma   //Formula

Therefore, I used the later one as mentioned in Stata's scobit manual. The code (after your corrections) did not work for auto.dta, but it did converged for another random dataset and the results were pretty similar to the scobit output. So, thank you for putting me on right track.

Based on your guidance and corrections, I extended the binary choice code for ordered categories. The code i am trying is:

clear all
sysuse auto
recode rep78 (1 = 3) (2 = 3), gen(rep788) //Just for demonstration
capture program drop ordscob
program ordscob
version 14.2
args lnf xb alpha12 gamma1         //alpha12 are thresholds for ord prob & gamma1 is the param controlling skeweness in logit.
local f  "(1 - 1/(1+exp(`xb'))^(-`gamma1'))"
quietly replace `lnf' = ln(        `f') if $ML_y1 == 0 //for cateogry y = 0 outcome
quietly replace `lnf' = ln((`f'+`alpha12')- (`f')) if $ML_y1== 1 //for category y = 1 outcome
quietly replace `lnf' = ln(1-(`f'+`alpha12')) if $ML_y1== 2 //for category y = 3 outcome
end
//constraint define 1 gamma1:_cons=1                        //constraining gamma1 to be equal to one so as to reduce it to ordered logit.
//ml model lf ordscob (b: rep788 = trunk weight length)(alpha12:) /gamma1, constraint(1)
ml model lf ordscob (b: rep788 = trunk weight length)(alpha12:) /gamma1
ml check
ml maximize, difficult tolerance(1e-6) ltolerance(1e-7)

When i run the above code, Stata gives the following error:

could not find feasible values.

As a check, you may notice that i constrained the skewness parameter to equal 1 so as to test if the code reduces to and correctly estimates ordered logit or otherwise. Doing so results in the following warning and error:

"(note: constraint number 1 caused error r(198))" & "could not find feasible values."

I am aware of the convergence challenges in skewed logit (for binary or ordinal outcomes). However, for now, i am just worried if i am making some obvious error in coding the likelihood for skewed ordered choice model? Can you please review the code and guide me in a right direction?

Many thanks!

[Maarten Buis]

The formula in the documentation is the same as the formula in the .ado file as 1/x^a=x^-a

[Siddiq Niaz]

Thank you very much, Maarten.

If anyone else can guide regarding the extension to ordered choices, it will be highly appreciated.

[Eli Freeman]

Hi Siddiq,

I'm still learning this material myself. Try applying something like this:

args lnf xb alpha1 gamma1
quietly replace `lnf' = ln((1+exp(`xb'))^(-`gamma1')) if $ML_y1 == 0
/*contribution to log-likelihood for outcome y=0 respondents*/
quietly replace `lnf' = ln((1+exp(`xb'-`alpha1'))^(-`gamma1')-(1+exp(`xb'))^(-`gamma1')) if $ML_y1== 1
/*contribution to log-likelihood for outcome y=1 respondents*/
quietly replace `lnf' = ln(1-(1+exp(`xb'-`alpha1'))^(-`gamma1')) if $ML_y1== 2
/*contribution to log-likelihood for outcome y=2 respondents*/

[Eli Freeman]

Siddiq,

So I was working on a separate project and realized where some of the problem might be.

clear all
sysuse auto
recode rep78 (1 = 3) (2 = 3), gen(rep788) //Just for demonstration

tabulate rep788
//outcomes are 3, 4, 5

capture program drop ordscob
program ordscob
version 14.1
args lnf xb a1 a2      //alpha12 are thresholds for ord prob & gamma1 is the param controlling skeweness in logit.
 quietly replace `lnf' = ln((1+exp(`xb'-`a1'))^(-1)) if $ML_y1 == 3
/*contribution to log-likelihood for outcome y=0 respondents*/
 quietly replace `lnf' = ln((1+exp(`xb'-`a2'))^(-1)-(1+exp(`xb'-`a1'))^(-1)) if $ML_y1== 4
/*contribution to log-likelihood for outcome y=1 respondents*/
 quietly replace `lnf' = ln(1-(1+exp(`xb'-`a2'))^(-1)) if $ML_y1== 5
/*contribution to log-likelihood for outcome y=2 respondents*/
end
//constraint define 1 gamma1:_cons=1                        //constraining gamma1 to be equal to one so as to reduce it to ordered logit.
//ml model lf ordscob (b: rep788 = trunk weight length)(alpha12:) /gamma1, constraint(1)
ml model lf ordscob (b: rep788 = trunk weight length, nocons) /cut1 /cut2
ml maximize, difficult tolerance(1e-6) ltolerance(1e-7)

//if gamma=1 the model is identical to an ordered logit. If gamma<>1 then the distribution is asymmetric
//the ologit should match the results above
ologit rep788 trunk weight length


//only change, we are trying to estimate gamma
clear all
sysuse auto
recode rep78 (1 = 3) (2 = 3), gen(rep788) //Just for demonstration
capture program drop ordscob
program ordscob
version 14.1
args lnf xb a1 a2 gamma1      //alpha12 are thresholds for ord prob & gamma1 is the param controlling skeweness in logit.
 quietly replace `lnf' = ln((1+exp(`xb'-`a1'))^(-`gamma1')) if $ML_y1 == 3
/*contribution to log-likelihood for outcome y=0 respondents*/
 quietly replace `lnf' = ln((1+exp(`xb'-`a2'))^(-`gamma1')-(1+exp(`xb'-`a1'))^(-`gamma1')) if $ML_y1== 4
/*contribution to log-likelihood for outcome y=1 respondents*/
 quietly replace `lnf' = ln(1-(1+exp(`xb'-`a2'))^(-`gamma1')) if $ML_y1== 5
/*contribution to log-likelihood for outcome y=2 respondents*/
end
//constraint define 1 gamma1:_cons=1                        //constraining gamma1 to be equal to one so as to reduce it to ordered logit.
//ml model lf ordscob (b: rep788 = trunk weight length)(alpha12:) /gamma1, constraint(1)
ml model lf ordscob (b: rep788 = trunk weight length, nocons) /cut1 /cut2 /gamma1
ml init 0 0 0 5 5 1, copy
ml maximize, difficult tolerance(1e-6) ltolerance(1e-7)

//does not converge

[Behram Wali]

quietly replace `lnf' = ln((`f'+`alpha12')- (`f')) if $ML_y1== 1 //for category y = 1 outcome quietly replace `lnf' = ln(1-(`f'+`alpha12')) if $ML_y1== 2 //for category y = 3 outcome

Siddiq, there can be more problems in your code but it seems that you wrote down the wrong log-likelihood. Instead of having, e.g. ln(P(xb - alpha)), your code is ln(P(xb) + alpha). The correct way would be to:

quietly replace `lnf' = ln((`f'-`alpha12')- (`f')) if $ML_y1== 1 //for category y = 1 outcome
quietly replace `lnf' = ln(1-(`f'-`alpha12')) if $ML_y1== 2 //for category y = 3 outcome