Bootstrap standard error of marginal effect for mixlogit

[Aurelie Bonein]

Dear all stata users,

I have data from a stated choice experiment with 162 individuals who answered from 5 sets of situations each.
Each set of situation requires from individuals to choose between a contract A, a contract B or the reference contract.
In each contract, there are 3 attributes : price, duration and loss. The level of attributes vary across contract.

I have conducted mixlogit regressions. To do so, I have used the User-written command mixlogit to fit mixed logit models (Hole, 2007). I have used ssc install mixlogit- to install -mixlogit; next I have computed the marginal effects using the mixlpred command.
I have read that it is possible to bootstrap the marginal effect in order to get the standard error.
I have read the help about the bootstrap command but I don't achieve to implement it after the mixlpred command.
TO compute the marginal effect, I have this command :

mixlogit decision asc level_price , group(id_occas) rand(dureeQ dureeQ2 perteQ) id(id) nrep(100)
mixlpred prep_base
preserve
quietly replace dureeQ=dureeQ+1 if decision==1
mixlpred prep_inc
generate diff=prep_inc - prep_base if decision==1
sum diff if decision==1
restore
With this command I get the mean of the marginal effect. How can I bootstrap it ?
I know that a usual bootstrap command is like this :
bootstrap (I don't know what to enter here),cluster(id) idcluster(newid) group(id_occas) reps(10): mixlogit decision asc level_price, group(id_occas) id(newid) rand(dureeQ dureeQ2 perteQ)

If someone can help me, I will greatly appreciate
Thank you

[Arne Risa Hole]

Aurelie

Here is an example. Note that this takes a long time to run even with only 10 bootstrap replications, which is why the bootstrap is typically only practical to implement with small datasets in this context.

use "http://fmwww.bc.edu/repec/bocode/t/traindata.dta"

* Create alternative id
set seed 12345
gen rnd = runiform()
bysort pid gid (rnd): gen alt = _n

capture program drop mixl_bstrap
program mixl_bstrap, eclass

    mixlogit y price, rand(contract local wknown tod seasonal) ///
    group(gid) id(newpid) nrep(500)

    preserve
    replace wknown = 0 if alt==1
    mixlpred p0, nrep(500)
    replace wknown = 1 if alt==1
    mixlpred p1, nrep(500)
    gen p_diff = p1-p0
    sum p_diff if alt==1
    restore

    matrix meff = r(mean)

    ereturn post meff
    ereturn local cmd "bootstrap"
end

bootstrap, cluster(pid) idcluster(newpid) group(gid) reps(10): mixl_bstrap

Arne

[Aurelie Bonein]

It works !
Thank you very much for your help
Best
Aurélie

[Arik Beremzon]

Dear Dr. Hole,

I am using a similar approach as you have provided above, yet I keep getting an error. The author of the script is my supervisor and we're both stuck. It does, however, work with his research.

My sample size is 162, 12 questions with 2 alternatives. Over 3000 observations total.

The error is:

insufficient observations to compute bootstrap standard errors
no results will be saved
r(2000);

capture program drop ame_mixlogit_sc3
program define ame_mixlogit_sc3, rclass

capture drop p0* p1* p_diff*
capture scalar drop diff*

capture preserve

mixlogit DCE, group(groupid) id(ResponseID) rand(price6 price10 time30 ///
    time50 taste_good taste_vgood health_neutral healthy)

global vars price6 price10 time30 time50 taste_good taste_vgood health_neutral healthy
    
foreach i in $vars {
    capture preserve
    replace `i'=0 if alt==1
    mixlpred p0`i', nrep(500)
    replace `i'=1 if alt==1
    mixlpred p1`i', nrep(500)
    gen p_diff`i' = p1`i' - p0`i'
    sum p_diff`i' if alt==1
    return sca diff_`i' = r(mean)
    capture restore
    }
capture restore
end
**************************************************************************************************
*            END PROGRAM AME_MIXLOGIT_SC3                                                         *
**************************************************************************************************
gen newid = ResponseID
set seed 1
gen rnd=runiform()
bysort ResponseID choiceset (rnd): gen alt = _n

*RUN BOOTSTRAP
capture preserve

bootstrap ame_price6=r(diff_price6) ame_price10=r(diff_price10) ame_time30=r(diff_time30) ame_time50=r(diff_time50) ///
    ame_taste_good=r(diff_taste_good) ame_taste_vgood=r(diff_taste_vgood) ame_health_neutral=r(diff_health_neutral) ///
    ame_healthy=r(diff_healthy), reps(1000) cluster(ResponseID) idcluster(newid) group(choiceset): ame_mixlogit_sc3
//estimates store est_sc3

capture restore

[Arne Risa Hole]

Arik

I suspect the error is related to the specification of the cluster(), idcluster() and group() options - see help bootstrap for more info about these options. Without knowing anything about your data, my guess is that the code should be:

capture program drop ame_mixlogit_sc3
program define ame_mixlogit_sc3, rclass

capture drop p0* p1* p_diff*
capture scalar drop diff*

capture preserve

mixlogit DCE, group(groupid) id(newid) rand(price6 price10 time30 ///
    time50 taste_good taste_vgood health_neutral healthy)

global vars price6 price10 time30 time50 taste_good taste_vgood health_neutral healthy
    
foreach i in $vars {
    capture preserve
    replace `i'=0 if alt==1
    mixlpred p0`i', nrep(500)
    replace `i'=1 if alt==1
    mixlpred p1`i', nrep(500)
    gen p_diff`i' = p1`i' - p0`i'
    sum p_diff`i' if alt==1
    return sca diff_`i' = r(mean)
    capture restore
    }
capture restore
end
**************************************************************************************************
*            END PROGRAM AME_MIXLOGIT_SC3                                                         *
**************************************************************************************************
gen newid = ResponseID
set seed 1
gen rnd=runiform()
bysort ResponseID choiceset (rnd): gen alt = _n

*RUN BOOTSTRAP
capture preserve

bootstrap ame_price6=r(diff_price6) ame_price10=r(diff_price10) ame_time30=r(diff_time30) ame_time50=r(diff_time50) ///
    ame_taste_good=r(diff_taste_good) ame_taste_vgood=r(diff_taste_vgood) ame_health_neutral=r(diff_health_neutral) ///
    ame_healthy=r(diff_healthy), reps(1000) cluster(ResponseID) idcluster(newid) group(groupid): ame_mixlogit_sc3
//estimates store est_sc3

capture restore

I am also not sure why you need all the capture preserve/capture restore commands, but I don't think they are causing the error.

Arne

[Anderson Tami]

Dear all,
I am having problems estimating the standard errors of the marginal effects using the mixlogit command. First of all, I want to explain what I am trying to do with the model. I want to estimate the socio economic factors that determine that a young person is on one of the next four situations: only attending to higher education (variable eduypar is then 2), only participating in the labor market (variable eduypar is then 3), doing both simultaneously (variable eduypar is then 4), and doing neither of both (variable eduypar is then 1).
Some of the independetnt variables are age (variable name is edad),
male (variable name is hombre),
labor experience as employed (variable name is exlabtot),
experience working by himself (variable name is excptot)
experience working with family without wage (variable name is exfamsinremtot)
dummies for schooling level of the father (variable names are edupadre2, edupadre3 and edupadre4)
dummies for schooling level of the mother (variable names are edumadre2, edumadre3 and edumadre4)
dummies for ocupation of the father (variable names are ocupadre1, ocupadre2, ocupadre3, and ocupadre5)
dummies for ocupation of the mother (variable names are ocumadre1, ocumadre2, ocumadre3, and ocumadre5)
a dummy if member of the family have tried (or are working) to be employed using the contacs with friends (relfam1)
a dummy if the young person lived before in another city (ep6074)
As i think that IIA assumption is not achieved, i want to estimate a mixed multinomial logit model. (This a cross sectional analysis)

My commands are
**
gen id=_n
expand 4
by id, sort: gen alt = _n
gen neet=(alt==1)
gen se = (alt == 2)
gen st = (alt == 3)
gen eyt = (alt == 4)
gen choice=(eduypar==alt)
global alist edad hombre qm qa edupadre2 edupadre3 edupadre4 edumadre2 edumadre3 edumadre4 ocupadre1 ocupadre2 ocupadre3 ocupadre5 ocumadre1 ocumadre2 ocumadre3 ocumadre5 excptot exfamsinremtot exlabtot relfam1 ep6074
foreach var of varlist $alist {
gen `var'_neet=`var'*neet
}
foreach var of varlist $alist {
gen `var'_se=`var'*se
}
foreach var of varlist $alist {
gen `var'_st=`var'*st
}
foreach var of varlist $alist {
gen `var'_eyt=`var'*eyt
}
mixlogit choice edad_se edad_st edad_eyt qm_se qm_st qm_eyt qa_se qa_st qa_eyt hombre_se hombre_st hombre_eyt edupadre2_se edupadre2_st edupadre2_eyt edupadre3_se edupadre3_st edupadre3_eyt edupadre4_se edupadre4_st edupadre4_eyt edumadre2_se edumadre2_st edumadre2_eyt edumadre3_se edumadre3_st edumadre3_eyt edumadre4_se edumadre4_st edumadre4_eyt ocupadre1_se ocupadre1_st ocupadre1_eyt ocupadre2_se ocupadre2_st ocupadre2_eyt ocupadre3_se ocupadre3_st ocupadre3_eyt ocupadre5_se ocupadre5_st ocupadre5_eyt ocumadre1_se ocumadre1_st ocumadre1_eyt ocumadre2_se ocumadre2_st ocumadre2_eyt ocumadre3_se ocumadre3_st ocumadre3_eyt ocumadre5_se ocumadre5_st ocumadre5_eyt excptot_se excptot_st excptot_eyt exfamsinremtot_se exfamsinremtot_st exfamsinremtot_eyt exlabtot_se exlabtot_st exlabtot_eyt relfam1_se relfam1_st relfam1_eyt ep6074_se ep6074_st ep6074_eyt, group(id) rand(se st eyt)

****

First, it was possible to obtain results only with a group variable (Group(id)) and without an id variable (id()) in the mixlogit command
The next thing i did was to obtain the marginal effects of each variable (It was troublesome, but i did it). Then I proved the bootstrap command with only one of the marginal effects
the command was the next. the variable mfxdep6074_eyt is the marginal effect of migration when the alternative is 4.

**

bootstrap mfxdep6074_eyt , reps(10): mixlogit choice edad_se edad_st edad_eyt qm_se qm_st qm_eyt qa_se qa_st qa_eyt hombre_se hombre_st hombre_eyt edupadre2_se edupadre2_st edupadre2_eyt edupadre3_se edupadre3_st edupadre3_eyt edupadre4_se edupadre4_st edupadre4_eyt edumadre2_se edumadre2_st edumadre2_eyt edumadre3_se edumadre3_st edumadre3_eyt edumadre4_se edumadre4_st edumadre4_eyt ocupadre1_se ocupadre1_st ocupadre1_eyt ocupadre2_se ocupadre2_st ocupadre2_eyt ocupadre3_se ocupadre3_st ocupadre3_eyt ocupadre5_se ocupadre5_st ocupadre5_eyt ocumadre1_se ocumadre1_st ocumadre1_eyt ocumadre2_se ocumadre2_st ocumadre2_eyt ocumadre3_se ocumadre3_st ocumadre3_eyt ocumadre5_se ocumadre5_st ocumadre5_eyt excptot_se excptot_st excptot_eyt exfamsinremtot_se exfamsinremtot_st exfamsinremtot_eyt exlabtot_se exlabtot_st exlabtot_eyt relfam1_se relfam1_st relfam1_eyt ep6074_se ep6074_st ep6074_eyt, group(id) rand(se st eyt)

***

After this, stata (I am using stata 14) shows the comment
"

Bootstrap replications (10)
1 ---+--- 2 ---+---	3 ---+-	- 4 ---+-	- 5
xxxxxxxxxx
insufficient observations to	compute	bootstrap	standard	errors
no results will be saved

"
I don´t know how to resolve this problem

I hope anyone can help me
Thank you
Anderson Tami

[Emmanuel Teitelbaum]

Arne Risa Hole

Thanks for posting this example and code! Recently, a Stata representative showed me a way to perform the same calculations using the new cmxtmixlogit command (which I will share below), but we were unable to agree on the interpretation of the results in the context of an experiment with unlabeled alternatives. Thus, my main question is:

Can we say that the marginal effect represents the average effect of the attribute level on the probability of choosing X concept? So to draw on the example that you posted here using the Train data, the quality of being well-known increases (on average) the likelihood that a respondent will choose an electricity provider by 13.47%?

I also have a couple of questions about the estimation of the marginal effects in your example:

1) I was wondering about the reason for randomizing the order of the attributes before calculating the marginal effects. Was this done to avoid order bias or something along those lines?

