A Latent Class Nested Logit Model for Rank-Ordered Data with Application to Cork Oak Reforestation

Abstract

We analyze stated ranking data collected from recreational visitors to the Alcornocales Natural Park (ANP) in Spain. The ANP is a large protected area which comprises mainly cork oak woodlands. The visitors ranked cork oak reforestation programs delivering different sets of environmental (reforestation technique, biodiversity, forest surface) and social (jobs and recreation sites created) outcomes. We specify a novel latent class nested logit model for rank-ordered data to estimate the distribution of willingness-to-pay for each outcome. Our modeling approach jointly exploits recent advances in discrete choice methods. The results suggest that prioritizing biodiversity would increase certainty over public support for a reforestation program. In addition, a substantial fraction of the visitor population are willing to pay more for the social outcomes than the environmental outcomes, whereas the existing reforestation subsidies are often justified by the environmental outcomes alone.

Appendices

Appendix 1: Summary of EM Algorithm

This appendix summarizes the expectation-maximization (EM) algorithm for estimating our latent class nested rank-ordered logit (LC-NROL) model. The use of the EM algorithm has been motivated by generic numerical difficulties associated with estimating finite mixture or latent class models via direct maximization of the sample log-likelihood function, not by any peculiar issue arising from estimating the LC-NROL model. While our own program was written in TSP International 5.1, it does not rely on any of TSP International’s specialized features. Programming the EM algorithm entails setting up a rudimentary loop over updating tasks (20) and (21) to be explained below, and can be implemented in any software package that allows the user to supply a self-written log-likelihood function e.g. Stata (Pacifico and Yoo 2013).

Unless specified otherwise, we follow the notations introduced in Sect. 3, and let $n=1,2,\ldots ,N$ index individuals and $c=1,2,\ldots ,C$ index classes. Our objective is to estimate three types of parameters: dissimilarity coefficient $\tau$, class-specific preference parameters $\varvec{\theta }=(\varvec{\theta }_{1},\varvec{\theta }_{2},\ldots ,\varvec{\theta }_{C})$ in the willingness-to-pay space, and the population share of each class $\varvec{\pi }=(\pi _{1},\pi _{2},\ldots ,\pi _{C-1})$. The population share of class $C, \pi _{C}$, is not a free parameter to be estimated since the class shares must add up to 1 and hence $\pi _{C}\equiv 1-\sum \nolimits _{c=1}^{C-1}\pi _{c}$. As discussed at the beginning of Sect. 4, the total number of classes C needs be set by the researcher: the common empirical practice is to estimate several specifications that vary in C, and choose one that results in the best Bayesian Information Criterion (BIC).

The sample log-likelihood function, $\ln L(\tau ,\varvec{\theta },\varvec{\pi })$, can be constructed in the usual manner by summing the log of each person’s likelihood, $L_{n}(\tau ,\varvec{\theta },\varvec{\pi })$ in Eq. (13):

$$\begin{aligned} \ln L(\tau ,\varvec{\theta },\varvec{\pi })=\sum \limits _{n=1}^{N}\ln L_{n}(\tau ,\varvec{\theta },\varvec{\pi })=\sum \limits _{n=1}^{N}\ln \left( \sum \limits _{c=1}^{C}\pi _{c}L_{n|c}(\tau ,\varvec{\theta }_{c})\right) . \end{aligned}$$

(18)

As explained in detail around Eq. (12), the kernel function $L_{n|c}(\tau ,\varvec{\theta }_{c})$ uses class c’s preference parameters to evaluate the baseline NROL model’s likelihood of observing person n’s actual sequence of rank orderings over choice scenarios. Since $L_{n|c}(\tau ,\varvec{\theta }_{c})$ has a closed-form expression, the sample log-likehood $\ln L(\tau ,\varvec{\theta },\varvec{\pi })$ also has a closed-form expression. It is therefore possible to proceed in the usual manner to compute the maximum likelihood estimates, say $\{\tau ^{ML},\varvec{\theta }^{ML},\varvec{\pi }^{ML}\}$, by applying gradient-based optimization methods (e.g. Newton method or quasi-Newton methods like BFGS) to maximize $\ln L(\tau ,\varvec{\theta },\varvec{\pi })$ with respect to $\{\tau ,\varvec{\theta },\varvec{\pi }\}$. But, as Bhat (1997) and Train (2008) note in the context of latent class multinomial logit models, maximizing the sample log-likelihood of a finite mixture model tends to be susceptible to convergence failures which often arise as a numerical optimizer gets trapped in flat regions of the sample log-likelihood function.

