Re: st: Simulating Multinomial Logit in Stata

[Maarten Buis <maartenlbuis@gmail.com>]

On Mon, Mar 11, 2013 at 6:06 PM, Ali Hashemi wrote:
> I'm trying to simulate multinomial logit data in Stata. In this
> particular case, each individual has to pick from a set of 10
> alternatives defined by a combination of two attributes (x1 and x2).
> My question is regarding the way I should create the choice variable.
> Here is what I currently do:
>
> 1) Calculate the probability of each alternative as exp(XB)/(1+exp(XB))
> 1) Calculate the cumulative probability of each alternative (CDF) for
> each individual
> 2) Generate a random number (r)  for each individual (which is fixed
> across all the alternatives that he faces)
> 3) The first alternative with the cumulative probability greater than
> this random number would the choice. (CDF>r)
>
>
> I would appreciate your comment on the above procedure. Also, I'm not
> clear why don't we simply calculate the probability of each
> alternative as exp(XB)/(1+exp(XB)) and then pick the alternative with
> highest probability? Your help is much appreciated.

Below is an example of how to simulate a multinomial process in Stata.
You can obviously do what you proposed, but -mlogit- will not recover
those parameters, as it estimates a different model.

In the example I have also added some code on how I would analyse the
results of such a simulation. That part requires -simsum- (SSC)
-qplot- (SJ) and simpplot (SSC). You can use -findit- to find those
programs and install them.

*------------------ begin example ------------------
clear all
set more off

program define sim
	drop _all
	set obs 500
	gen x1 = rnormal()
	gen x2 = rnormal()

	gen xb1 = ln(.5) + ln(1.5)*x1 + ln(.5)*x2
	gen xb2 = ln(2) + ln(.33)*x1 + ln(2)*x2

	gen p1 = 1        / ( 1 + exp(xb1) + exp(xb2) )
	gen p2 = exp(xb1) / ( 1 + exp(xb1) + exp(xb2) )
	gen p3 = exp(xb2) / ( 1 + exp(xb1) + exp(xb2) )

	gen u = runiform()

	gen y = cond(u < p1, 1, ///
			cond(u < p1 + p2, 2, 3 ))

	mlogit y x1 x2, base(1) rrr
end

simulate  b1_2= _b[2:x1]  b2_2= _b[2:x2]  b0_2= _b[2:_cons] ///
         se1_2=_se[2:x1] se2_2=_se[2:x2] se0_2=_se[2:_cons] ///
          b1_3= _b[3:x1]  b2_3= _b[3:x2]  b0_3= _b[3:_cons] ///
         se1_3=_se[3:x1] se2_3=_se[3:x2] se0_3=_se[3:_cons] ///
         , reps(20000) : sim

// get an overview of how well -mlogit- performed		
simsum b1_2, se(se1_2) true(ln(1.5) ) mcse
simsum b2_2, se(se2_2) true(ln(0.5) ) mcse
simsum b0_2, se(se0_2) true(ln(0.5) ) mcse
simsum b1_3, se(se1_3) true(ln(0.33)) mcse
simsum b2_3, se(se2_3) true(ln(2)   ) mcse
simsum b0_3, se(se0_3) true(ln(2)   ) mcse

// take a closer look at the sampling distributions of the coefficients
qplot b1_2 b2_2 b0_2,   trscale(invnormal(@)) ///
    yline(`=ln(1.5)' `=ln(.5)'  ) xline(0)
qplot b1_3 b2_3 b0_3,   trscale(invnormal(@)) ///
    yline(`=ln(2)'    `=ln(.33)') xline(0)

// take a closer look at the p-values
gen p1_2 = 2*normal(- abs(( b1_2 - ln(1.5 ) ) / se1_2) )
gen p2_2 = 2*normal(- abs(( b2_2 - ln(0.5 ) ) / se2_2) )
gen p0_2 = 2*normal(- abs(( b0_2 - ln(0.5 ) ) / se0_2) )

gen p1_3 = 2*normal(- abs(( b1_3 - ln(0.33) ) / se1_3) )
gen p2_3 = 2*normal(- abs(( b2_3 - ln(2   ) ) / se2_3) )
gen p0_3 = 2*normal(- abs(( b0_3 - ln(2   ) ) / se0_3) )

gen id = _n
reshape long p1_ p2_ p0_, i(id) j(eq)

label define eq  2 "2 versus 1" ///
                 3 "3 versus 1"
label value eq eq
label var p0_ "constant"
label var p1_ "x1"
label var p2_ "x2"

simpplot p0_ p1_ p2_, by(eq)       ///
   scheme(s2color) legend(cols(3)) ///
   ylabel(,angle(horizontal))
*------------------- end example -------------------
(For more on examples I sent to the Statalist see:
http://www.maartenbuis.nl/example_faq )


---------------------------------
Maarten L. Buis
WZB
Reichpietschufer 50
10785 Berlin
Germany

http://www.maartenbuis.nl
---------------------------------
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/