2) How would you handle a "none" option in this context? Let's say you had three unlabeled alternatives and a "none" option. Would you randomize the none option along with the other alternatives and calculate the marginal effects from there? Or would you randomize the three unlabeled alternatives, leaving the none option out of the marginal effects calculation?

3) Using this method, we are calculating the marginal effect for one of the alternatives, but leaving the other three out. I have an intuition that says we should make use of all of the data in the calculation. Is there a way to take an average of the marginal effects across all of the alternatives? Or is that simply unnecessary, or impossible?

As promised, here is the code using the new cmxtmixlogit command:

use "http://fmwww.bc.edu/repec/bocode/t/traindata.dta", clear

set seed 12345
gen rnd = runiform()
bysort pid gid (rnd): gen alt = _n
cmset pid gid alt
cmxtmixlogit y price, r(contract local i.wknown tod seasonal) ///
             intpoints(100) intburn(15) intm(halton) nocons
margins r.wknown, outcome(1) alternative(1)

[Marc Ruggiano]

Hi Emmanuel Teitelbaum, did you ever get any additional information regarding your question #3? I have several questions about using the new CM commands and margins in choice experiments with unlabeled alternatives. I don't see a forum or board on this. Do you know of a Stata expert that could address some of the questions that seem to be going unanswered in various threads, possibly including this one? Thanks, Marc

