The benefits of incorporating utility dependencies in finite mixture probit models

Abstract

We propose an application of a new finite mixture multinomial conditional probit (FM-MNCP) model that accommodates preference heterogeneity and explicitly accounts for utility dependencies between choice alternatives considering both local and background contrast effects. The latter is accomplished by using a one-factor structure for segment-specific covariance matrices allowing for nonzero off-diagonal covariance elements. We compare the model to a finite mixture multinomial independent probit (FM-MNIP) model that as well accommodates heterogeneity but assumes independence. That way, we address the potential benefits of a model that additionally accounts for dependencies over a model that accommodates heterogeneity only. Our model comparison is based on empirical data for smoothies and is assessed in terms of fit, holdout validation, and market share predictions. One of the main findings of our empirical study is that allowing for utility dependencies may counterbalance the effects of considering heterogeneity, and vice versa. Additional findings from a simulation study indicate that the FM-MNCP model outperforms the FM-MNIP model with respect to parameter recovery.

Appendix

1.1 Model identification

To ensure identification of the estimated models we proceeded as follows: First of all, we used 20 choice sets per respondent consisting of pairwise disjunct alternatives (except for the base alternative) to prevent respondent fatigue on the one hand, but to provide as much respondent information as possible for model estimation on the other hand. Eighteen choice sets were used for estimation, while two choice sets served as holdouts for validation. Having three real alternatives plus the base (none) alternative in each choice set, $20 \cdot 3 +1 = 61$ pairwise different alternatives were displayed to a respondent over the whole choice task. Following Haaijer et al. (1998), the number of covariance parameters in the FM-MNCP model must not exceed the total number of alternatives as a necessary condition for identification. According to this requirement, the chosen study design would therefore enable the identification of the FM-MNCP model for up to five segments. More specifically, since with three attributes at four levels each a total of ten covariance parameters need to be estimated per segment (formally one less than the number of levels for each attribute plus one for the base option) the inequality $18 \cdot 3+1 \ge 10 \cdot G$ holds for up to $G=5$ segments. As sufficient condition for model identification Hesse matrices must further be checked for full rank and negative definiteness. We computed the Hesse matrices for the estimated models, respectively, with the result that Hessian matrices were negative definite in all cases. Therefore, since both the necessary and the sufficient conditions were satisfied, identification of the estimated models could be guaranteed.

For conjoint studies with a larger number of attributes and/or a larger number of latent segments, and therefore a larger number of covariance parameters to be estimated in the FM-MNCP model, the necessary condition for model identification seems more critical as it is. On the one hand, the researcher can increase the number of alternatives per choice set which lies typically between 3 and 5 (excluding the none option). This way, the number of different alternatives involved in the whole choice task can be increased. On the other hand, as it is common for the specification of interaction terms between attributes in practical applications, it may also be possible to limit the occurrence of context effects to specific pairs of attributes. In this case, the number of covariance parameters can be reduced. In addition, a recent meta-analysis on applications of finite mixture models indicates that commonly a moderate number of segments (on average three to four) is sufficient to describe empirical data adequately (Tuma and Decker 2013).

1.2 Model estimation

The EM algorithm consists of I iterations and incorporates a Gibbs sampling with R iterations in the E-step and the BFGS algorithm in the M-Step.

According to established stopping criteria for the EM algorithm, we terminate the EM algorithm in iteration i if changes in parameters are less than an a priori specified bound (cf. Wedel and Kamakura 2000, p. 88 or Train 2008, p. 44). Using several different starting points, every estimation was started with an initial vector $\theta _0$, whose $\Theta$ values are drawn component by component from the standardized normal distribution. Several pretests revealed an appropriate number of $R=800$ iterations for the Gibbs sampling algorithm in the E-Step and a maximum number of 5000 iterations for the BFGS algorithm in the M-Step. We used the BFGS algorithm provided in the “optim”-function coded in R. The EM algorithm was stopped if none of the estimated parameters changed more than 0.03 over the last three iterations.

1.2.1 Estimation: E-step

The utility difference vectors $z_k$ are $(M-1)\cdot P$-variate Gaussian distributed. A closer examination reveals that this Gaussian distribution is truncated on $(-\infty , 0)^{(M-1)\cdot P}$. According to Xu and Craig (2009), p. 4, the E-step consists of a Gibbs sampling that cycles through the conditional distributions. Specifically, we run R iterations of the following two steps in the ith iteration of the EM algorithm:

