Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
76,576 characters · 14 sections · 25 citation commands
Get me out of this hole: a profile likelihood approach to identifying and avoiding inferior local optima in choice models
The field of choice modelling has seen substantial growth over the last four decades, with the models now used widely across disciplines as diverse as transport, health, marketing and environmental research. Along with the growth in the number and diversity of applications, there has been a shift towards ever more complex model types, especially in relation to capturing complex correlation patterns and accommodating preference heterogeneity. Examples include generalized extreme value models (with nested logit as an important special case), continuous and finite mixture models (exemplified through mixed logit and latent class), hybrid choice models, and extensions such as multiple discrete continuous models and models using the Berry-Levinson-Pakes (BLP) approach (a.k.a. random coefficient logit demand systems).
Increased complexity of models is almost inevitably accompanied by a variety of practical challenges when computing parameter estimates for models. The context considered here is maximum likelihood estimation (MLE), which requires solving a non-linear optimization problem, but our results and discussion are readily applicable to any optimization-based estimators (e.g., other extremum estimators such as non-linear least squares or generalized method of moments). Many analysts are often unaware of (or underappreciate) the nuances and implications of performing this type of estimation, uncritically accepting at face value results produced by whatever estimation software they are using. Clearly understanding when a search algorithm has successfully converged to a valid solution, and/or whether there might be identification issues, is critically important. The ability of analysts to properly address all such concerns is paramount, especially when the results are to be used for decision making or policy analysis.
These issues arise for a wide range of complex econometric and statistical models, and the concerns expressed above are not necessarily new, a notable example being 10.1257/000282803322157133. They systematically enumerate these in detail, provide specific examples of analyst shortcomings, opine on the underlying causes, recommend a procedure for analysts to follow (which we discuss below), and demonstrate it via a case study.
One especially challenging situation is when there is the possibility of multiple local optima, the main topic of this paper. Returning to choice models, a natural starting point for any discussion is the multinomial logit (MNL) model, which represents the simplest (and essentially prototypical) case. In contrast to other families of models where the simplest special case is linear regression, MLE for the MNL model requires solving a non-linear optimization problem, which in the early years was far from trivial. It is well known that for an MNL model with a linear-in-parameters utility function specification, i.e. $V_i=\beta'x_i$, where $\beta$ is a vector of parameters and $x_i$ a vector of attributes for alternative $i$, the log-likelihood function is globally concave, and if an MLE exists inside the parameter space\footnote{For any dataset, there is no guarantee that the MLE will exist, particularly in small samples. Assuming there is no theoretical over-specification due to rank deficiency in the explanatory variables (as in linear regression, e.g. a missing normalisation for a categorical variable), the observed data must also not include any patterns whereby some or all discrete choices can be explained deterministically by the explanatory variables (e.g. an interaction effect for a subgroup that never chooses a specific alternative). For such cases, one or more of the parameters can improve the log-likelihood by approaching plus-or-minus infinity, yielding non-unique maximizers on the boundary of the parameter space.}, it will be unique, that is, the log-likelihood will be strictly concave at the solution, and moreover, will be a unique global maximum.
The MLE can then be computed using any well-implemented optimization algorithm that is guaranteed to converge to a local optimum (which in this case is also the global optimum) from an arbitrary starting value. For example, so-called quasi-Newton methods that use line searches or trust regions fit this description. In numerical optimization, these algorithms are called globally convergent (which should not be confused with `finding the global optimum')\footnote{It should be noted that the Newton-Raphson iteration approach described in many econometrics texts does not satisfy this requirement, that is, it is not guaranteed to converge to a local optimum from an arbitrary starting point.}. Moreover, specific versions of these algorithms ensure fast local convergence behaviour in contrast to others, e.g., gradient descent, or the expectation-maximisation (EM) algorithm. Examples of locally fast methods include BFGS [Broyden-Fletcher-Goldfarb-Shanno] and BGW [Bunch-Gay-Welsch] - for details, see Bunch2024.
The theoretical guarantee of a globally concave log-likelihood function does not apply once we move away from a linear-in-parameters MNL, whether considering MNL with utility specifications that are not linear in parameters, or more complex choice models. The question of whether the log-likelihood functions for these extended models remain globally concave becomes immediately murky, in the sense that it becomes very difficult to say one way or the other whether this is true by direct mathematical proof. There is thus a risk that there may be multiple local optima, which presents both theoretical and computational challenges.
Choice modellers routinely acknowledge the potential risk of multiple local optima, most notably in the context of latent class models, but most studies do little about it, often for example still relying on a single set of parameter starting values. The lack of published papers highlighting the pervasiveness of local optima with many advanced model specifications has likely added to the scepticism that some analysts express in relation to the risk of convergence to inferior solutions, and the faith many put in their estimation results.
Unfortunately, the econometrics/statistical literature offers no easy answers. The theory of maximum likelihood focuses on specifying conditions under which the estimator is consistent and asymptotically normally distributed, with the additional feature of being efficient (minimum variance). Econometrics texts typically consider MLEs of two types for purposes of proving their properties. One type (Type 1) is defined in terms of a global maximum, and the other (Type 2) in terms of local maxima that satisfy first-order conditions (a.k.a. the likelihood equations) on the interior of the parameter space. The first is problematic in the absence of global concavity, because reliably identifying a global maximum is essentially impossible to guarantee. Results for local maxima are useful for proving asymptotic normality (which are not available for Type 1) and guarantee that one of the local maxima is consistent. However, this is of limited usefulness since it provides no means of identifying a consistent maximum -- see, e.g., 222.
In fact alma990016651030403126 succinctly describes the grim situation:
A strategy that has evolved in the literature (particularly for mixture models) consistent with this description, given the availability of globally convergent search algorithms, is to randomly generate multiple starting values (in some cases, a relatively large number), find the local maximizer associated with each starting value, and choose the maximizer with the largest log-likelihood. Many starting values may end up converging to the same local maximizer, and those that do are said to be in the same `basin of attraction' for that local maximizer.
Although the approach of using randomly generated starting values is widely used, no known approach can guarantee that the global maximum will be found, and it is therefore difficult to tell how effective the approach might be for any particular problem. Alternative approaches employ methods from the global optimization literature, including heuristic approaches such as simulated annealing, tabu search, genetic algorithms, etc. These algorithms attempt to ensure that, for any current iterate in the search, there is some positive probability of `escaping' from whatever `neighbourhood' the iterate is in (that is, the basin of attraction to a specific local optimizer) to improve the chances of finding the global optimizer. Such approaches typically cannot guarantee finding the solution in finite time. A brief introduction to the heuristic approaches mentioned above (as well as others) is included in Bierlaire2010, who report that these methods have rarely been applied to econometric modelling in general, and discrete choice models in particular. In this spirit, they develop a heuristic method that attempts to “prematurely interrupt the local search if the iterates are converging to a local [optimum] that has already been visited, or if they are reaching an area where no significant improvement can be expected.” Perhaps due to its complexity, and despite implementation in Apollo hess_palma_apollo, the algorithm appears to have been little used in practice.
An alternative optimization approach (again, often used for mixture models), is the expectation-maximisation (EM) algorithm Train2009. The feature of EM that allows it to work in this context is that it is globally convergent, because it follows a path that ensures improvement in the log-likelihood at each step. However, its local convergence rate is extremely slow, and it can sometimes be difficult to tell when it has converged. One issue with EM is that some analysts are confused about its properties, thinking that it improves their chances of finding a global optimum from a single starting point.
As mentioned earlier, despite being routinely acknowledged as a problem, there is a lack of published results on the prevalence of local optima. A potential contributing factor is that for very complex models, there are multiple computational challenges that occur simultaneously, most notably identification and conditioning issues, that can introduce ambiguity regarding the nature of the results. A factor that further complicates the situation is that many of the recently developed complex choice models use simulation approaches to approximate integrals, as in mixed logit and hybrid choice models. The implications of this have not been well studied, but models like these can easily lead to cases where estimation software returns `solutions' that are highly ill-conditioned and even unidentified. The fact that, e.g., changing the number of random draws can yield different solutions further complicates the situation when it comes to investigating issues related to multiple optima.
For all of the above, there is a critical need for a highly rigorous characterization of the results produced by estimation software, a point emphasized by 10.1257/000282803322157133. They suggest a four-step process that we paraphrase here:
The first two steps require some familiarity with the technical details of non-linear optimization -- see Bunch2024 for details. The third step is probably more familiar to analysts, but it is important to use software that specifically performs it. Quasi-Newton methods appropriate for complex models use Hessian approximation approaches to avoid computing the full Hessian (which can be time-consuming), and thus it is not available at the end of the search. It is important to understand that any Hessian approximation that might be available at the conclusion of a search is insufficient for this purpose, so this step must be implemented separately.
The fourth step above is directly quoted from 10.1257/000282803322157133, and is a specific application of profile likelihood concepts that have long been employed in the literature, and which form the basis of the contribution of the present paper. As 10.1257/000282803322157133 note, this type of investigation was recommended to address issues arising from complex models in time series analysis, including the possibility of `multimodality.' In the present paper, we specifically put forward the use of profile likelihood to systematically examine the behaviour of the log-likelihood in the neighbourhood around an identified local maximizer, and to identify better optima that might exist in alternative `neighbourhoods' adjacent to the current one. This thus provides modellers with a method to a) better diagnose the properties of a solution they have obtained, and b) assist with finding a better optimum. What we propose is a more structured approach than relying on multiple (random) starting values and hoping that they converge to different solutions, and can be readily implemented using the same estimation tools.
Before proceeding, a caveat is in order. This paper is about introducing an approach that helps an analyst determine if the local optima that are most readily obtained using standard estimation are in fact inferior. The paper is not concerned with the details of what might cause local optima, or how an analyst may reduce the risk of local optima by changing model specification, implementation, or the estimation approach. Our starting point is that a user makes a decision on model structure, utility specification, starting values for parameters, and implementation details (including the number of draws for simulation-based estimation). Subject to these decisions, they obtain a solution. Our approach allows the analyst to then test, conditional on these decisions, whether there are better local optima in adjacent neighbourhoods. Of course, an analyst can then also use our approach to investigate whether these risks are lower with different model specification and implementation decisions, for example.
Another important point to raise is the distinction between local optima that are identical solutions and those that are not. For example, with latent class models, there exists the possibility of class switching. With $K$ classes, there are $K!$ local optima that obtain the exact same log-likelihood, each with a different ordering of the classes. A similar notion arises when using simulation-based estimation for models with continuous random heterogeneity, where, if the draws used for $\xi_k$ in $\beta_k=\mu_k+\sigma_k\xi_k$ are symmetrical around zero\footnote{This is not generally the case unless a user relies on antithetic draws.}, $2^K$ local optima with identical log-likelihoods exist for all combinations of positive and negative values for $\sigma_k$. The focus of the present paper is not on such local optima that are identical but on local optima that lead to different model results (e.g. different estimates of willingness-to-pay). These different local optima will also imply differences in log-likelihood, albeit that these differences may be (much) smaller than expected, a point we return to later in the paper.
The remainder of this paper is organised as follows. Section (ref) introduces the profile likelihood approach for identifying and avoiding inferior local optima. Section (ref) presents a case study on a well-used dataset, highlighting the applicability of the approach for both latent class and mixed logit. Finally, Section (ref) offers some conclusions and ideas for future research.
In this section, we review the key concepts and definitions that are relevant for the remainder of the paper. Our focus is on maximum likelihood, but, as noted in the introduction, similar approaches can be applied to other extremum estimators.
Assume that we have a sample of $N$ people, with $Y_n$ giving the sequence of observed choices for person $n$. The analyst specifies a model $L(Y|\beta)$ giving the probability of observing $Y$ conditional on $\beta$ (and some explanatory variables, which for convenience are suppressed), where $\beta$ (of length $K$) is a vector of model parameters with $\beta=\left[\beta_1,\hdots,\beta_K\right]$. The log-likelihood function for this model is then given by:
where $L_{n}\left(Y_n\mid\beta\right)$ is the likelihood of the observed sequence of choices for person $n$, and the subscript $N$ in $LL_N(\beta)$ denotes the sample size. As noted in the introduction, there are two estimation concepts in the literature: (1) the value $\hat{\beta}$ that globally maximizes $LL_N(\beta)$ over a parameter space $B$, or (2) a value $\hat{\beta}$ that satisfies the likelihood equations (first-order conditions) and is a strong local maximum (on the interior of the parameter space). The standard estimation approaches discussed in the introduction employ the second concept.
MLE analysis identifies a local maximum and typically performs inference using so-called Wald statistics based on a sample estimate of an asymptotic variance-covariance matrix. Under the conditions that (1) the model is correctly specified and (2) there is a unique, true parameter vector $\beta^*$, there are multiple options available for how to compute an estimate of the variance-covariance matrix. Because the MLE is an extremum estimator, the limit distribution of $\sqrt{N}(\hat{\beta}-\beta^*)$ has variance matrix $\Omega=A_0^{-1}B_0A_0^{-1}$ (the so-called sandwich form), and so $\hat{\beta}$ has asymptotic variance matrix $\Omega=N^{-1}A_0^{-1}B_0A_0^{-1}$ where\footnote{An accessible reference for this material is alma990024527500403126. A caveat to add is that these derivations are for the case of one choice per individual.}:
For maximum likelihood, the above conditions mean that $B_0$ is the Fisher information matrix ($I$), and that the Information Matrix (IM) equality $-A_0 = B_0$ holds. Moreover, the estimates are minimum variance.
Using the earlier expression for $\Omega$, the IM equality implies that $-A_0^{-1}$ and $B_0^{-1}$ are asymptotically equivalent expressions for $\Omega$. There are multiple options for consistent sample estimates of $A_0$ and $B_0$ evaluated at $\hat{\beta}$, but most estimation codes use the Hessian, or possibly the sum of Outer Products (OP) of partial derivatives, respectively:
Estimation software can vary in which $\hat{\Omega}$(s) are available:
A standard practice is to perform inference under the assumption that $\beta$ is asymptotically normal, with mean $\hat{\beta}$ and variance $\hat{\Omega}$, where many possible tests are available using the Delta method Daly2012a. The simplest case is to perform a test of the null hypothesis $\beta_k=0$ using an estimated standard error $\hat{\sigma}_k=\sqrt{\hat{\Omega}_{kk}}$ to compute the asymptotic $t$-statistic. Similarly, a $1-\alpha$ two-sided confidence interval for $\beta_k$ would be estimated as $\hat{\beta}\text{ } \pm\text{ }z_{{\alpha}/{2}}\hat{\sigma}_k$, with $z$ giving the appropriate critical value from a $N\left(0,1\right)$ distribution. These approaches are associated with Wald (as in `Wald statistic' or `Wald confidence interval').
To the degree that the MLE has achieved asymptotic behaviour, a quadratic model constructed using Hessian information at the solution will yield a good approximation to the log-likelihood function over a large enough area in the parameter space to yield correct inferences. However, for complex models using non-linear estimation, there is a risk that treating the MLE as though it has attained asymptotic behaviour could be a poor assumption, as also discussed by 406. For this reason, the idea of examining the shape and behaviour of the likelihood function itself has been extensively explored in the literature, and this leads us directly to the profile likelihood approach that we discuss in the next section.
The idea of `profiling the likelihood' to `assess the adequacy of the quadratic approximation' is step 4 of 10.1257/000282803322157133's recommended procedure, and they provide details on how to do this using graphical methods. They also remark that this type of investigation could reveal a lack of monotonicity in the log-likelihood, which would be a strong indicator that multiple optima exist.
To understand the notion of a profile likelihood approach, we revisit likelihood ratio test procedures. Consider the null hypothesis $H_0: \beta_k = \tilde{\beta}_k$. The likelihood ratio test statistic is
i.e., $LRT$ is distributed $\chi^{2}$ with $1$ degrees of freedom (due to the single constraint on $\beta_k$). Since we are going to limit our attention to the behaviour of the log-likelihood when varying one parameter at a time, it is useful to define the profile log-likelihood by
This can be used to define a likelihood-based confidence interval
that could be compared to the Wald-based CI described in the prior section.
However, note that implementing this requires some computational effort. In our case, we are primarily interested in examining the behaviour of the log-likelihood when $\tilde{\beta}_k$ is varied over some interval around $\hat{\beta}_k$. A useful benchmark is the estimated standard error ($\hat{\sigma}_k$). One obvious comparison involves the upper and lower bounds of the Wald 95% confidence interval given by $\tilde{\beta}_k=\hat{\beta}_{k}\pm 1.96\sigma_{k}$. For either the upper or lower bound, the expected drop in the value of the log-likelihood would be half of $\chi^2_{1,\alpha=0.05}$ ($\frac{3.84}{2}$), if $\hat{\beta}$ were behaving asymptotically.
Of course, the main point is that the small-sample MLE properties for complex models will often be different, and the log-likelihood function could deviate substantially from quadratic behaviour. Recognizing this, 10.1257/000282803322157133 suggest a starting point of $\bar{\beta}_k=\hat{\beta}_{k}\pm 4\sigma_{k}$, which can be readily modified based on initial findings.
The estimation of our base model gives us $LL_{base}$, with MLE $\hat{\beta}_{base}=\left[\hat{\beta}_{1,base},\hdots,\hat{\beta}_{K,base}\right]$ and associated standard errors $\hat{\sigma}_{base}=\left[\hat{\sigma}_{1,base},\hdots,\hat{\sigma}_{K,base}\right]$. The profile likelihood approach now consists of three steps.
\paragraph{Step A1: Define candidate values\\}
In this step, we produce $M$ values for $\beta_k$ going either side of the MLE $\hat{\beta}_k$, $\forall k$. Specifically, we obtain candidate value $m$ for parameter $k$ as:
where $\gamma$ is a vector of $M$ evenly spaced points between $\gamma_a$ and $\gamma_b$. The choice of $M$ and the width $\gamma_b-\gamma_a$ are analyst decisions that we return to in our case study below, where we can already say that this involves a trade-off between robustness and computational cost. A decision is also needed by the analyst on whether to use classical or robust standard errors in Equation (ref).
\paragraph{Step A2: Estimate constrained models\\}
In this second step, $M$ constrained models are estimated for each of the $K$ parameters. Specifically, in model $m$ for parameter $k$, we fix $\beta_k=\tilde{\beta}_{k,m}$ and freely estimate the remaining $K-1$ parameters, starting from $\hat{\beta}_l,\,\forall l\neq k$. We denote the log-likelihood for the resulting constrained model as $\widetilde{LL}_{m_k}$, with $\hat{\beta}_{m_k}=\left[\hat{\beta}_{1,m_k},\hdots,\tilde{\beta}_{k,m},\hdots,\hat{\beta}_{K,m_k}\right]$, i.e. a vector of freely estimated parameters except for the $k^{th}$ element which is constrained to $\tilde{\beta}_{k,m}$.
\paragraph{Step A3: Evaluate results\\}
The final step compares the model fit for the base model, i.e. $LL_{base}$, and the fits for the additional $K\cdot M$ models estimated in step A2, i.e. $\widetilde{LL}_{m_k},\,\forall k,m$.
First, if $$LL_{base}>\widetilde{LL}_{m_k},\,\forall k,m,$$ there is some “guarantee” that estimation has converged to “good” local optimum, especially with large $M$ and wide $\gamma_b-\gamma_a$.
Second, if $$LL_{base}-\widetilde{LL}_{m_k}\sim 1.92,\,\text{for}\mid\gamma_{m_k}\mid\sim 1.96,$$ then the asymptotic normality properties hold, at least in an area close to the MLE.
Otherwise, if we find cases where the likelihood of candidate points exceed the base likelihood, then we know our original estimates formed a local optimum only and that we have found an improved value. The following section describes an iterative procedure to find better local optima.
Section (ref) describes the profile likelihood approach for identifying superior local optima in the area around the original MLE solution. If step A3 of that approach identifies such solutions, then these can be used as a starting point for further improving the model.
We propose a simple iterative approach for this, based on a number of steps.
\paragraph{Step B1: Select new starting points\\}
Any solutions where $\widetilde{LL}_{m_k}>LL_{base}$ become new starting points, say $\beta_{round_2,s}$ for solution $s$ out of $S$ better local optima.
\paragraph{Step B2: Unconstrained estimation\\}
Step A2 in Section (ref) imposes the constraint that $\beta_k=\tilde{\beta}_{k,m}$, meaning that $\widetilde{LL}_{m_k}$ is the result of a constrained estimation. By using the resulting estimates as starting values but removing the constraint on $\beta_k$, a further improvement is possible. This step is now repeated for all $S$ starting points identified in step B1 above, giving a new set of solutions, with model $s$ having a log-likelihood of $LL_{round_2,s}$ and MLE $\hat{\beta}_{round_2,s}$.
\paragraph{Step B3: Filtering of solutions\\}
It is likely that the unconstrained estimation in step B2 has converged to the same solution in some of the $S$ runs, and in this step, the resulting solutions are filtered to produce a unique set of $S^*$ new solutions.
\paragraph{Step B4: Iterate\\}
After completing step B3, an analyst will have one or more solutions that are better than $LL_{base}$. However, the properties of these are again unknown in relation to the existence of better local optima in the near vicinity, and we now return to step A1 with the new solutions as starting point, using e.g. $\hat{\beta}_{round_2,s^*}$ as a starting value.
For maximum robustness, this process should be repeated with all $S^*$ solutions found in step B3. For more computationally demanding models such as Mixed Logit, an analyst may take the pragmatic decision to use the best solution from step B3 as the new starting point for A1.
The data used for our case study comes from a stated choice (SC) survey of public transport route choice conducted in Switzerland 814. A total of $388$ people were faced with $9$ choices each between two public transport routes, both using train (leading to $3,492$ observations in the data). The two alternatives are described on the basis of travel time (tt), travel cost (tc), headway (hw, time between successive trains) and the number of interchanges (ch).
The data used in our paper have been subjected to extensive analysis by others, and it should be noted that there are additional variables (e.g. income, trip purpose) that can be used to capture observable heterogeneity, as well as other interaction effects. However, for this case study we limit ourselves to the three model specifications that follow. The first (MNL) is a natural starting point for illustrative and comparison purposes. The other two represent alternative specifications for capturing the effect of unobserved heterogeneity.
For our profile likelihood implementation, we worked with $\gamma_a=-4$, $\gamma_b=4$, and $M=51$, ensuring that the middle point equates to the base MLE solution. We also used robust standard errors in Equation (ref), thus further widening the range.
All models were estimated using Apollo hess_palma_apollo using the BGW optimiser bgw.
Table (ref) reports the estimation results for our MNL model. We see that all four coefficients are negative, in line with expectations, and with robust $t$-ratios that allow us to reject the null hypothesis of no effect at any reasonable level of significance for all four LOS variables. We also report the implied valuations of travel time (VTT), headway (VHC) and interchanges (VCH), where the VHW is lower than VTT, as would be expected, while the monetary value of avoiding an interchange in this case is around the same as that of reducing in vehicle travel time by $19$ minutes. Asymptotic $t$-ratios for these valuations are obtained via the Delta method Daly2012a.
The maximum Hessian eigenvalue is negative (indicating convergence to a maximum) and far from zero. We also report the reciprocal of the conditional number for the Hessian\footnote{The condition number of a matrix ($A$) is an estimate of how ill-conditioned $A$ is with regard to numerical computations such as matrix inversion or solving a linear system $Ax = b$. In this context, the term “ill-conditioned” means that a small change in $A$ (or in $b$, in the case of solving a linear system) will produce a very large change in the outcome of the numerical computation. That is, the results are completely unreliable. The most extreme case of ill-conditioning is a singular matrix, which has a condition number of infinity. This is why, when we want to estimate conditioning, we compute an estimate of the reciprocal of the condition number, which for a singular matrix would be zero. The smaller the reciprocal condition number, the more ill-conditioned the matrix is. Conditioning of the Hessian matrix is closely linked to whether the solution is identified or not. There is no hard and fast rule, but one thing to look for is whether the reciprocal condition number is less than the square root of machine $\epsilon$, i.e. the smallest floating-point number such that $1+\epsilon$ is distinguishable from $1$.}, which is substantially larger than the square root of machine precision ($\sqrt{\epsilon}=1.490116e-08$). Both of these findings suggest that the solution is a numerically stable and well-defined optimum.
\afterpage{
}
\paragraph{Profile likelihood: steps A$1$ to A$3$\\}
For a linear-in-parameters MNL model, we of course know that the log-likelihood function is globally concave. To illustrate the impact of this on a profile likelihood analysis, we report the findings of step A$3$ from Section (ref) in Figure (ref). This shows that the LL drops monotonically on either side of the MLE for all parameters. In addition, we see that the asymptotic normality properties hold, where a move away from the MLE by $\pm 1.96\sigma_k$ leads to a drop in LL by 1.92 units, in line with the discussions in Section (ref).
We next move to our latent class (LC) model, with the base estimation results reported in Table (ref). We used small shifts around MNL results as starting values for the parameters in the two classes. We see that the coefficients for the four LOS variables are negative in both classes, where the robust $t$-ratios further indicate that all $8$ estimates are different from zero above the $99\%$ level of confidence using a one-sided test. Class $1$ is characterised by higher monetary valuations for travel time, headways and interchanges. This class is also larger than class $2$. In comparison with MNL, the (mean) population VTT is $33\%$ higher, the VHW is $6\%$ lower, and the VCH is $11\%$ higher. The class structure of course then also leads to heterogeneity in the population valuations, as shown in Table (ref).
The LL and Bayesian Information Criterion (BIC) for LC are clearly higher than for the MNL model in Table (ref), and a Ben-Akiva & Swait test benakiva_swait conclusively rejects the null hypothesis of no difference between the models ($p<10^{-38}$). The maximum Hessian eigenvalue is still negative (indicating convergence to a maximum) and although it and the reciprocal of the condition number are closer to zero than they were for MNL, an analyst would be justified in concluding that the estimation has converged to a proper local optimum.
\afterpage{
}
\paragraph{Profile likelihood: steps A$1$ to A$3$\\}
We now proceed with the profile likelihood analysis for our LC model. Steps A$1$ and A$2$ implied the setup and estimation of $450$ constrained models\footnote{We do not include in this count the $9$ models where $\gamma_m=0$, which exist due to setting $M=51$ with $\gamma_a$ and $\gamma_b$ symmetrical around $0$. These $9$ models are in fact at the MLE and do not require a new estimation.} involving shifts away from $\hat{\beta}_{base}$. Turning to step A$3$, across these $450$ models, we obtain an improvement in LL for $99$ models, compared to $LL_{base}$. As shown in Figure (ref), better local optima are identified for all parameters except for $\beta_{tc}^{\left(1\right)}$ and $\beta_{hw}^{\left(1\right)}$. For the remaining $7$ parameters, at least $10$ new local optima are discovered out of the $50$ constrained models for each parameter.
\paragraph{Profile likelihood: steps B$1$ to B$3$\\}
Given that step A$3$ identified solutions superior to the base estimation results, we next proceed with the iterative algorithm from Section (ref). Step B$1$ uses the estimates from the $99$ better solutions found in A$3$ to create starting values for the unconstrained estimation in step B$2$. In step B$3$, we then filter the results from these $99$ estimation runs, which in our case yields $4$ new unique solutions, labelled as $s_1^*$ to $s_4^*$.
\afterpage{
}
The estimation results for the new solutions are shown in Table (ref). All four models obtain estimates for the impact of the four LOS variables that are negative in both classes, and clearly different from zero at common levels of significance testing. Models $s_1^*$ and $s_2^*$ obtain solutions that are very similar in fit and also in implied monetary valuations. Models $s_3^*$ and $s_4^*$ are also close to each other in fit, but with larger differences in implied valuations. Looking at the convergence information, we see that for all four models, the findings indicate a proper optimum.
\paragraph{Profile likelihood: step B$4$\\}
The profile likelihood approach identified four unique new optima when we applied it to our base solution. There is at this stage no guarantee yet that there are not other better local optima near these solutions, and we thus apply step B$4$, which uses the results from step B$3$ as starting values for a further round of profile likelihood analysis, starting with A$1$.
Figure (ref) shows the results for the profile likelihood analysis starting from solution $s_1^*$. While there is evidence of non-monotonic decreases in the LL away from the optimum, any local optima identified are not better than the original solution $s_1^*$. Crucially, compared to Figure (ref), the shape of the log-likelihood function around the $s_1^*$ solution more closely exhibits asymptotic normality properties.
Turning to solution $s_2^*$ in Figure (ref), we identify $20$ new solutions better than $LL_{s_2^*}$. After applying steps B$1$ to B$3$ on these solutions, these reduce to a single solution, which is identical to $s_1^*$, implying that this process again leads to solution $s_1^*$.
For solution $s_3^*$ (cf. Figure (ref)), a similar story emerges. We identify $32$ new local optima in step A$3$, which step B$3$ reduces to two unique solutions that in this case are identical to $s_1^*$ and $s_2^*$. We already know that starting the iterative profile likelihood algorithm with solution $s_2^*$ leads to $s_1^*$, meaning that we again can conclude the process here.
Turning finally to $s_4^*$ (cf. Figure (ref)), we identify $88$ new local optima in step A$3$, which step $B3$ reduces to $2$ unique solutions, which are identical to $s_1^*$ and $s_3^*$. We can thus again conclude the process as we know that starting from $s_3^*$ will lead to $s_1^*$.
\afterpage{
}
\afterpage{
}
\paragraph{Differences across local optima\\}
The four different new solutions offer improvements over the base solution $LL_{base}$ by between $13.72$ and $25.73$. These differences in LL are under 2%, which contrasts with the much larger differences we see in the implied monetary valuations.
\afterpage{
}
In our case, the profile likelihood eventually leads to solution $s_1^*$ also when starting from $s_2^*$, $s_3^*$ or $s_4^*$, so our focus is on comparing the results for that model to the base model results. For VTT, we see a drop in the population mean by $10.07\%$, while the standard deviation drops by $21.13\%$. On the other hand, for VHW, we see increases by $30.94\%$ for the mean and $93.69\%$ for the standard deviation. Finally, for VCH, the mean drops by $10.93\%$, while the standard deviation drops by $15.89\%$.
For both estimates, the log-likelihood profiles were reasonably well-behaved for all parameters over the nominal $95\%$ confidence regions in Figures (ref) and (ref), respectively. A direct comparison of estimates along with their asymptotic $95\%$ confidence intervals are provided in Figure (ref). With the exception of the standard deviation for VHW, the values for the $s_1^*$ solution fall within the $95\%$ confidence interval for the original solution. Notwithstanding this, the substantial differences in the actual values would lead to very different policy implications. This illustrates the unfortunate situation described in the earlier quote from alma990016651030403126, but the analyst would appear to have little choice but to adopt solution $s_1^*$.
Additionally, it should be highlighted again that $s_1^*$ exhibits properties that are more in line with asymptotic normality (contrasting Figure (ref) and Figure (ref)), which suggests that the profile likelihood approach has led not just to a slightly better log-likelihood, but also a better behaved solution.
One final additional point worth raising is that the most striking difference between the $\hat{\beta}_{base}$ and $s_1^*$ solution is the much more negative cost sensitivity picked up for one class ($\beta_{tc,2}$) in the $s_1^*$ solution. This solution is much further away from 0 and outside the space that most reasonable starting value searches would cover. At the same time, the presence of several other strong (but inferior) local optima between $0$ and this solution means that estimation processes starting closer to $0$ are likely to get stuck in one of these other local optima. It is only by spreading the search more widely as done by our profile likelihood approach that we end up in solution $s_1^*$.
We finally move to our Mixed Multinomial Logit (MMNL) model, with the base estimation results reported in the left hand side of Table (ref). As our $\beta$ coefficients follow a multivariate Lognormal distribution, the estimated parameters relate to the multivariate Normal distribution of the logarithm of $\beta$, where, in the case of $\beta_{tc}$, it is the logarithm of $-\beta_{tc}$, given the use of a negative Lognormal. We used the logarithms of the MNL estimates as starting values for the means of the logarithms of $\beta$, with $0$ values for the Cholesky matrix parameters, i.e. starting at the MNL solution, which is common practice.
We first look at the mean parameters. We see that $\mu_{\log(\beta_{\text{vtt}})}$, $\mu_{\log(-\beta_{\text{tc}})}$ and $\mu_{\log(\beta_{\text{vhw}})}$ are negative and clearly different from zero, indicating that the median of the value of time parameter ($\beta_{\text{vtt}}$) and the value of headway parameter ($\beta_{\text{vhw}}$) are smaller than $1$ (after applying the exponential), while the median for the cost (scale) coefficient $\beta_{\text{tc}}$ is less negative than $-1$ (after applying the negative exponential). On the other hand, $\mu_{\log(\beta_{\text{vch}})}$ is positive, meaning that for the value of interchange parameter ($\beta_{\text{vch}}$), the median is greater than $1$. The $c$ terms relate to the Cholesky decomposition. The parameters are all different from zero, including the terms relating to covariances, justifying the use of a multivariate distribution.
Compared to the best of the LC models (cf. solution $s_1^*$ in Table (ref)), we see that mean valuations for travel time, headway and interchanges increase by $6.36\%$, $11.87\%$ and $29.79\%$, respectively. For the heterogeneity (in terms of standard deviations), the increase is much larger, at $156.78\%$, $412.32\%$ and $410.56\%$ - this finding is not surprising when contrasting a continuous Lognormal distribution with a two-value discrete distribution. We also see strong positive correlations between the three different valuations, which are in line with expectation. Additionally, we see small positive correlations between $\beta_{tc}$ and the three valuation terms.
The Bayesian Information Criterion (BIC) for MMNL is clearly higher than for the best LC model (cf. Table (ref)), and a Ben-Akiva & Swait test clearly rejects the null hypothesis of no difference between the models ($p<10^{-64}$). The maximum Hessian eigenvalue is more negative than for LC, and the reciprocal of the condition number is further from zero than for LC. An analyst would again be justified in concluding that the estimation has successfully converged to a proper optimum.
\paragraph{Profile likelihood: steps A$1$ to A$3$\\}
We now proceed with the profile likelihood analysis for our MMNL model. Steps A$1$ and A$2$ implied the setup and estimation of $700$ constrained models involving shifts away from $\hat{\beta}_{base}$. Turning to step A$3$, across these $700$ models, we obtain an improvement in LL for $210$ models, compared to $LL_{base}$. As shown in Figure (ref), better local optima are identified for all parameters. In this case, the profile likelihood approach also shows a notable lack of asymptotic behaviour, and a very unstable pattern overall, more so than for LC.
\paragraph{Profile likelihood: steps B$1$ to B$3$\\}
Given that step A$3$ identified solutions superior to the base estimation results, we next proceed with the iterative algorithm from Section (ref). Step B$1$ uses the estimates from the $210$ better solutions found in A$3$ to create starting values for the unconstrained estimation in step B$2$. In step B$3$, we then filter the results from these $210$ estimation runs, which in our case yields $76$ new unique solutions, labelled as $s_1^*$ to $s_{76}^*$.
We focus on the best of these $76$ solutions in terms of model fit, given by $s_{22}^*$, for which we show the results on the right hand side of Table (ref). Compared to the findings for LC, the differences in fit between the base model and the best profile likelihood solution are much smaller (1.77 units in LL). All parameter estimates retain the same sign as for the base model, and the fact that $\mu_{\log(-\beta_{\text{tc}})}$ is no longer different from zero simply implies that the median of $\beta_{tc}$ is not different from $-1$. Looking at the convergence information, the findings again indicate convergence to a proper optimum.
\paragraph{Profile likelihood: step B$4$\\}
The profile likelihood approach identified $76$ unique new optima when we applied it to our base solution. There is at this stage no guarantee yet that there are not other better local optima near these solutions, and we thus apply step B$4$. In contrast with the LC model, given the higher estimation costs, we use the version of the algorithm that uses only the best solution from step B$3$ as starting values (i.e. $s_{22}^*$) for a further round of profile likelihood analysis, starting with A$1$.
Figure (ref) shows the results for the profile likelihood analysis starting from solution $s_{22}^*$. While there is still some evidence of non-monotonic decreases in the LL away from the optimum, any local optima identified are not better than the original solution $s_{22}^*$. Crucially, compared to Figure (ref), the shape of the log-likelihood function around the $s_{22}^*$ solution much more closely exhibits asymptotic normality properties. This difference highlights the benefits of the proposed approach not just for finding better local optima, but also ones that have better asymptotic properties.
\afterpage{
}
\paragraph{Differences across local optima\\}
We finally contrast the results between the base solution and $s_{22}^*$ in terms of the implied monetary valuations. Alongside the values in Table (ref), we again show the confidence intervals for the population mean and standard deviations of the valuations in Figure (ref), and for the correlations between random terms in Figure (ref). While the differences for VTT and VCH are negligible, we see a drop in the mean VHW by $8.93\%$, while the standard deviation is reduced by $18.16\%$. The correlations reduce in magnitude except for $\rho\left(\beta_{vtt},\beta_{vhw}\right)$ and $\rho\left(\beta_{vtt},\beta_{vch}\right)$. While the new valuations remain within the confidence intervals of the original valuations, the differences for VHW would again have the potential for differences in policy implications. Additionally, it should be highlighted again that $s_{22}^*$ exhibits properties that are more in line with asymptotic normality (contrasting Figure (ref) and Figure (ref)). This could suggest that the profile likelihood approach has led to a better behaved solution, but definitely highlights the fact that the original estimates may not be as stable as they appear. An analyst may also find reassurance in the fact that the confidence intervals are narrower for all three mean valuations, for two out of the three standard deviations for valuations (with the exception of VHW), and for all correlations except $\rho\left(\beta_{tc},\beta_{vhw}\right)$.
\afterpage{
}
With the wide availability of more flexible estimation software packages, researchers have been increasingly using choice models that are much more complex than standard MNL models to capture a variety of phenomena, especially relating to preference heterogeneity. In many cases, it is not practical to determine whether the log-likelihood function might be globally concave, and there is a risk of multiple local optima. This has been recognized as a challenging problem across many areas of quantitative social science and in statistics, and there appears to be no clearly agreed upon approach to addressing it. This is particularly well-known for various classes of mixture models, and a variety of procedures have been developed that rely on randomly generated starting values and the EM algorithm and/or quasi-Newton methods. There has also been some application of concepts from the global optimization literature.
Our area of specific interest is choice modelling, where this issue has received limited attention even when researchers recognize the risk of multiple local optima. In this paper we have developed and demonstrated an approach based on the concept of profile likelihood that systematically explores the parameter space in a way that identifies potentially better solutions in `nearby neighbourhoods' of an original solution.
We demonstrate the potential benefit of our approach in applying both latent class and mixed logit specifications to a well-used stated choice dataset. In both cases, multiple local optima are shown to clearly exist after an initial search based on `reasonable' starting values achieves a proper local solution. In the case of latent class, we initially find $4$ new solutions, while in the case of mixed logit, the number is much higher, at $76$ new solutions. If an analyst was to spread the net wider, an even larger number may be found, highlighting the pervasiveness of the issue.
The profile likelihood procedure produces a sequence of improved solutions that may or may not have been found by using randomly generated starting values. In both cases, the process culminates in a final solution that has the largest log-likelihood from among a competing set of improved candidates, and that also does not lead to any additional improvements after examining profile likelihoods over a relatively wide range for each parameter. Moreover, the contours of the profile likelihoods show more `stability' and behaviour that is `more asymptotic' than those found in earlier iterations of the procedure. We found this observation to be encouraging, and potentially worth pursuing in future work.
In the case of the latent class models, we see that multiple solutions with quite similar log-likelihood values can have very different policy implications, where the best solution from a set of four has notably different willingness-to-pay estimates than the initial solution. Conversely, two of the four with the largest log-likelihoods produce WTP values that are similar, so concerns about `which one is the true consistent estimate' are in this case minimized. The analogous mixed logit results showed more similar results between the initial and final solutions, but with the caveat that the standard error estimates for both are relatively large. In both models, it is clear that the differences in parameter estimates for different local optima can be much larger than what would be expected given the differences in log-likelihood.
While analyst may be concerned about the computational cost of the profile likelihood approach, efficient implementation can help manage this - for example, any pre-processing of the model such as computing analytical derivatives does not need to be repeated for each new trial value for $\beta_k$. The higher computational cost should also be contrasted with what we posit to be increased robustness compared to simply using a set of random starting values - especially the findings for the latent class model suggest that the range that a typical analyst would try might not be wide enough. Furthermore, the rigorous nature of the profile likelihood approach also has the added advantage of providing analysts with direct insights into the shape of the log-likelihood function near their solution, and to what extent this adheres to the often taken for granted property of asymptotic normality.
In addition to demonstrating the potential advantages of this approach, an important underlying theme of this work is that researchers need to be much more rigorous and take more care when producing model estimates using non-linear optimization-based estimators, and review the key steps to consider in doing so. This theme is an underpinning to the implementation our method. That is, we took great care to use very reliable (and fast) estimation software with well-implemented stopping rules that provide diagnostic warnings if potential problems with a solution were detected.
Specifically, computing the likelihood profiles required repeated estimations to be performed over many points, where both speed and reliability are important. We monitored all solutions for negative diagnostic flags, and none were generated. However, for any complex choice model it is possible in some cases that a warning may not be raised, yet a more complete analysis of the full Hessian at a potential solution could still indicate problems. In any case, no solution should be adopted before performing such an analysis. In our case, we performed complete Hessian analyses for those solutions that were candidates to be improved optima, and found that they were valid.
A related issue is the problem of ill-conditioning, which also relates to an analysis of the full Hessian. Formally, any solution where the Hessian is positive definite is ‘identified’; however, as the largest (least negative) eigenvalue gets closer and closer to zero, there should be increased concern about whether the current specification is a prudent choice or can be improved upon. This is a worthy topic for investigation in its own right, and while it may nominally be outside the scope of this study, it is one more computation-related challenge that can also arise with the type of complex models that potentially give rise to multiple optima, and should not be ignored.
In closing, it should be noted again that the present paper was concerned with the identification of local optima and the avoidance of poor local optima. We did not investigate what causes local optima, whether different starting values help, or what the role of the estimation procedure is, including the implementation of simulation-based estimation (e.g. number and type of draws). These are important topics for future research, which are however all reliant on an approach that helps identify different local optima, which this paper has provided.
Stephane Hess acknowledges the support of the European Research Council through the advanced grant 101020940-SYNERGY.