[Hong Il Yoo]

Marc Ruggiano and Emmanuel Teitelbaum: As Emmanuel had already guessed, this is "simply unnecessary." By construction choice probabilities add up to 1 across all alternatives and changes in probabilities such as marginal effects add up to zero across all alternatives. The average of marginal effects across all alternatives is therefore zero.

[Emmanuel Teitelbaum]

Marc Ruggiano: I had an extensive back and forth with Joerg Luedicke, Senior Social Scientist and Statistician through [email protected]. I was looking for something like an Average Marginal Component Effect (AMCE) that political scientists (e.g. Hainmuller, Hopkins and Yamamoto) use in their standard OLS analysis of conjoint data. Joerge told me that it is not possible to get an "average marginal effect" for these kinds of models because it would not be "uniquely identifiable in a technical sense." In the end I went in a completely different direction, learning R and using hierarchical bayes for my analysis.

[Marc Ruggiano]

Thank you Hong Il Yoo and Emmanuel Teitelbaum. I am new to choice models and Stata's CM commands, and appreciate your guidance. Given your comments, here are a couple of questions.

After running the code from Manny's post above:

use "http://fmwww.bc.edu/repec/bocode/t/traindata.dta", clear
set seed 12345
gen rnd = runiform()
bysort pid gid (rnd): gen alt = _n
cmset pid gid alt
cmxtmixlogit y price, r(contract local i.wknown tod seasonal) ///
                     intpoints(100) intburn(15) intm(halton) nocons