• 1. Let $g_k\in \left\{ 1,\ldots ,G\right\}$ describe the membership of the kth choice pattern to segment g. $g_k$ is multinomial distributed and can therefore be drawn from a multinomial distribution with the following G parameters: For $g=1,\ldots ,G$

  $$\begin{aligned} p_{g_k=g}^{(ir)}=\frac{\pi ^{(i-1)}_{g^*} \cdot \phi _{(M-1)P}\left( z_k^{(i-1)}; A_k X \beta _{g^*}^{(i-1)}, A_k\varSigma _{g^*}^{(i-1)}A_k^\mathrm{T}\right) }{\sum \nolimits _{g=1}^G\pi ^{(i-1)}_{g} \cdot \phi _{(M-1)P}\left( z_k^{(i-1)}; A_k X \beta _{g}^{(i-1)}, A_k\varSigma _{g}^{(i-1)}A_k^\mathrm{T}\right) }, \end{aligned}$$

  (10)

  and $\sum _{g=1}^{G} p_{g_k=g}^{(ir)} =1$ for given i and r.

• 2. Based on the segment memberships $g_k^{(ir)}$ resulting from (1.), the corresponding utility difference vectors $z_k^{(ir)}$ can be determined. The difference vectors $z_k^{(ir)}=(z_{k1}^{(ir)},\ldots , z_{k((M-1)P)}^{(ir)} )^\mathrm{T}$ are drawn component by component q, $q=1,\ldots ,(M-1)\cdot P$, from an univariate truncated Gaussian distribution $TN(\mu _q^{(ir)}, \sigma _q^{(ir)})$ on $(-\infty ,0)$, where the computation of $\mu _q^{(ir)}$ and $\sigma _q^{(ir)}$ follows the approach of Xu and Craig (2009):

  $$\begin{aligned} \sigma _q^{(ir)}= \frac{1}{\left( A_k\varSigma _g^{(i-1)}A_k^\mathrm{T}\right) ^{-1}_{[q,q]}}, \end{aligned}$$

  (11)

  where [q, q] denotes the qth diagonal element of the $(M-1)\cdot P \times (M-1)\cdot P$ matrix $(A_k\varSigma _g^{(i-1)}A_k^\mathrm{T})^{-1}$ and

  $$\begin{aligned} \mu _q^{(ir)}= \left( A_k X \beta _g^{(i-1)}\right) _{[q]}-\sigma _q^{(ir)}\left( A_k\varSigma _g^{(i-1)}A_k^\mathrm{T}\right) ^{-1}_{[q,-q]} (z_{k[-q]}^{(i-1)}-\left( A_k X \beta _g^{(i-1)})_{[-q]}\right) ,\!\!\!\!\!\nonumber \\ \end{aligned}$$

  (12)

  where [q] denotes the qth component of the vector, $[-q]$ the corresponding vector without the qth component, and $[q,-q]$ the qth row of the corresponding matrix without the qth element.

After R iterations of the Gibbs sampling procedure and under consideration of a burn-in phase of R / 2 iterations (also see Xu and Craig 2009), the conditional expectations of the segment memberships and the utility difference vectors in iteration i of the EM algorithm are calculated as follows (for all $g=1,\ldots ,G$):

$$\begin{aligned} E\left( g_k^{(i)} = g|z_k, y_k, \beta _g, \varSigma _g\right)= & {} \frac{2}{R} \cdot \sum \limits _{r=R/2}^R g_k^{(ir)}\ and \nonumber \\ E\left( z_k^{(i)}|g_k = g, y_k, \beta _g, \varSigma _g\right)= & {} \frac{2}{R} \cdot \sum \limits _{r=R/2}^R z_k^{(ir)} \end{aligned}$$

(13)

Hence, the relative segment masses $\pi _g^{(i)}$ can be derived as

$$\begin{aligned} \pi _g^{(i)}= \frac{1}{J} \cdot \sum \limits _{k=1}^{M^P}n_k\cdot E\left( g_k^{(i)}=g|z_k, y_k, \beta _g, \varSigma _g\right) . \end{aligned}$$

(14)

1.2.2 Estimation: M-step

The conditional expectations $E(z_k^{(i)}|g_k=g, y_k, \beta _g, \varSigma _g)$ of the utility difference vectors $z_k$ and the relative segment masses $\pi _g^{(i)}$ build the input for the M-step and are used as fixed quantities during the maximization of the log likelihood. Based on those quantities, the BFGS algorithm determines the segment-specific part-worth vectors $\beta _{g}^{(i)}$ (and covariance vectors $\sigma _g^{(i)}$) in iteration i of the EM algorithm.

1.3 Simulation study