It turns out that parametric values $\{\tau ^{ML},\varvec{\theta }^{ML},\varvec{\pi }^{ML}\}$ that maximize $\ln L(\tau ,\varvec{\theta },\varvec{\pi })$ should, in theory, also maximize another objective function $Q(\tau ,\varvec{\theta },\varvec{\pi })$ and vice versa. The EM algorithm is an estimation strategy that aims at obtaining $\{\tau ^{ML},\varvec{\theta }^{ML},\varvec{\pi }^{ML}\}$ by maximizing this alternative objective function which is specified as

$$\begin{aligned} Q(\tau ,\varvec{\theta },\varvec{\pi })= & {} \sum \limits _{n=1}^{N}\sum _{c=1}^{C}h_{nc}(\tau ,\varvec{\theta },\varvec{\pi })\times \left( \ln \pi _{c}+\ln L_{n|c}(\tau ,\varvec{\theta }_{c})\right) \nonumber \\= & {} \sum \limits _{n=1}^{N}\sum _{c=1}^{C}h_{nc}(\tau ,\varvec{\theta },\varvec{\pi })\ln \pi _{c}+\sum \limits _{n=1}^{N}\sum _{c=1}^{C} h_{nc}(\tau ,\varvec{\theta },\varvec{\pi })\ln L_{n|c}(\tau ,\varvec{\theta }_{c}). \end{aligned}$$

(19)