margins r.wknown, outcome(1) alternative(1)

I run the following.

. *Example 1
. margins

Predictive margins        Number    of    obs     =    4,780
Model VCE    : OIM

Expression   : Pr(alt), predict()

                    
Delta-method
Margin   Std. Err.    z    P>z        [95% Conf.    Interval]
                    
_outcome 
1     .2398743   .0007613    315.10    0.000        .2383823    .2413664
2     .2533285   .0005554    456.13    0.000        .2522399    .254417
3     .2512891    .000736    341.43    0.000        .2498466    .2527316
4      .255508   .0005877    434.77    0.000        .2543562    .2566599

And then the following:

. *Example 2
. margins r.wknown

Contrasts of predictive margins                 Number of obs     =      4,780
Model VCE    : OIM

Expression   : Pr(alt), predict()

-------------------------------------------------------
                    |         df        chi2     P>chi2
--------------------+----------------------------------
wknown@_outcome#alt |
      (1 vs 0) 1 1  |          1      116.35     0.0000
      (1 vs 0) 1 2  |          1      116.86     0.0000
      (1 vs 0) 1 3  |          1      115.13     0.0000
      (1 vs 0) 1 4  |          1      116.49     0.0000
      (1 vs 0) 2 1  |          1      115.77     0.0000
      (1 vs 0) 2 2  |          1      115.32     0.0000
      (1 vs 0) 2 3  |          1      113.85     0.0000
      (1 vs 0) 2 4  |          1      113.41     0.0000
      (1 vs 0) 3 1  |          1      115.45     0.0000
      (1 vs 0) 3 2  |          1      114.19     0.0000
      (1 vs 0) 3 3  |          1      114.69     0.0000
      (1 vs 0) 3 4  |          1      114.32     0.0000
      (1 vs 0) 4 1  |          1      117.49     0.0000
      (1 vs 0) 4 2  |          1      114.67     0.0000
      (1 vs 0) 4 3  |          1      114.86     0.0000
      (1 vs 0) 4 4  |          1      114.84     0.0000
             Joint  |         11      281.95     0.0000