To gain further evidence on the performance of FM-MNCP versus FM-MNIP models, we conducted a simulation study and compared both models with respect to model fit, predictive validity, as well as parameter recovery. Parameter recovery as a further performance criterion can only be addressed with simulated data since the true parameters are known here (as opposed to empirical data). The simulation study manipulates commonly used experimental factors for finite mixture models, i.e., number of segments, separation between segments, and relative segment masses (cf. Vriens et al. 1996; Andrews et al. 2002b; Andrews and Currim 2003). In addition, we included the covariance structure (zero versus nonzero structure) as another experimental factor. A zero covariance structure implies full independence between alternatives and therefore coincides with the independence assumption of the FM-MNIP model.

The simulated data comprise individual choice patterns for 600 respondents as well as artificially generated “true” part-worth utilities and covariance vectors (the latter only for the FM-MNCP model). The data generation process leans on the approaches of Vriens et al. (1996) and Andrews and Currim (2003). For model comparison, we computed the model-specific means of several performance measures both per factor level and across factor levels. Model fit is measured by the log likelihood, while parameter recovery is assessed by the root mean square error RMSE$(\beta ) =(1/(G\cdot S) \sum _{g=1}^G \sum _{s=1}^S (\beta _{gs} - \widehat{\beta }_{gs})^2 )^{1/2}$ between the true segment-specific part-worths $\beta _{g}$ and the estimated segment-specific part-worths $\widehat{\beta }_{g}$. Predictive validity is evaluated by the root mean square error $\hbox {RMSE}(V) =(1/(J\cdot W\cdot M) \sum _{j=1}^J \sum _{w=1}^W \sum _{m=1}^M(V_{jwm} - \widehat{V}_{jwm})^2 )^{1/2}$ between the true deterministic utilities of holdout alternatives $V_{jwm}$ and the estimated ones $\widehat{V}_{jwm}$ (cf. Andrews et al. 2002b; Andrews and Currim 2003), where W denotes the number of holdout choice sets consisting of M alternatives each. To be comparable to our empirical study, the choice tasks consisted of 18 choice sets for estimation and 2 holdouts for validation, 4 alternatives including the none option per choice set, and 3 attributes with four levels for designing the alternatives. With our four experimental factors at two factor levels each and one additional replication for each of the resulting 16 treatments, we obtained a total of 32 observations for each of the performance statistics and model type.

Table 8 displays the means of all three performance measures both at the individual factor level and overall across factor levels.

Table 8 Simulation study: model-specific means of performance measures for model fit, parameter recovery, and predictive validity

Across factor levels, the FM-MNCP model significantly ($\alpha \le 5\%$) outperforms the FM-MNIP model with respect to model fit ($\hbox {LL}^\mathrm{MNCP}=-10{,}396.289>-10{,}458.533=\hbox {LL}^\mathrm{MNIP}$) and parameter recovery ($\hbox {RMSE}(\beta )^\mathrm{MNCP}=0.529 < 0.671= \hbox {RMSE}(\beta )^\mathrm{IP}$) (see row “overall”). Concerning predictive accuracy the difference in the root mean square error statistic is not significant between models, although the FM-MNCP model shows a better (overall) $\hbox {RMSE}(V)$ value ($\hbox {RMSE}(V)^\mathrm{MNCP}=29.209 < 31.513 = \hbox {RMSE}(V)^\mathrm{IP}$).

Considering the performance of the models at individual factor levels, the main results can be described as follows: For two segments, the FM-MNCP model provides a significantly better predictive accuracy than the FM-MNIP model ($\alpha \le 5\%$). Parameter recovery is always better under the FM-MNCP model, and differences to the FM-MNIP model in $\hbox {RMSE}(\beta )$ are significant at all factor levels except for a zero covariance structure and equal segment masses ($\alpha \le 10\%$). For three segments as well as for a large separation of segments the FM-MNCP model significantly ($\alpha \le 10\%$) outperformed the FM-MNIP model with regard to model fit, too.

In contrast to our empirical study, the simulation study is based on comparisons of FM-MNCP and FM-MNIP models that were estimated for the same number of segments. For example, if the true number of segments was two, we estimated and compared the two-segment FM-MNCP solution to the two-segment FM-MNIP solution. This procedure is common practice in simulation studies (see, e.g., Andrews et al. 2002b). Obviously, the FM-MNIP model performs worse than the FM-MNCP model in many situations when the number of segments is held fixed a priori. Under those conditions, the FM-MNIP model seems not able to overcome its independence assumption by accounting for heterogeneity only. Probably, the differences between the two models may turn out still larger with a higher power of the design (i.e., when more replications under each treatment are used). For a discussion of the simulation study in full length, see Paetz and Steiner (2014) or Paetz and Steiner (2015).

1.4 Utility functions

Fig. 1

Segment-specific utility functions for the one- and two-segment solutions of the FM-MNCP and FM-MNIP models