$h_{nc}(\tau ,\varvec{\theta },\varvec{\pi })$ refers to person n’s posterior probability of membership in class c or $h_{nc}$ defined in Eq. (14), but we change the notation slightly here to emphasize that it is derived from the set of parameters being estimated. $Q(\tau ,\varvec{\theta },\varvec{\pi })$ may be interpreted as an expected complete data log-likelihood that views a set of indicators, $1[\varvec{\theta }_{n}=\varvec{\theta }_{c}]$ for $n=1,2,\ldots ,N$ and $c=1,2,\ldots ,C$, as missing data. Maximizing $Q(\tau ,\varvec{\theta },\varvec{\pi })$ is computationally easier and more numerically stable than maximizing $\ln L(\tau ,\varvec{\theta },\varvec{\pi })$ directly because, as we shall summarize shortly, $Q(\tau ,\varvec{\theta },\varvec{\pi })$ can be maximized with respect to $\{\tau ,\varvec{\theta }\}$ and $\varvec{\pi }$ in two separate tasks. While the use of the EM algorithm to estimate a finite mixture model is common outside discrete choice modeling too, it is Bhat (1997) who introduced this estimation strategy into the discrete choice modeling literature. Train (2008, (2009) masterfully summarizes the conceptual foundations and operational aspects of the EM algorithm for discrete choice models.

Our implementation of the EM algorithm builds on Bhat (1997) and Train (2008, (2009). Specifically, let superscript s denote candidate estimates obtained at the $s^{th}$ iteration of this algorithm. Then, at iteration $s+1$, the estimates are updated as follows.

$$\begin{aligned} \{\tau ^{s+1},\varvec{\theta }^{s+1}\}= & {} \arg \max _{\{\tau ,\varvec{\theta }\}}\sum \limits _{n=1}^{N}\sum _{c=1}^{C}h_{nc}(\tau ^{s},\varvec{\theta }^{s},\varvec{\pi }^{s})\ln L_{n|c}(\tau ,\varvec{\theta }_{c}) \end{aligned}$$

(20)

$$\begin{aligned} \varvec{\pi }^{s+1}= & {} \arg \max _{\varvec{\pi }}\sum \limits _{n=1}^{N}\sum _{c=1}^{C}h_{nc}(\tau ^{s},\varvec{\theta }^{s},\varvec{\pi }^{s})\ln \pi _{c} \end{aligned}$$

(21)

Each person’s posterior class membership probabilities are evaluated at the sth estimates, thereby influencing computation of the $(s+1)$th estimates only as any other type of known sampling weight would be. Both (20) and (21) therefore represent relatively simple maximization tasks. The algebraic structure of task (20) is just like that of estimating the baseline NROL model (not the LC-NROL model) using C years of data, allowing for year-specific coefficients and accounting for sampling weights: it can be readily solved by a maximum likelihood estimation program coded for the baseline NROL model. Our own program for estimating the baseline NROL model initially uses the BHHH method, a quasi-Newton method which is the default optimizer for maximum likelihood estimation in TSP International, and double-checks convergence by supplying the BHHH solution as starting values for executing the Newton method. The second updating task (21) is even easier to solve since it does not require any numerical optimization. An analytic solution to (21) can be derived and coded as

$$\begin{aligned} \pi _{c}^{s+1}=\frac{\sum _{n=1}^{N}h_{nc}(\tau ^{s},\varvec{\theta }^{s},\varvec{\pi }^{s})}{\sum _{n=1}^{N}\sum _{l=1}^{C} h_{nl}(\tau ^{s},\varvec{\theta }^{s},\varvec{\pi }^{s})}\quad \text {for}\quad c=1,2,\ldots ,C \end{aligned}$$

(22)

using that $\pi _{C}\equiv 1-\sum \nolimits _{c=1}^{C-1}\pi _{c}$.

Once starting values are provided for initial estimates at $s=0$, the EM algorithm proceeds by repeatedly updating the candidate estimates as above until $\Delta \ln L^{s+1}=\ln L(\tau ^{s+1},\varvec{\theta }^{s+1},\varvec{\pi }^{s+1})-\ln L(\tau ^{s},\varvec{\theta }^{s},\varvec{\pi }^{s})$ is small enough. Our own program uses the following set of starting values. To initialize the dissimilarity coefficient and class shares, we assume no within-nest correlation and equal class sizes i.e. $\tau ^{0}=1$ and $\pi _{c}^{0}=1/C$ for all $c=1,2,\ldots ,C$. To initialize class-specific preference parameters, $\varvec{\theta }^{0}=(\varvec{\theta }_{1}^{0},\varvec{\theta }_{2}^{0},\ldots ,\varvec{\theta }_{C}^{0})$ , we randomly partition the sample into equally sized C subsamples and estimate the ROL model on each subsample: then, the ROL estimates from the $c^{th}\,$subsample are used as starting values for class c’s parameters, $\varvec{\theta }_{c}^{0}$. Our program takes the $(s+1)^{th}$ estimates as the final estimates if $\Delta \ln L^{s+1}$ is smaller than $0.0001\,\%$ of $\ln L(\tau ^{s},\varvec{\theta }^{s},\varvec{\pi }^{s})$: call these final estimates $\{\tau ^{EM},\varvec{\theta }^{EM},\varvec{\pi }^{EM}\}$.

There are two drawbacks inherent in the EM algorithm as implemented here. First, as one may infer from the updating tasks (20) and (21), it does not produce valid standard errors of the final estimates, $\{\tau ^{EM},\varvec{\theta }^{EM},\varvec{\pi }^{EM}\}$. Second, while $\{\tau ^{EM},\varvec{\theta }^{EM},\varvec{\pi }^{EM}\}$ are equivalent to $\{\tau ^{ML},\varvec{\theta }^{ML},\varvec{\pi }^{ML}\}$ in theory, they may diverge in practice since the EM algorithm may declare convergence prematurely, or even in case the model is empirically unidentified: its stopping criterion is based on $\Delta \ln L^{s+1}$ and does not execute checks on the gradient and Hessian of $\ln L(\tau ^{s+1},\varvec{\theta }^{s+1},\varvec{\pi }^{s+1})$. To obtain standard errors and double-check convergence, therefore, we have followed the hybrid estimation strategy of Bhat (1997) that uses $\{\tau ^{EM},\varvec{\theta }^{EM},\varvec{\pi }^{EM}\}$ as starting values for direct maximization of $\ln L(\tau ,\varvec{\theta },\varvec{\pi })$. The main manuscript reports the resulting estimates $\{\tau ^{ML},\varvec{\theta }^{ML},\varvec{\pi }^{ML}\}$ and associated standard errors. In our experience, this direct maximization step always achieves convergence within a very small number of iterations, since the use of a stringent stopping criterion ($0.0001\,\%$ change in the sample log-likelihood) ensures that $\{\tau ^{EM},\varvec{\theta }^{EM},\varvec{\pi }^{EM}\}$ are almost identical to $\{\tau ^{ML},\varvec{\theta }^{ML},\varvec{\pi }^{ML}\}$ in practice as well as in theory. Since starting values $\{\tau ^{EM},\varvec{\theta }^{EM},\varvec{\pi }^{EM}\}$ tend to be close to the final solution $\{\tau ^{ML},\varvec{\theta }^{ML},\varvec{\pi }^{ML}\}$, the use of a quasi-Newton method does not bring in practical benefits and our own program immediately uses the Newton method for this direct maximization step.

Appendix 2

See Tables 6 and 7.

Table 6 Comparison of LC-NROL with mixed ROL estimation results

Table 7 Individual characteristics and definitions