Help with maximum log likelihood estimator

[Bruno Nobrega]

I'm trying to replicate a research that proposed the following logistic equation to estimate travel mode choice.

From what i've read, I should be using the "program" command to describe my equation and then use the model maximize to estimate the thetas. However, I've been trying to run the following code, but the model does not converge. It starts with a positive log-likelihood and when it maximizes it starts growing to infinite.
Is there a better way of estimating thetas?
Did I do something wrong?

program modoescol
 version 13
 args lnf theta1 theta2 theta3 theta4
 tempvar c o d e
 quietly gen double `c' = exp(-`theta1')
 quietly gen double `o'= exp(-`theta2')
 quietly gen double `d' = `theta3'
 quietly gen double `e' = `theta4'
 quietly replace `lnf' = `d'/(1+`d') + ((`c'/(`c'+`o')) *(1/(1+`d')))
 
 end
ml model lf modoescol (modo = custoporrenda_carro Tempo_carro ) (custoporrenda_onibus Tempo_ônibus ) () ()
ml check
ml search
ml maximize

the dataset is somewhat large, but if it helps I made a subset of the data for testing

* Example generated by -dataex-. For more info, type help dataex
clear
input int morador_id byte modo double(Tempo_carro curtoporrenda_carro Tempo_ônibus curtoporrenda_onibus)
  2 0 18.71666666666667  .03614832236531986                42 .0011515151515151514
11 1 25.98333333333333  .03577375505892256             55.35 .0011515151515151514
12 1 31.76666666666667  .03547564759259259 92.33333333333333 .0033030303030303033
19 1               4.4  .02213177321212121             10.75 .0006909090909090909
26 1              17.1  .18115827499999998             33.65  .005757575757575757
27 0              35.9  .03526259095959596 127.7333333333333 .0036363636363636364
33 1 14.83333333333333  .18174246254208754 40.61666666666667  .008333333333333333
37 1 10.43333333333333  .10972588419191918 49.86666666666667  .007454545454545454
52 1 14.06666666666667  .10916403323232322 35.46666666666667                 .005
53 1 15.71666666666667  .06806805045770202             40.15              .003125
61 1 4.766666666666667  .06912635344065657 24.31666666666667 .0015340909090909092
62 1 9.316666666666666  .06868660197285353 22.13333333333333              .003125
63 0 26.46666666666667  .06702907720959596 89.86666666666666              .003125
75 0              31.6  .10645271575757574 77.88333333333334  .008454545454545456
83 1 19.23333333333333  .02167301417171717             40.25 .0006909090909090909
99 1             21.75  .10797589886363634 72.81666666666666  .007454545454545454
100 1             12.65  .10938310356060606 60.16666666666666  .005909090909090909
109 0             34.25  .10604292537878787              64.8                 .005
110 1 8.383333333333333  .11004289184343434             39.35  .008454545454545456
136 1 6.983333333333333 .015751340696248197 13.18333333333333 .0004935064935064935
139 0             27.55  .03569299972222222              68.7 .0016666666666666668
144 0 24.41666666666667 .035854510395622895 58.93333333333333  .002818181818181818
146 1              10.9  .06853357490530303             27.45 .0015340909090909092
148 1 32.38333333333333  .06645723922032828 69.66666666666667  .005284090909090909
149 0 16.68333333333333 .021751879489898988 38.26666666666667                 .001
150 1 18.61666666666667 .021692086176767676              52.4 .0006909090909090909
158 0               6.2  .18396752979797978              17.1  .008333333333333333
162 0 12.13333333333333   .1094629997979798              29.7  .008454545454545456
165 0              7.85  .36708455126262624 31.23333333333333   .02303030303030303
167 0 44.28333333333333   .0348304639983165 89.06666666666666  .002818181818181818
171 0 12.26666666666667   .0364807938047138 17.86666666666667 .0016666666666666668
175 0 24.88333333333333  .06718210427714646 49.96666666666667              .003125
193 1              18.9 .021683323363636365 60.68333333333333  .001690909090909091
196 0              24.8  .06719015833333332             63.25              .003125
222 0             31.65  .06652811491477272 70.08333333333333  .005284090909090909
223 1 24.73333333333333   .1791909375420875 71.08333333333333  .008333333333333333
226 0 26.31666666666667  .06704357451073231             55.95  .005284090909090909
236 0              27.1  .06696786638257576 57.58333333333334   .00840909090909091
242 0              15.4  .18159641565656565 35.28333333333333  .004090909090909091
263 1 19.93333333333333  .03608560811447811 69.76666666666667  .002818181818181818
264 1 17.21666666666667 .036225641304713806              44.8 .0011515151515151514
266 1              23.2  .06734479621212121             45.95  .005284090909090909
274 0               3.4  .06925843996212121 10.26666666666667  .002159090909090909
277 0 8.516666666666667  .06876392091224746             21.25  .005284090909090909
284 1 31.78333333333333  .06651522842487373 63.63333333333333              .003125
285 0 34.98333333333333  .06620595266729798 73.01666666666667               .00625
286 1 12.36666666666667  .06839182351641414 26.03333333333333  .002159090909090909
287 1             29.65  .06672141226325758             69.75               .00625
291 1             35.15  .06618984455492424             76.85  .005284090909090909
299 1 13.36666666666667  .06829517484217172 25.58333333333333  .002159090909090909
306 0 33.51666666666667  .06634770405618687 66.48333333333333              .003125
310 0 15.66666666666667  .06807288289141414 37.76666666666667  .002159090909090909
312 0 26.16666666666667  .06705807181186868 62.66666666666666              .003125
315 1 31.11666666666667  .06657966087436869              62.8  .005284090909090909
373 1 9.266666666666667  .03663543168350168 18.93333333333333 .0016666666666666668
378 1              12.9  .10934444409090908              27.9 .0034545454545454545
379 0             10.55  .18284640517676767 27.03333333333333  .005757575757575757
393 1 17.16666666666667   .0362282186026936 36.66666666666666 .0011515151515151514
394 1 38.16666666666666  .03514575345117845 99.96666666666667  .004151515151515152
397 1 11.43333333333333  .03652374877104377              21.2 .0016666666666666668
402 1 3.566666666666667  .06924233184974747 12.91666666666667  .002159090909090909
403 1 18.06666666666667   .0678409260732323 58.23333333333333  .005284090909090909
405 1 4.933333333333334  .06911024532828282 28.21666666666667  .004318181818181818
406 0 33.71666666666667  .06632837432133838 65.18333333333334  .004659090909090909
407 0 12.76666666666667  .06835316404671717 31.56666666666667  .002159090909090909
413 0 22.53333333333333 .035951588619528616 55.43333333333333 .0024848484848484847
414 0             15.75  .03630124204545455             48.25 .0016666666666666668
416 0              32.6  .03543269262626263 80.03333333333333  .002818181818181818
419 0              12.4    .036473921010101 32.78333333333333 .0011515151515151514
421 1             15.15  .03633216962121212             35.05 .0011515151515151514
431 1              10.4  .10973103878787878 30.01666666666667 .0034545454545454545
434 0 14.03333333333333  .10916918782828282             42.55                 .005
436 1              14.9  .18172528055555554              30.5  .005757575757575757
437 0 26.43333333333333  .06703229883207071 49.61666666666667              .003125
438 0 32.18333333333333  .06647656895517677             58.75              .003125
439 0 33.31666666666667  .06636703379103534 56.51666666666667              .003125
443 1 17.23333333333333    .067921466635101 39.18333333333333  .002159090909090909
469 0              29.5                   .              61.6                    .
470 0              12.2                   . 37.36666666666667                    .
471 0 14.63333333333333                   .             35.75                    .
474 1 33.53333333333333  .03538458306397306 99.96666666666667  .004484848484848485
475 1 16.28333333333333  .10882125260101008 49.88333333333333                 .005
477 1 5.133333333333334  .02210909298989899             14.55 .0006909090909090909
486 1 23.06666666666667   .0673576827020202              53.4              .003125
503 1             11.35  .02191682656060606              18.4                 .001
518 1              11.2  .03653577616161616 23.68333333333333 .0011515151515151514
522 0 31.86666666666667  .17735246498316498 73.88333333333334  .014090909090909091
525 1 20.36666666666667  .18031635765993265 44.36666666666667  .005757575757575757
540 1              15.7  .10891145803030303 28.41666666666667 .0034545454545454545
543 1             35.55  .10584189613636363 81.13333333333334                 .005
548 0              10.4  .18288506464646462              23.1  .005757575757575757
561 0              21.1  .18012735580808079 65.51666666666667   .01818181818181818
562 1             30.65  .06662476358901515 87.88333333333334  .005284090909090909
563 1 31.98333333333333  .06649589869002524 70.18333333333334              .003125
564 0 26.18333333333333  .06705646100063131             52.95  .002159090909090909
568 0 26.96666666666667  .35723068198653196             77.75  .024848484848484845
572 0             23.35  .17954746376262626 55.71666666666667  .016666666666666666
574 1 25.08333333333333 .035820146422558916              57.7 .0033333333333333335
579 0 18.91666666666667  .18069006586700334 33.83333333333334  .005757575757575757
end

Thanks in advance

[Joao Santos Silva]

Bruno Nobrega,

You are trying to maximize the probability, not the log-likelihood function (which is a function of the probability).

Best wishes,

Joao

[Joseph Coveney]

I'm trying to replicate a research that proposed the following logistic equation to estimate travel mode choice. . . . Is there a better way of estimating thetas?

You could try something like the following.

I used a Student's t prior for everything in a preliminary examination, but if you want maximum likelihood estimates with the full dataset, then you can substitute a flat prior.

Complete do-file and log file attached for your convenience: be sure to check that I got the equation in correctly.

version 18.0

clear *

quietly input int morador_id byte modo double(Tempo_carro curtoporrenda_carro Tempo_ônibus curtoporrenda_onibus)
<redacted for brevity>
end

#delimit ;
rename (morador_id modo Tempo_carro curtoporrenda_carro Tempo_ônibus
    curtoporrenda_onibus) (mid out tca cpc tbu cpb);

bayesmh out = ({d} / (1 + {d}) + (1 / (1 + {d})) *
        (exp({drf: tca cpc, xb})  / ((exp({drf:})) + exp({tra: tbu cpb, xb})))),
    likelihood(logistic)
    prior({d} {drf:tca cpc _cons} {tra: tbu cpb _cons}, t(0, 2.5, 7))
    rseed(1842550798) thin(10);
#delimit cr

set scheme s2color
local o ylabel( , angle(horizontal) nogrid)
local i 1
foreach type in trace ac histogram {
    foreach parm in "drf:tca" "drf:cpc" "drf:_cons" "tra:tbu" "tra:cpb" ///
            "tra:_cons" d {
        bayesgraph `type' {`parm'}, `o' name(A`i')
        local graphs `graphs' A`i'
        local ++i
    }
}
graph combine `graphs', rows(7) colfirst ysize(5.5) xsize(4)
quietly graph export driving.png

exit

I thinned in order to prettify the caterpillars and autocorrelograms.

[Joseph Coveney]

be sure to check that I got the equation in correctly.

Actually, in the way I put that equation in the syntax it defines the linear predictor of a logistic (generalized linear) regression model. As the embedded figure indicates, and as Joao reminds, it is intended to define the probability per se.

So, you'll either have to back-transform the parameters somehow, or else model the probability directly, for example, as in the following modification.

bayesmh out,
    likelihood(dbinomial(({d} / (1 + {d}) + (1 / (1 + {d})) *
        (exp({drf: tca cpc, xb})  / ((exp({drf:})) + exp({tra: tbu cpb, xb})))), 1))
    prior({d} {drf:tca cpc _cons} {tra: tbu cpb _cons}, t(0, 2.5, 7));

[Joseph Coveney]

. . . the model does not converge. It starts with a positive log-likelihood and when it maximizes it starts growing to infinite.

In that equation whose image you embedded in your post, do those two exponentiated linear predictors have constants (intercepts)? Or is d intended to represent a common constant for the equation?

When I set up the maximum likelihood estimator (see code below) to fit constants (intercepts) for each of the two linear predictors, it requires about 600 iterations to converge and the coefficient for d approaches -1 as if it's not really identified.


ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿobsÿ=ÿÿÿÿÿ96
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿWaldÿchi2(2)ÿÿ=ÿÿ19.67
Logÿlikelihoodÿ=ÿ-47.824802ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿProbÿ>ÿchi2ÿÿÿ=ÿ0.0001

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿoutÿ|ÿCoefficientÿÿStd.ÿerr.ÿÿÿÿÿÿzÿÿÿÿP>|z|ÿÿÿÿÿ[95%ÿconf.ÿinterval]
-------------+----------------------------------------------------------------
carÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿtcaÿ|ÿÿ-.1079771ÿÿÿ.0268219ÿÿÿÿ-4.03ÿÿÿ0.000ÿÿÿÿ-.1605469ÿÿÿ-.0554072
ÿÿÿÿÿÿÿÿÿccaÿ|ÿÿÿ5.408441ÿÿÿ2.980806ÿÿÿÿÿ1.81ÿÿÿ0.070ÿÿÿÿ-.4338306ÿÿÿÿ11.25071
ÿÿÿÿÿÿÿ_consÿ|ÿÿ-5.988818ÿÿÿ2.024793ÿÿÿÿ-2.96ÿÿÿ0.003ÿÿÿÿÿ-9.95734ÿÿÿ-2.020296
-------------+----------------------------------------------------------------
busÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿtbuÿ|ÿÿ-.0486883ÿÿÿ.0137139ÿÿÿÿ-3.55ÿÿÿ0.000ÿÿÿÿ-.0755671ÿÿÿ-.0218095
ÿÿÿÿÿÿÿÿÿcbuÿ|ÿÿÿ424.4245ÿÿÿ94.15121ÿÿÿÿÿ4.51ÿÿÿ0.000ÿÿÿÿÿ239.8915ÿÿÿÿ608.9575
ÿÿÿÿÿÿÿ_consÿ|ÿÿ-20.84686ÿÿÿ2.030301ÿÿÿ-10.27ÿÿÿ0.000ÿÿÿÿ-24.82618ÿÿÿ-16.86755
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿ/dÿ|ÿÿ-.9999964ÿÿÿ6.08e-06ÿ-1.6e+05ÿÿÿ0.000ÿÿÿÿ-1.000008ÿÿÿ-.9999845
------------------------------------------------------------------------------

But when I remove the two constants (intercepts) from the exponentiated linear predictors and allow d to be the sole constant for the nonlinear equation, then it converges quickly to what seems a more reasonable maximum for d. (You can remove the constants by using the noconstant option for each equation when setting up ml model or you can remove them by setting up constraints. I chose the latter here.)


ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿobsÿ=ÿÿÿÿÿ96
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿWaldÿchi2(2)ÿÿ=ÿÿÿ4.73
Logÿlikelihoodÿ=ÿ-61.536976ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿProbÿ>ÿchi2ÿÿÿ=ÿ0.0941

ÿ(ÿ1)ÿÿ[car]_consÿ=ÿ0
ÿ(ÿ2)ÿÿ[bus]_consÿ=ÿ0
------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿoutÿ|ÿCoefficientÿÿStd.ÿerr.ÿÿÿÿÿÿzÿÿÿÿP>|z|ÿÿÿÿÿ[95%ÿconf.ÿinterval]
-------------+----------------------------------------------------------------
carÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿtcaÿ|ÿÿÿÿ-.12178ÿÿÿ.0613992ÿÿÿÿ-1.98ÿÿÿ0.047ÿÿÿÿ-.2421202ÿÿÿ-.0014398
ÿÿÿÿÿÿÿÿÿccaÿ|ÿÿÿ7.434001ÿÿÿÿ7.65093ÿÿÿÿÿ0.97ÿÿÿ0.331ÿÿÿÿ-7.561545ÿÿÿÿ22.42955
ÿÿÿÿÿÿÿ_consÿ|ÿÿÿÿÿÿÿÿÿÿ0ÿÿ(omitted)
-------------+----------------------------------------------------------------
busÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿtbuÿ|ÿÿ-.0516099ÿÿÿ.0241575ÿÿÿÿ-2.14ÿÿÿ0.033ÿÿÿÿ-.0989577ÿÿÿ-.0042621
ÿÿÿÿÿÿÿÿÿcbuÿ|ÿÿÿ225.3094ÿÿÿ139.3781ÿÿÿÿÿ1.62ÿÿÿ0.106ÿÿÿÿÿ-47.8666ÿÿÿÿ498.4854
ÿÿÿÿÿÿÿ_consÿ|ÿÿÿÿÿÿÿÿÿÿ0ÿÿ(omitted)
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿ/dÿ|ÿÿÿ.2127019ÿÿÿ.3351028ÿÿÿÿÿ0.63ÿÿÿ0.526ÿÿÿÿ-.4440875ÿÿÿÿ.8694913
------------------------------------------------------------------------------

The evaluator program for the equation with a binomial likelihood is given below.

version 18.0

clear *

quietly input int morador_id byte modo ///
    double(Tempo_carro curtoporrenda_carro Tempo_ônibus curtoporrenda_onibus)
<redacted for brevity?
end

*
* Begin here
*
rename ///
    (modo Tempo_carro curtoporrenda_carro Tempo_ônibus curtoporrenda_onibus) ///
    (out tca cca tbu cbu)

program define modoescol
    version 18.0
    args lnf theta1 theta2 theta3 theta4

    tempvar c o d e
    generate double `c' = exp(`theta1')
    generate double `o' = exp(`theta2')
    generate double `d' = `theta3'

    generate double `e' = `d' / (1 + `d') + `c' / (`c'+`o') / (1 + `d')

    quietly replace `lnf' = ln(`e') if $ML_y1
    quietly replace `lnf' = ln(1 - `e') if !$ML_y1

end

constraint define 1 _b[car:_cons] = 0
constraint define 2 _b[bus:_cons] = 0
ml model lf modoescol (car: out = tca cca) ///
    (bus: tbu cbu) (d:, freeparm), constraints(1/2)
ml check
ml search
ml maximize // , nolog

exit

Do-file and log file (constrained model as shown above) attached.