-------------------------------------------------------

---------------------------------------------------------------------
                    |            Delta-method
                    |   Contrast   Std. Err.     [95% Conf. Interval]
--------------------+------------------------------------------------
wknown@_outcome#alt |
      (1 vs 0) 1 1  |   .1346791   .0124861      .1102068    .1591513
      (1 vs 0) 1 2  |  -.0450467    .004167     -.0532138   -.0368795
      (1 vs 0) 1 3  |   -.044357    .004134     -.0524596   -.0362544
      (1 vs 0) 1 4  |  -.0450941    .004178     -.0532828   -.0369053
      (1 vs 0) 2 1  |  -.0452978     .00421     -.0535493   -.0370463
      (1 vs 0) 2 2  |   .1374194   .0127967      .1123383    .1625004
      (1 vs 0) 2 3  |  -.0465884   .0043663     -.0551462   -.0380305
      (1 vs 0) 2 4  |  -.0462908   .0043467     -.0548102   -.0377714
      (1 vs 0) 3 1  |  -.0444699   .0041388     -.0525819    -.036358
      (1 vs 0) 3 2  |  -.0462961   .0043323     -.0547873   -.0378049
      (1 vs 0) 3 3  |   .1369797   .0127909        .11191    .1620494
      (1 vs 0) 3 4  |  -.0459722   .0042996     -.0543994    -.037545
      (1 vs 0) 4 1  |  -.0449113   .0041435     -.0530323   -.0367903
      (1 vs 0) 4 2  |  -.0460766   .0043029     -.0545101   -.0376431
      (1 vs 0) 4 3  |  -.0460343   .0042954     -.0544532   -.0376155
      (1 vs 0) 4 4  |   .1373571   .0128174      .1122354    .1624787
