Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

Antwort: st: Re: mixed logit with simulated ml


From   Arne Uhlendorff <auhlendorff@diw.de>
To   statalist@hsphsun2.harvard.edu
Subject   Antwort: st: Re: mixed logit with simulated ml
Date   Sun, 14 Aug 2005 17:15:54 +0200

Stephen,

thank you for your answer. I corrected that, but this was not the main 
problem. 

the main problem was that we did not add up the log likelihood 
contributions within one individual (persnr) per random draw before 
dividing it by the number of draws. I corrected that mistake and now the 
results are comparable to the results we get from gllamm. 


owner-statalist@hsphsun2.harvard.edu schrieb am 13.08.2005 11:15:43:

> Arne Uhlenforff <auhlendorff@diw.de> wants to estimate a MNL model 
> with a random intercept using simulated ML (rather than using -gllamm-).
> 
> I haven't looked at your program in detail, but one thing that strikes 
me 
> is that I think you are not treating the pseudo-random uniform draws 
> correctly.  You want 50 replications/draws per ML iteration, but are 
> creating them anew every iteration -- you shouldn't. They should be 
> created just once.  (See e.g. the gory details of the code underlying 
> -mvprobit- or associated SJ article.)  I don't know if this is having 
> knock-on effects elsewhere.
> 
> good luck
> Stephen
> 
> ------------------------------
> 
> Date: Fri, 12 Aug 2005 15:58:35 +0200
> From: Arne Uhlendorff <auhlendorff@diw.de>
> Subject: st: mixed logit with simulated ml
> 
> dear all,
> 
> 
> we try to estimate a multionmial logit model with a random intercept 
using
> simulated ml.
> 
> 
> Yet, our program has the problem that we never can identify unoberserved
> heterogeneity, although we know
> that it is present. (we estimated the model using gllamm)
> In other words,  the likelihood has the same value as when estimating 
with
> the normal mlogit command,
> whithout unobserved heterogeneity.
> 
> 
> Here is our program and the log file, do you have an idea what we are
> doing wrong?
> 
> thank you very much
> peter and arne
> 
> 
> sort persnr;
> qui by persnr: gen t=_n;
> label define tbylabel 1 "base"  2  "two"  3  "three";
> label values tby tbylabel;
> #d cr
> 
> mlogit tby sex, base(1)
> 
> matrix Init= e(b)
> 
> cap prog drop multinom_sim_d0
> program define multinom_sim_d0
>    args todo b lnf
>    tempvar etha2 etha3 random lj pi1 pi2 pi3 spi1 spi2 spi3 mpi1 mpi2 
mpi3
>    tempname lnsig sigma
>    mleval `etha2' = `b', eq(1)
>    mleval `etha3' = `b', eq(2)
>    mleval `lnsig' = `b', eq(3) scalar
> 
>    qui scalar `sigma'=(exp(`lnsig'))
>    qui gen double `random' = 0
> 
> qui  gen double `pi1'= 0
> qui  gen double `pi2'= 0
> qui  gen double `pi3'= 0
> qui  gen double `spi1'=0
> qui  gen double `spi2'=0
> qui  gen double `spi3'=0
> 
> set seed 123456789
> local repl=50
> local r=1
> while `r' <= `repl'  {
> qui replace `random' = (invnorm(uniform()))*`sigma' if t==1
> qui by persnr: replace `random' = `random'[1]
> 
> qui replace `pi1'= 1/(1 + exp(`etha2' + `random')+exp(`etha3' + 
`random'))
> qui replace `pi2'= exp(`etha2' + `random')*`pi1'
> qui replace `pi3'= exp(`etha3' + `random')*`pi1'
> 
> qui replace `spi1'=`spi1'+(`pi1')
> qui replace `spi2'=`spi2'+(`pi2')
> qui replace `spi3'=`spi3'+(`pi3')
> local r=`r'+1
>          }
> 
> qui gen double `mpi1'=`spi1'/`repl'
> qui gen double `mpi2'=`spi2'/`repl'
> qui gen double `mpi3'=`spi3'/`repl'
> 
>    qui gen double `lj' = `mpi1' if ($ML_y1==1)
>    qui replace `lj' = `mpi2' if ($ML_y1==2)
>    qui replace `lj' = `mpi3' if ($ML_y1==3)
>    qui mlsum `lnf'=ln(`lj')
>    if (`todo'==0|`lnf'>=.) exit
> 
> end
> 
> ml model d0 multinom_sim_d0     ( two: tby = sex )  /*
> */                                      ( three: tby = sex ) (baumi:)
> 
> 
> matrix start = (Init)
> *ml init start baumi:_cons=0.5
> ml maximize, difficult trace
> 
> 
> 
/***************************************************************************************************/
> LOG FILE
> 
> . mlogit tby sex, base(1)
> 
> Iteration 0:   log likelihood = -1349.1298
> Iteration 1:   log likelihood = -1332.0463
> Iteration 2:   log likelihood = -1331.9206
> Iteration 3:   log likelihood = -1331.9206
> 
> Multinomial logistic regression                   Number of obs   = 1313
>                                                    LR chi2(2)      = 
34.42
>                                                    Prob > chi2     = 
0.0000
> Log likelihood = -1331.9206                       Pseudo R2       = 
0.0128
> 
> - 
> 
------------------------------------------------------------------------------
>           tby |      Coef.   Std. Err.      z    P>|z|     [95% Conf.
> Interval]
> - 
> 
-------------+----------------------------------------------------------------
> two          |
>           sex |   .4227768   .1368457     3.09   0.002     .1545643
> .6909894
>         _cons |   .5381711   .0885475     6.08   0.000     .3646211
> .7117211
> - 
> 
-------------+----------------------------------------------------------------
> three        |
>           sex |   .9436537   .1632866     5.78   0.000     .6236178 
1.26369
>         _cons |  -.5460938   .1161788    -4.70   0.000    -.7737999
> - -.3183876
> - 
> 
------------------------------------------------------------------------------
> (Outcome tby==base is the comparison group)
> 
> .
> . matrix Init= e(b)
> 
> .
> . cap prog drop multinom_sim_d0
> 
> . program define multinom_sim_d0
>    1.   args todo b lnf
>    2.   tempvar etha2 etha3 random lj pi1 pi2 pi3 spi1 spi2 spi3 mpi1 
mpi2
> mpi3
>    3.   tempname lnsig sigma
>    4.   mleval `etha2' = `b', eq(1)
>    5.   mleval `etha3' = `b', eq(2)
>    6.   mleval `lnsig' = `b', eq(3) scalar
>    7.
> .   qui scalar `sigma'=(exp(`lnsig'))
>    8.   qui gen double `random' = 0
>    9.
> . qui  gen double `pi1'= 0
>   10. qui  gen double `pi2'= 0
>   11. qui  gen double `pi3'= 0
>   12. qui  gen double `spi1'=0
>   13. qui  gen double `spi2'=0
>   14. qui  gen double `spi3'=0
>   15.
> . set seed 123456789
>   16. local repl=50
>   17. local r=1
>   18. while `r' <= `repl'  {
>   19. qui replace `random' = (invnorm(uniform()))*`sigma' if t==1
>   20. qui by persnr: replace `random' = `random'[1]
>   21.
> . qui replace `pi1'= 1/(1 + exp(`etha2' + `random')+exp(`etha3' +
> `random'))
>   22. qui replace `pi2'= exp(`etha2' + `random')*`pi1'
>   23. qui replace `pi3'= exp(`etha3' + `random')*`pi1'
>   24.
> . qui replace `spi1'=`spi1'+(`pi1')
>   25. qui replace `spi2'=`spi2'+(`pi2')
>   26. qui replace `spi3'=`spi3'+(`pi3')
>   27. local r=`r'+1
>   28.         }
>   29.
> . qui gen double `mpi1'=`spi1'/`repl'
>   30. qui gen double `mpi2'=`spi2'/`repl'
>   31. qui gen double `mpi3'=`spi3'/`repl'
>   32.
> .   qui gen double `lj' = `mpi1' if ($ML_y1==1)
>   33.   qui replace `lj' = `mpi2' if ($ML_y1==2)
>   34.   qui replace `lj' = `mpi3' if ($ML_y1==3)
>   35.   qui mlsum `lnf'=ln(`lj')
>   36.   if (`todo'==0|`lnf'>=.) exit
>   37.
> . end
> 
> .
> . ml model d0 multinom_sim_d0     ( two: tby = sex )  /*
> > */                                      ( three: tby = sex ) (baumi:)
> 
> .
> .
> . matrix start = (Init)
> 
> . *ml init start baumi:_cons=0.5
> . ml maximize, difficult trace
> 
> initial:       log likelihood = -1461.4353
> trying nonzero initial values ++
> alternative:   log likelihood = -1438.4468
> rescaling entire vector ..
> rescale:       log likelihood = -1438.4468
> rescaling equations .+.+.++++++.
> rescaling equations ......
> rescale eq:    log likelihood = -1353.8803
> - 
> 
------------------------------------------------------------------------------
> Iteration 0:
> Coefficient vector:
>           two:      two:    three:    three:    baumi:
>           sex     _cons       sex     _cons     _cons
> r1         0         1         0       .25  .0078125
> 
>                                                     log likelihood =
> - -1353.8803
>                                                                   (not
> concave)
> - 
> 
------------------------------------------------------------------------------
> Iteration 1:
> Coefficient vector:
>            two:       two:     three:     three:     baumi:
>            sex      _cons        sex      _cons      _cons
> r1   .5790007   .2912365    1.10536  -.7637936  -.9995876
> 
>                                                     log likelihood =
> - -1339.133
> - 
> 
------------------------------------------------------------------------------
> ...
> - 
> 
------------------------------------------------------------------------------
> Iteration 10:
> Coefficient vector:
>            two:       two:     three:     three:     baumi:
>            sex      _cons        sex      _cons      _cons
> r1   .4227916   .5382313   .9436733  -.5460403  -14.37365
> 
> numerical derivatives are approximate
> nearby values are missing
>                                                     log likelihood =
> - -1331.9206
> - 
> 
------------------------------------------------------------------------------
> 
>                                                    Number of obs   = 
1313
>                                                    Wald chi2(1)    = 
9.55
> Log likelihood = -1331.9206                       Prob > chi2     = 
0.0020
> 
> - 
> 
------------------------------------------------------------------------------
>               |      Coef.   Std. Err.      z    P>|z|     [95% Conf.
> Interval]
> - 
> 
-------------+----------------------------------------------------------------
> two          |
>           sex |   .4227916   .1368477     3.09   0.002      .154575
> .6910082
>         _cons |   .5382313   .0885484     6.08   0.000     .3646797
> .7117829
> - 
> 
-------------+----------------------------------------------------------------
> three        |
>           sex |   .9436733   .1632881     5.78   0.000     .6236345
> 1.263712
>         _cons |  -.5460403   .1161791    -4.70   0.000    -.7737472
> - -.3183333
> - 
> 
-------------+----------------------------------------------------------------
> baumi        |
>         _cons |  -14.37365   684.2106    -0.02   0.983    -1355.402
> 1326.654
> - 
> 
------------------------------------------------------------------------------
> 
> 
> Stephen
> =============================================
> Professor Stephen P. Jenkins <stephenj@essex.ac.uk>
> Institute for Social and Economic Research (ISER)
> University of Essex, Colchester CO4 3SQ, UK
> Phone: +44 1206 873374.  Fax: +44 1206 873151.
> http://www.iser.essex.ac.uk
> 
> *
> *   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/
> 

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



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index