st: Re: mixed logit with simulated ml

[Jenkins S P <stephenj@essex.ac.uk>]

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/