---------------------------------------------------------------------

Couple of questions:

In an unlabeled choice experiment such as this one:

1. What explains the difference in choice probability among the four highlighted rows in Example 1?
2. How do I interpret the four slightly different choice probabilities when there is no meaning to the number label applied to each alt, other than order of presentation? I don't believe these are meaningful in an unlabeled experimen, other than perhaps as indicators of order bias or potential design issues.
3. What explains the difference in the contrasts for wknown in the four highlighted rows in Example 2?
4. How do I interpret the four slightly different contrast values when, again, there is no meaning to the number label applied to each alt, other than order of presentation? I believe these are meaningful values, and I am not sure whether to interpret these as a range (i.e., being well know increases probability of choice by 13.47-13.74%) or whether and how to select one of them over the other three (as indicated by the Stata rep.)

Thank you for your thoughts on this. And if you prefer I start a new thread for this, please let me know.

Marc

[Hong Il Yoo]

Marc Ruggiano: -cmxtmixlogit- was introduced after I had stopped working with mixed logit models, so I don't understand the syntax diagram and output. Having said, is there any reason why you're keen to compute the marginal effects? In all unlabelled choice experiment studies that I have read and authored, the focus is on the WTP estimates rather than the marginal effects. If you'd like to carry out a marginal effect analysis in a meaningful manner, I think you'll first have to think about a policy relevant choice scenario and deliberately select combinations of attribute-levels that describe alternatives in that scenario. Then you'll have to compute the marginal effects for that choice scenario, rather than compute the average marginal effects across actual choice scenarios used in your experiment.

[Emmanuel Teitelbaum]

Marc Ruggiano : I agree with Hong Il Yoo that you cannot get an average marginal effect using mixed logit. Instead you have to base your post-estimation analysis on relevant scenarios. If you really want average marginal effects, then the easiest thing is to cite Hainmuller et. al. and just run OLS.

There is a nice article by Lancsar et. al. et. al. that goes over how to analyze discrete choice experiments with various methods including mixed logit models. They use the old mixlogit command in Stata and show how to do marginal effects for various scenarios in a study of patient drug choice. They provide replication data online that you can use to try and replicate their analysis using the new cmxtmixlogit command with the margins post-estimation command. That might give you a nice idea of what is possible with it.

I would say that mixed logit models and hierarchical bayes are cool if you want to test specific real-world scenarios (e.g. what happens if I change the product color from red to blue)? But they are not well-suited for estimating average effects.
----------------
Hainmueller, J., Hopkins, D., & Yamamoto, T. (2014). Causal Inference in Conjoint Analysis: Understanding Multidimensional Choices via Stated Preference Experiments. Political Analysis, 22(1), 1-30. doi:10.1093/pan/mpt024


Lancsar, E., Fiebig, D.G. & Hole, A.R. Discrete Choice Experiments: A Guide to Model Specification, Estimation and Software. PharmacoEconomics 35, 697–716 (2017). https://doi.org/10.1007/s40273-017-0506-4

[Marc Ruggiano]

Hong Il Yoo Emmanuel Teitelbaum I am still working through this topic, but am making progress. Thank you both for you guidance - it has been helpful!

[Francois Dessart]

Dear all,
I am using the new choice model mixed logit command in Stata to analyse data from a DCE. More specifically, the cmxtmixlogit command.
Would you have any guidance as to how to estimate willingness to pay (WTP) using this type of model? There is nothing mentioned about WTP in the Stata cmxtmixlogit guide.
I would be grateful if you could direct me to any material that could help me.