Scalable Estimation of Multinomial Response Models with Random Consideration Sets

Scalable Estimation of Multinomial Response Models with Random Consideration Sets

\def\spacingset#1{ {#1}} \spacingset{1}

\newtheorem{lemma}{Lemma} \newtheorem{claim}{Claim} \newtheorem{assumption}{Assumption} \newtheorem{prop}{Proposition} \newtheorem{theorem}{Theorem}

{

}

}

{\it Keywords: Multinomial response, Bayesian computation, Dirichlet process mixture, Markov chain Monte Carlo, Metropolis-Hastings algorithm, Posterior consistency}

\spacingset{1.88}

Introduction

A common assumption when fitting unordered multinomial response models, whether applied to cross-sectional or longitudinal data, is that the responses stem from the same set of \( J \) mutually exclusive categories across all subjects. However, this assumption may be questionable, especially when modeling the choices made by human subjects. For example, in fields such as economics and marketing, it is recognized that individuals may select from only a subset of the available alternatives, termed the “consideration set" (Manski1977; HonkaHortacsuWildenbeest2019handbook). Neglecting this heterogeneity in the consideration sets can result in biased parameter estimates in the model (BronnenbergVanhonacker1996JMR; ChiangChibNarasimhan1998; Goeree2008; DraganskaKlapper2011; DeLosSantos2018IJIO; Morozov2021MrkSci; CrawfordGriffithIaria2021). Such biases are problematic because these models are typically employed to understand the impact of covariates on outcomes and inform decision making.

In order to fix ideas, let $\mathcal{C}_i$ represent the latent consideration set for subject $i$. When $J$ alternatives are available, $\mathcal{C}_i$ is a subset of $\{1, \ldots, J\}$, and there are $J^* = 2^J - 1$ possible consideration sets. A priori, $\mathcal{C}_i$ is assumed to be drawn from a probability mass function $\Pr(\mathcal{C}_i = c)$. When $J$ is small, the direct approach proposed by ChiangChibNarasimhan1998 is effective. In this approach, all possible consideration sets $1, 2, \ldots, J^*$ are enumerated and assigned unknown probabilities $\pi_1, \pi_2, \ldots, \pi_{J^*}$, which can be estimated using MCMC methods under a Dirichlet prior. However, when $J$ is large, the model has traditionally been estimated under the assumption that the distribution over consideration sets is determined by $J$ independent attention probabilities. In this framework, it is assumed that each alternative appears independently in any given consideration set BenAkivaBoccara1995, Goeree2008, ManziniMariotti2014, KawaguchiUetakeWatanabe2021, AbaluckAdams2021. Specifically, let $q_{ij}$ denote the probability that subject $i$ considers the alternative $j$ for $j = 1, \ldots, J$. The probability that $\mathcal{C}_i = c$ given $\boldsymbol{q}_i=(q_{i1},\ldots,q_{iJ})'$ is then modeled as: \[ \Pr\left( \mathcal{C}_i = c \mid \boldsymbol{q}_i \right) = \prod_{j \in c} q_{ij} \prod_{j \notin c} \big(1 - q_{ij}\big). \] Although this model is appealing for handling the large $J$ case, the distribution over consideration sets is unrealistic and leads to model misspecification CrawfordGriffithIaria2021.

In another approach, the consideration sets are modeled as vectors of 0-1 binary variables vanNieropPaap2010. This vector is then modeled by a multivariate probit model (AlbertChib1993JASA, ChibGreenberg1998MVP). Although this can generate correlation of items in consideration sets, inference is challenging because the number of parameters in the correlation matrix of the multivariate probit model increases quadratically in $J$.

Given the significant interest in incorporating consideration set heterogeneity in various fields - such as marketing vanNieropPaap2010, Ching2014simple, KawaguchiUetakeWatanabe2021, Turlo2025discrete, economics Goeree2008, Ching2009price, KashaevPeerEffect2019peer, AgarwalSomaini2022Restud, transportation science SwaitBenAkiva1987TransportationB, PaletiTransportationLetters2021, and psychology FritsPsychology2022 - there is a pressing need to develop a scalable estimation approach for estimating such generalized multinomial response models. The importance of accounting for consideration set heterogeneity becomes even more critical as $J$ increases, which is precisely the case that current methods struggle to address. The method we propose is based on two key components. The first component is a representation of the probability masses $\pi_{1}, \pi_{2}, \ldots, \pi_{J^{\ast }}$ in terms of a weighted average of products of item-specific inclusion $q_j$ and exclusion $1-q_j$ probabilities, which is based on a result from DunsonXing2009. We refer to this approach as a mixture of independent consideration models. To simulate the latent consideration sets, we introduce a straightforward and intuitive Metropolis-Hastings algorithm. It is important to highlight that, in this context, the consideration sets are latent, unlike in Dunson and Xing (2009), where the categorical variables are observed. This difference necessitates additional steps in both the theoretical derivations and the computational procedure. Another crucial feature of the method is the sparsity of the posterior distribution of the consideration sets, which occurs because sets that do not include the actual choices made by a subject must have a posterior probability of zero ChiangChibNarasimhan1998. The scalability of the proposed approach is demonstrated through an application to marketing data involving $J=101$ brands.

We establish two key theoretical results. First, under regularity conditions, as the number of subjects increases, we demonstrate that the posterior distribution of the marginal response probabilities is consistent. Second, under certain additional identification assumptions, the posterior distribution of the model parameters also achieves consistency.

In general, this paper contributes to the expanding literature on high-dimensional demand estimation in statistics and marketing: (BraunMcAuliffe2010variational; ChiongShum2019MS; SmithAllenby2019JASA; LoaizaNibbering2022ScalableProbit_JBES; Jiang2024high_MS; WangIaria2024_JEEA; Ershov2024RAND; Amano2018large). To incorporate latent consideration sets, it is necessary to generalize the standard multinomial response model by conditioning the distribution of responses on a latent subject-specific consideration set, which is drawn from the power set of \(\{1,2,\ldots,J\}\). This results in a mixture of multinomial models based on a probability distribution over consideration sets. However, the exponential size of this power set renders the estimation of this mixture of multinomial response models computationally infeasible in general. Moreover, the proposed method can be interpreted as a generalized multinomial logit (MNL) model, with “structural zeros” incorporated in the first layer of its hierarchical structure. In the field of biostatistics, methodologies have been extensively explored to estimate microbial compositions that account for the sparsity due to excessive zero counts (e.g.\ Aitchison1982; Martin2015; Liu2020empirical; Cao2020multisample; Paulson2013differential; Chen2016two; Tang2019zero). More recently, Zeng2023zero introduced a zero-inflated probabilistic PCA model designed for high-dimensional, sparse microbiome data sets. Although our paper focuses on a different problem, the proposed method has the potential to be applied in similar contexts, as we discuss in the concluding section.

The remainder of the article is structured as follows. Section (ref) introduces the model, while Section (ref) presents the theoretical results. Section (ref) discusses posterior inference and computational methods. Section (ref) reports numerical simulations, and Section (ref) applies the methodology to a marketing dataset. Finally, the concluding section explores the broader implications of the proposed framework.

The approach

Suppose that we have panel (longitudinal) data with $n$ a priori independent subjects that contains multinomial (polychotomous) responses from a set $\mathcal{J}=\{1,\ldots ,J\}$ of $J$ mutually exclusive nominal categories/items as well as some covariates. Let $Y_{it}\in \mathcal{J}$ be the measured response for unit $i$ at time $t$, where $i=1,\ldots,n$ and $t=1,\ldots,T_i$. Let $\bm w_{it}=\{ \bm w_{ijt}\}_{j \in \mathcal{J}}$, where $\bm w_{ijt}$ is the vector of covariates characterizing the category $j$ for subject $i$ at time $t$. Each subject $i$ is associated with a latent consideration set $\mathcal{C}_i$, which is a subset of the entire set of alternatives $\mathcal{J}$. We model the distribution of the observed outcomes using a hierarchical approach. Specifically, we first specify the marginal distribution of the consideration sets and then define the response distribution conditional on a given consideration set. In this framework, we make the following assumptions.

Assumption 1: Consideration sets $\mathcal{C}_i$ vary over subjects but not over time, and the distribution over consideration sets, denoted by $\pi_c = \Pr(\mathcal{C}_i = c)$ for $c \in \mathcal{C}$, the set of all possible consideration sets minus the empty set, is free of covariates.

The assumption of time invariance is relatively mild and aids in inference. It also plays a role in the identification of model parameters. Covariates can be included in the model for consideration sets, but, as noted by ChiangChibNarasimhan1998, a covariate-dependent model is difficult to specify without increasing the risk of model mis-specification.

Assumption 2: For each $j \in \mathcal{J}$, the responses $Y_{it}$ of subject $i$ given $\mathcal{C}_i$ and random effects $\bm b_i$ are independent over time and follow the multinomial logit model.

Based on Assumptions 1 and 2, the generalized multinomial logit model of interest has the hierarchical form:

for $i = 1, \ldots, n$, where $\bm \pi = \{\pi_c : c \in \mathcal{C}, 0 \leq \pi_c \leq 1, \sum_{c \in \mathcal{C} } \pi_c =1 \}$ denotes the collection of probabilities associated with all possible consideration sets, and $\bm b_i$ are random effects normally and independently distributed across subjects with zero mean and unknown covariance matrix $\bm{D}$. The covariates are denoted by $\bm{w}_{it} = \{ \bm{x}_{ijt}, \bm{z}_{ijt} \}_{ j \in \mathcal{J}}$, where $\bm{x}_{ijt} \in \mathbb{R}^{d_{x}}$ and $\bm{z}_{ijt} \in \mathbb{R}^{d_{z}}$. Stage 1 can be interpreted as introducing another layer of random effects, where heterogeneity arises from the random consideration sets.

Letting $\Pr (\bm Y_{i} \vert \bm \theta, \bm w_{i}, \mathcal{C}_{i}=c)$ denote the distribution of outcomes $\bm Y_i = (Y_{1i},\ldots,Y_{iT_i})$ of subject $i$ marginalized over the random effects given covariates $\bm w_i=\{\bm w_{i1},\ldots,\bm w_{iT_i}\}$, the distribution of responses takes the finite mixture form: \[ \Pr (\bm Y_{i} \vert \bm \theta, \bm w_{i}) = \sum_{c\in \mathcal{C}} \pi_c \Pr (\bm Y_{i} \vert \bm \theta, \bm w_{i}, \mathcal{C}_i=c). \] This can be seen as a generalized multinomial logit response model.

Assumptions 1 and 2 imply time-invariant consideration sets, conditional independence of responses, and full support of the conditional response probabilities given consideration sets. These conditions, along with additional assumptions detailed below, establish the point identification of the model parameters AguiarKashaev2024identification in the model that excludes random effects. Furthermore, in Theorem 2 of Section (ref), we show the posterior consistency of the parameters in this case.

The latent consideration sets

To fix notation, let $\mathcal{C}$ represent the collection of all possible consideration sets, which corresponds to the power set of $\mathcal{J} = \{1, \ldots, J\}$, excluding the empty set. The consideration set for subject $i$ is indicated by $\mathcal{C}_i = c$, where $c \in \mathcal{C}$. For example, when $J = 3$, $\mathcal{C} = \{\{ 1 \}, \{2 \}, \{3\}, \{1,2\}, \{1,3\},\{2,3\}, \text{and} \, \{1,2,3\} \}$, and $c$ is one of these elements. Furthermore, by $\bm C_i=(C_{i1},\ldots,C_{iJ})'$, we mean a $J \times 1$ multivariate binary vector where $C_{ij}=1$ if category $j$ is in the consideration set, and $0$ otherwise. In the example of $J=3$, $\mathcal{C}_i=\{1\}$ is equivalent to $\bm C_i=(1,0,0)'$ and $\mathcal{C}_i=\{1,3\}$ is equivalent to $\bm C_i=(1,0,1)'$ etc. In the following, we use the two notations interchangeably depending on the context. Researchers sometimes include an outside option in the model that is always considered by each subject. We can incorporate this into our framework by adding a $(J+1)$th category and fixing $C_{iJ+1}=1$ for all $i$. Our goal is to put a probability distribution on $\mathcal{C}$ that is rich enough to accommodate dependencies while maintaining scalability.

Dimensionality reduction via tensor decomposition

We now review the factor decomposition technique that we employ to specify the distribution over consideration sets. DunsonXing2009 consider modeling large contingency tables that, for example, represent DNA sequences, each of which is defined as a collection of $J$ categorical variables, each having $d_j$ possible values $j=1,\ldots,J$, where $J$ is large. A realization of the contingency table can be expressed as a vector $(a_1,\ldots,a_J)'$, where $a_j \in \{1,\ldots,d_j\}$ for $j=1,\ldots,J$. The true distribution of the contingency tables is a probability tensor $\bm \pi= \{ \pi_{a_1a_2\cdots a_J}, a_j=1,\ldots,d_j, j=1,\ldots,J \}$, where $0 \leq \pi_{a_1a_2\cdots a_J} \leq 1$ and $\sum_{a_1=1}^{d_1}\cdots \sum_{a_J=1}^{d_J} \pi_{a_1a_2\cdots a_J} = 1$. Note that consideration sets can be seen as contingency tables with $d_j=2$ for all $j$. Generally, there are a large number of elements in the tensor $\bm \pi$, $d_1\times \cdots \times d_J$, when $J$ is large. DunsonXing2009 show that $\bm \pi$ can be expressed as a finite mixture of rank 1 tensors. We describe this result for the special case that corresponds to modeling consideration sets.

This result states that a mixture of $K$ independent consideration models can model an arbitrary distribution over the $J^{\ast} = 2^J - 1$ possible consideration sets. Within each component $h$, items are included in or excluded from a consideration set $c$ according to an independent consideration model defined by a vector of attention probabilities $\bm q_h=(q_{h1},\ldots,q_{hJ})'$. Therefore, the number of parameters needed to model the probabilities in $\bm{\pi}$ is reduced from $J^{\ast}$ to $K \times J + (K - 1)$, which scales linearly with $J$.

Infinite mixture of independent consideration models

Building on this result, we model the $J$-dimensional latent vectors \( \{ \bm{C}_i \} \) as a mixture of independent probabilities. Since the number of components \( K \) in (ref) is unknown, we follow DunsonXing2009 and use a Dirichlet process (DP) prior Ferguson1973BayesianNonparametrics to induce an infinite mixture model. One key difference from DunsonXing2009 is that their categorical variables (contingency tables) are observed, while the corresponding consideration sets are latent. This difference leads to differences in the theoretical analysis (Section (ref)) and in the posterior simulation approach (Section (ref)).

In our approach, we do not estimate $K$. This is because existing methods for consistently estimating $K$, such as those proposed by KwonMbakop2021AoS, may not be applicable when the variables modeled by the mixture are latent. Posterior consistency in our framework only requires that the prior on $K$ has positive mass for all positive integers. Posterior inferences on model parameters and their functions (e.g., predictions) automatically account for uncertainty regarding the value of $K$.

Assume that $\{\bm C_i\}$ is i.i.d. with density $f({}\cdot{} \vert G) = \int \prod_{j=1}^J q_j^{C_{ij}} \left( 1-q_j\right)^{1-C_{ij}} dG(\bm{q})$. The discrete mixing distribution $G$ is modeled by a DP prior with a concentration parameter $\alpha$ and a specified base probability measure $G_0$ that depends on a hyperparameter $ \underline{\bm \phi}_q$. Equivalently, by using the stick breaking construction (Sethuraman1994constructiveDP), we have the following representation: $\bm C_i$'s are i.i.d.\ with the density for the infinite mixture of independent consideration models:

where $\bm c_i=(c_{i1},\ldots,c_{iJ})'$, $\omega_1= V_1, \omega_h = V_h \prod_{\ell<h} (1-V_\ell), h=2,\ldots,\infty,$ $V_h\overset{iid}{\sim}\text{Beta}(1,\alpha)$, and $\bm q_h \overset{iid}{\sim} G_0(\ \cdot \ \vert \underline{\bm \phi}_q), h=1,\ldots,\infty$, with $\bm q_h=(q_{h1},\ldots,q_{hJ})'$ being the vector of attention probabilities specific to the component $h$. A priori, the first few weights dominate and cover most of the probability mass, which are then adjusted by the data. Although the model (ref) includes infinitely many components, typically only a small number of distinct values for $\bm{q}_h$ are imputed.

For the baseline distribution $G_{0}$, we assume that $q_{hj} \sim G_{0j}$ independently for $j=1,\ldots,J$ and $h=1,\ldots,\infty$. Specifically, we assume that $q_{hj}\sim \text{Beta}(\underline{a}_{q_j},\underline{b}_{q_j})$, independently over $j=1,\ldots,J$, for $h=1,\ldots,\infty$, and we define $\underline{\bm \phi}_q=(\underline{\bm a}_q,\underline{\bm b}_q)$ with $\underline{\bm a}_q =(\underline{a}_{q_1} ,\ldots,\underline{a}_{q_J})'$ and $\underline{\bm b}_q =(\underline{b}_{q_1},\ldots,\underline{b}_{q_J})'$. Note that $\underline{\bm \phi}_q=(\underline{\bm a}_q,\underline{\bm b}_q)$ are the hyperparameters chosen by the user. We discuss this in more detail in the Supplementary Material. We complete the model specification by assuming the prior distribution for the DP concentration parameter $ \alpha \sim \text{Gamma}(\underline{a}_\alpha, \underline{b}_\alpha), $ where $(\underline{a}_\alpha, \underline{b}_\alpha)$ are the hyperparameters chosen by the user. For smaller values of $\alpha$, $\omega_h$ decreases toward zero more rapidly as $h$ increases, so that the prior favors a sparse representation with most of the weight on a few components. We allow the data to inform about $\alpha$ and, therefore, an appropriate degree of sparsity.

Theoretical results

We establish two key results. For simplicity, let $T_i=T$, $\forall i$ and suppose that $T\geq 1$ is fixed and $n\to \infty$. In Theorem (ref) we show that the posterior of the marginal response probabilities is consistent, and in Theorem (ref), we show that the posterior of the model parameters is consistent when $T$ is large enough and the model does not include random effects.

Let $\bm \theta=\{ \bm \beta, \bm D \}$ denote the parameters in the response model. Also, recall that the distribution over the consideration sets is denoted by $\bm \pi=\{ \pi_c: c \in \mathcal{C}\}$, where $0\leq \pi_c \leq 1$ and $\sum_{c \in \mathcal{C}}\pi_c=1$. Define the probability that the sequence of items $\bm y=(y_1,\ldots,y_T)'\in \mathcal{J}^T$ is chosen conditional on covariates $\bm w_i =\{ \bm w_{i1},\dots,\bm w_{iT} \}$ taking some specific value $\bm w =\{\bm w_1,\ldots,\bm w_T\} \in \mathbb{R}^{TJ(d_x+d_z)}$: \[ p_{\bm{\theta}, \bm{\pi} } (\bm y\vert \bm w) \equiv \sum_{c\in \mathcal{C}} \pi_{c} \Pr\left( \bm Y_{i} = \bm y \vert \bm \theta, \bm w, c \bm \right), \] where the response probability given a consideration set $c$ is \[ \Pr (\bm Y_{i}= \bm y \vert \bm \theta, \bm w, c) = \int \prod_{t=1}^T \Pr(Y_{it} = y_t \mid \bm{\beta}, \bm{w}_{t}, \mathcal{C}_i=c, \bm{b}_i) \phi(\bm b_i \vert \bm 0, \bm D) d \bm b_i , \] where the integrand is defined in (ref). The data set contains responses $\bm y_i=\{y_{it}\}$ and covariates $\bm w_i=\{\bm w_{it} \}$ and we let $\bm D^n=\{(\bm y_i,\bm w_i): i=1,\ldots,n \}$. The covariates $\bm w_i $ are i.i.d.\ and follow an unknown distribution with density $g^*$ with support $\mathcal{W}\subset \mathbb{R}^{TJ(d_x+d_z)}$. We do not model the covariate distribution. Conditional on covariates, responses are generated from the collection of the data-generating response probabilities $\bm p^* =\{p_{\bm{\theta}^*, \bm{\pi}^* } (\bm y \vert \bm w) \}_{\bm y \in \mathcal{J}^T,\bm w \in \mathcal{W}}$, where $\bm \theta^*$ denotes the true response model parameter and $\bm \pi^*=\{ \pi^*_c: c \in \mathcal{C}\}$ denotes the true probability mass function over consideration sets. We emphasize that $\bm \pi^*$ does not have to be a finite mixture. The joint probability measure implied by $\bm p^*$ and $g^*$ is denoted by $F_0$. For $\varepsilon>0$, define a Kullback-Leibler neighborhood of $\bm p^*$ as \[ KL_\varepsilon(\bm p^*) = \left\{ (\bm \theta, \bm \pi): \int \sum_{\bm y \in \mathcal{J}^T} \log \left( \frac{ p_{\bm{\theta}^*, \bm{\pi}^* } (\bm y\vert \bm w)}{ p_{\bm{\theta}, \bm{\pi} } (\bm y\vert \bm w) } \right) p_{\bm{\theta}^*, \bm{\pi}^* } (\bm y\vert \bm w) g^* (\bm w) d\bm w < \varepsilon \right\}. \] It is essentially a set of $(\bm \theta,\bm \pi)$ that makes $p_{\bm{\theta}, \bm{\pi}}$ close to $p_{\bm{\theta}^*, \bm{\pi}^*}$.

Given a $K\in \mathbb{Z}^+$, define $\bm \phi_{1:K}=\{\omega_h,\bm q_h: h=1,\ldots,K\}$, the collection of all component-specific parameters, where $\bm q_h=(q_{h1},\ldots, q_{hJ})'$. Note that by Lemma (ref), there exist $\{K, \tilde{\bm \phi}_{1:K}\}$, which may not be unique, such that $ \pi_c^* =\sum_{h=1}^{K} \tilde{\omega}_h \left\{ \prod_{j \in c} \tilde{q}_{h j} \prod_{ j \notin c} \big(1-\tilde{q}_{h j} \big) \right\},$ $\text{ for all } c\in \mathcal{C}, $ and the KL divergence is zero at $\{\bm \theta^*, K, \tilde{\bm \phi}_{1:K} \}$. In the following lemma, we establish that the KL divergence can be made arbitrarily small in sufficiently small neighborhoods of $(\bm \theta^*, \tilde{\bm \phi}_{1:K})$. Define the model induced probability for a consideration set $c\in \mathcal{C}$: $ \pi(c\vert K, \bm \phi_{1:K}) = \sum_{h=1}^K \omega_h \prod_{j\in c} q_{hj} \prod_{j \notin c}(1- q_{hj}), $ and the model induced marginal response probability as \[ p(\bm y \vert \bm w; \bm \theta, K, \bm \phi_{1:K}) =\sum_{c\in \mathcal{C}} \pi(c \vert K, \bm \phi_{1:K}) \Pr(\bm Y_{i}= \bm y \vert \bm \theta, \bm w, c). \]

The proof can be found in the Appendix. Let $\Pi(\cdot)$ denote the prior for the response model parameter $\bm \theta$ and the distribution of consideration sets $\bm \pi$.

This result shows that the model-induced response probability in the limit converges to the true data-generating process. A similar result is proved in DunsonXing2009, Theorem 2, but for the case in which the categorical variables are observed and there are no covariates. Because our setup relaxes both of those conditions, we have a more involved proof that involves the KL-divergence (Lemma (ref)). NoretsShimizu2024 also establish a related result for semiparametric dynamic discrete choice models, but our proof strategy is different, due to the random effects, continuous covariates, and a different model. Last, the compactness assumption (ii) is common in Bayesian nonparametric estimation, and condition (iii) of Theorem 1 is satisfied by our DP prior for $\omega_h$'s and the Beta prior for $q_{hj}$'s, following Dunson and Xing (2009).

We now address the possibility that multiple parameter pairs \((\bm{\theta}, \bm{\pi})\) may be consistent with the true response probabilities. This relates to the issue of partial identification MasatliogluNakajimaOzbay2012revealed, CattaneoMaMasatliogluSuleymanov2020, BarseghyanCoughlinMolinariTeitelbaum2021, Lu2022, where point identification holds only under specific conditions DardanoniManziniMariottiTyson2020ECMA, AbaluckAdams2021, BarseghyanMolinariThirkettle2021discrete. Following AguiarKashaev2024identification, we impose the assumption that the panel is sufficiently long and that random effects are absent. Under these conditions, we show that the two sources of variation in responses—differences in utility and differences in consideration sets—can be separately identified. Formally, we show in the next theorem that the posterior distribution contracts to within an arbitrarily small ball around $(\bm \beta^*, \bm \pi^*)$ under the distance function $d((\bm \beta, \bm \pi), (\bm \beta', \bm \pi'))=\max \{||\bm \pi - \bm \pi'||_1, ||\bm \beta-\bm \beta'||_2 \}$.

We remark that in Theorem 2 we suppose a model without random effects, though we use random effects in our modeling. The complication in having both is that latent consideration sets in our model operate similarly to random effects and introduce dependence across time. Disentangling these two sources of dependence at a theoretical level requires a stronger condition on $T$, though the precise details are not straightforward to establish. We leave this extension for future work. Nonetheless, the numerical experiments in the Supplementary Material indicate that the convergence described in the theorem holds more generally, as we observe convergence to the true values even in the presence of random effects.

Inference

Let $\bm{Y}_i=(Y_{i1},\ldots,Y_{iT_i})'$ and $\bm{y}_i=(y_{i1},\ldots,y_{iT_i})'$ be the sequence of random responses made by unit $i$ over $T_i$ periods and its observed counterpart. Define

where $\bm w_i = \{\bm w_{i1},\ldots, \bm w_{iT_i}\}$ and $\Pr(Y_{it}=y_{it} \vert \bm \beta, \bm{b}_i, \bm w_{it},\mathcal{C}_i )$ is

Note that $\bm C_i$ is the conditioning variable on the left side of (ref), while $\mathcal{C}_i$ is on the right side. Although the two objects represent the same information, the $J$-dimensional vector $\bm C_i$ is easier to use when we discuss posterior sampling of individual consideration sets. Hence, we use $\bm C_i$ to define the individual's contribution to the likelihood. Let $\bm Y=\{\bm{Y}_1,\ldots,\bm{Y}_n\}$ and $\bm y=\{\bm{y}_1,\ldots,\bm{y}_n\}$ denote the random and observed sequences of the responses made by all units, and let $\bm W=\{\bm w_1,\ldots,\bm w_n\}$ be the observed covariates. Then the likelihood conditional on the common fixed-effects $\bm \beta$, the random effects $\bm b = (\bm{b}_1,\ldots,\bm{b}_n)'$, the covariates $\bm W$, and the latent consideration sets $\bm C=(\bm{C}_1,\ldots,\bm{C}_n)$ is given by

We complete the model by specifying standard prior distributions for the parameters in the response model: $\bm{\beta} \sim \mathcal{N}_{d_x}\left( \bm{0}, \underline{\bm{V}}_\beta \right)$ and $\bm{D}^{-1} \sim \text{Wishart} \left( \underline{v}, \underline{\bm{R}}\right)$, indepdently, a normal distribution for $\bm{\beta}$, and an inverse Wishart distribution for $\bm D$ with degrees-of-freedom parameter $\underline{v}$ and scale matrix $\underline{\bm{R}}$. The hyperparameters $( \underline{\bm{V}}_\beta, \underline{v}, \underline{\bm{R}})$ are chosen by the user.

Posterior distribution

For the mixture model on the latent consideration sets $\bm C=(\bm{C}_1,\ldots,\bm{C}_n)$, let $S_i \in\{1,2,\ldots\}$ be the latent cluster assignment such that $C_{ij} \vert S_i=h\sim \text{Bernoulli}(q_{hj})$, independently $j=1,\ldots,J$, for $i=1,\ldots,n$. We have the latent consideration sets $\bm C$, the common fixed-effects $\bm \beta$, the random effects $\bm b$, the corresponding covariance matrix $\bm D$, the DP parameters $\bm V =(V_1,V_2,\ldots)$ as well as $\bm Q=(\bm q_1,\bm q_2,\ldots)$, the DP cluster assignment variables $\bm S =(S_1,\ldots,S_n)$, and the DP concentration parameter $\alpha$. Let $\pi(\cdot)$ denote the prior density. Then, from the Bayes theorem, the posterior density of interest is

where the first term is given by (ref) and only the last term is associated with the DP prior.

We sample from the posterior distribution using a tailored Markov Chain Monte Carlo (MCMC) algorithm. The method is designed for scalability and consists of simple and intuitive steps. Posterior inference is then based on the sampled values

where $G$ is the number of MCMC draws beyond a suitable burn-in period.

Simulation of consideration sets

We now focus on sampling the conditional distribution of consideration sets. The other steps in the MCMC simulation follow from standard calculations and are given in the Supplementary Material. From Equation (ref), the full conditional distribution of $\bm C_i$ is

where the proportionality sign is with respect to $\bm C_i$, and the first term is defined in (ref). Importantly, consideration sets that exclude any observed response made by subject $i$ receive zero posterior probability (see Table (ref) for an example). This is because the first term on the left-hand side of (ref) is zero for these consideration sets. This desirable feature of our approach is based on ChiangChibNarasimhan1998. In contrast, in many existing methods, every consideration set receives a strictly positive probability, as pointed out by CrawfordGriffithIaria2021. Now, due to the independence structure in (ref) over $j=1,\ldots,J$, \[ \pi( C_{ij} \vert \bm C_{i}\setminus \{ j \}, \bm \beta,\bm b_i,\bm q_{S_i}, S_i,\bm y_i, \bm w_i ) \propto \ p\big(\bm Y_{i}=\bm y_{i} \big\vert \bm \beta,\bm b_i, \bm w_i, \bm C_i \big) \cdot q_{S_ij}^{C_{ij}} (1-q_{S_ij})^{1-C_{ij}} , \] where $\bm C_{i}\setminus \{ j \}$ denotes $\bm C_i$ without the coordinate $j$. To sample from this distribution, we employ the Metropolis-Hastings (M-H) algorithm ChibGreenberg1995understandingMH. An effective implementation of this approach is detailed in Algorithm (ref). \FloatBarrier

In Step 1 of Algorithm (ref), we generate a proposal from a one-dimensional Bernoulli distribution. In Step 2, given the current state $\bm{C}_i^{(0)}$ and the proposed state $\bm{C}_i^{(1)}$, the acceptance probability is computed as the ratio of the likelihood contributions for subject $i$. This Metropolis-Hastings step is valid because the likelihood $p\big(\bm Y_i = \bm y_i \mid \bm \beta, \bm b_i, \bm w_i, \bm C_i \big)$ is uniformly bounded. See ChibGreenberg1995understandingMH (p.\ 330, the third algorithm) for more discussion. In practice, we update the states in a random order within each MCMC iteration. In addition, the computational burden is minimized by parallelizing the loop on the $n$ subjects.

Finally, the proposed Metropolis-Hastings step exhibits an important sparsity property. Suppose that an alternative $j$ was not chosen by the subject $i$ in any period (otherwise, it must be in the consideration set for $i$ and $C_{ij}=1$). Depending on the current $C_{ij}^{(g)}$, and the proposed $\tilde{C}_{ij}$, there are four possible moves in the M-H step. First, if $\tilde{C}_{ij}=C_{ij}^{(g)}=1$ or $\tilde{C}_{ij}=C_{ij}^{(g)}=0$, then the proposed value is accepted with probability one. Second, if $\tilde{C}_{ij}=0$ and $C_{ij}^{(g)}=1$, then the proposed value is also accepted with probability one. In other words, the algorithm “prefers" a smaller consideration set. This sparsity-inducing property is proven below. Lastly, when the proposed consideration set adds an alternative $j$ that is not in the current consideration set, that is, $\tilde{C}_{ij}=1$ and $C_{ij}^{(g)}=0$, the acceptance probability is between 0 and 1 and is determined by the likelihood ratio.

Numerical illustration

We illustrate posterior probabilities of consideration sets on synthetic panel data with $n=100$ subjects observed over $T\in \{1,2,\ldots,15 \}$ time periods. We let $J=4$ and give the $2^{J}-1=15$ consideration sets in the first column of Table (ref). In the table we report the posterior probabilities of each possible consideration set for a randomly chosen subject $i$ whose true consideration set is $\mathcal{C}_{i}^*=\{1,3,4 \}$.

The first column ($T=1$) shows the results for the initial period given the observed outcome of $1$. Consideration sets that do not include item 1 have a posterior probability of zero. As $T$ increases, the posterior concentrates on the true consideration set $\{1,3,4\}$.

Monte Carlo Simulation

We demonstrate the sampling performance of the proposed approach through simulation studies first with $J=4$ alternatives, where it is possible to enumerate all the support points in $\bm \pi$, and then extend the study to a high-dimensional case with $J = 100$. The goal is to empirically validate the findings of Theorem 2 and demonstrate that the proposed approach can effectively assess consideration dependence. In the Supplementary Material, we conduct additional experiments under autocorrelated covariates, random effects, and time-varying true consideration sets. In general, the experiments show posterior consistency in the estimation of $\bm\theta=(\bm \beta, \bm D)$ and $\bm \pi$, and that the restrictive approach with $K=1$ produces larger root mean squared errors and biases.

$J = 4$

We let $T_i = T$ for all $i$. In one case we set $T=5$ and in the other $T=15$. The latter satisfies the length condition of Theorem (ref). In simulating the data, we first specify the distribution of the consideration sets $\bm \pi^*=\{ \pi^*_c = \Pr(\mathcal{C}_i=c): c \in \mathcal{C}\}$. We induce dependence in product consideration by letting the first two and last two products have a relatively high probability of being considered together: $\pi^*_{\{1,2\}}=\pi^*_{\{3,4\}}=0.25$. As motivation, the first two products might represent non-vegetarian options, and the last two vegetarian. The other 13 consideration sets $c \in \mathcal{C}$ are given a probability of 0.0385. Figure (ref) shows $\bm \pi^*$ in red. Given this $\bm \pi^*$, we generate the true consideration sets $ \mathcal{C}_i^*$, for $i=1,\ldots,n$. We then generate outcomes from the logit model with $V_{ijt} =\delta^*_j+ \beta^* x_{ijt}, $ letting $(\delta_1^*,\delta_2^*,\delta_3^*,\delta_4^*)'=(1.0, 0.5, -1.0, 0)'$ and $\beta^*=1$, and $x_{ijt} \overset{iid}{\sim} N(0,1)$. We let $n\in\{50,100\}$.

We compare the performance between the proposed infinite mixture of independent consideration models ($K=\infty$) and the model that assumes independent consideration ($K=1$) over 200 replicated data sets. The results are given in Table (ref) where we report the root mean squared error (RMSE) for the response parameter $\bm \beta=(\delta_1,\delta_2,\delta_3, \beta)$ as well as the $L_1$ norm between the posterior mean and the truth for the distribution of the consideration sets $\bm \pi$ (L1-error), their Monte Carlo errors (MCE), the posterior standard deviation (SD), the empirical standard deviation (ESD), the empirical coverage of the equal-tailed 95% credible intervals (Cov), and the computational time. The MCE quantifies the precision for the performance criterion. The MCEs are negligible, allowing for valid comparisons based on the $200$ replications. As $n$ increases, the posterior of $\bm \beta$ and $\bm \pi$ contracts to the true values, even when $T = 5$, indicated by the smaller RMSEs and L1-errors as well as SDs. In contrast, when $K=1$, we do not observe sufficient evidence of posterior consistency. The RMSEs and L1-error are much larger in some cases than those under $K=\infty$, due to misspecification. When $T$ increases to $15$, a value that satisfies the identifying condition in Theorem (ref), the RMSEs/L1-errors/SDs become smaller for both $K=\infty$ and $K=1$, but for $K=1$, they are larger, and there are distortions in the coverage. Finally, our approach ($K=\infty$) delivers good coverages in general. The SDs are similar to ESDs, indicating that the posterior standard deviations provide a good representation of the sampling variability of the posterior means.

The vertical axes of Figure (ref) list the 15 consideration sets, with the true distribution of the consideration sets, $\bm \pi^*$, highlighted in red. Each panel of the figure displays the posterior mean (solid with dots, blue) along with the 95% credible intervals (dashed, blue), based on one realized data set. The first two panels illustrate that under the proposed approach ($K=\infty$), as the sample size $n$ increases, the discrepancy between the posterior mean and the true distribution diminishes. In contrast, the right two panels show that when $K=1$, even as $n$ increases, the posterior does not adequately converge to the truth. This is because the model does not account for the true consideration dependence.

$J=100$

We now consider a high-dimensional scenario with $J=100$ alternatives. One mechanism by which the dependence of consideration among categories can be induced is through multiple latent subpopulations of subjects having different probabilities of consideration. Within a subpopulation, considerations are independent across categories. However, marginalizing out the latent subpopulation indicator, one obtains dependence in those category considerations. We generate the data with two subpopulations. To generate the true consideration set of a given subject, we used a Bernoulli distribution with attention probability 0.05 for each category except for categories 10, 30, 50, 70, and 90 for the first subpopulation ($i=1,\ldots,n/2$) where the attention probability was set to 0.8. For the remaining subjects in the second subpopulation ($i=n/2+1,\ldots,n$), the Bernoulli probability was set at 0.05 except for categories 20, 40, 60, 80, and 100 where the probability was set to 0.8. Conditional on the true consideration sets, we generated the responses as in the case with $J=4$ with $\delta_j^*=0, \ j=1,\ldots,J-1$.

Because in this case there are $2^{100}-1$ support points in $\bm \pi$, it is not possible to show the entire distribution as in the case of $J=4$. Also, there are 99 $\delta_j$'s to estimate. Hence, in Table (ref), we report the results only for the slope $\beta$ as well as $\delta_{97}, \delta_{98},$ and $\delta_{99}$. The results are based on 200 replications. The general observations from the small $J$ simulation still hold: as $n$ increases, the RMSEs/SDs of $\bm \beta$ become smaller even when $T=5$, supporting posterior consistency. For $T = 200$, a span that approximately satisfies the identification condition in Theorem (ref), our approach also results in good coverages.

\FloatBarrier

\FloatBarrier

Application to Cereal Consumption in Midwest

In this section, we apply our approach to a manually constructed longitudinal data set that includes $J=101$ cereal brands, a size that is significantly beyond the feasibility of existing methods. For comparison, $J$ was 4 in ChiangChibNarasimhan1998, 10 in vanNieropPaap2010, and 5 in AguiarKashaev2024identification. We constructed the data set by integrating Nielsen Consumer Panel data with Retail Scanner Data, focusing on weekly shopping trips in 2019 in stores operated by a single anonymous retailer primarily based in the United States Midwest. Although data from 2020 are available, we chose to use the most recent pre-pandemic year to avoid potential biases introduced by pandemic-related shopping behavior. This particular retailer was selected because it consistently stocked more than 100 cereal brands throughout the sample period. Furthermore, we limited our analysis to a single retailer to prevent inconsistencies in brand definitions between different retailers, which would have required speculative alignment of brand names from various sources. The final data set includes $J=101$ brands and $n=1880$ households, covering 25,849 purchases in 239 stores during the 52-week period in 2019. See Figure (ref), for the locations of these stores with relative purchase volumes, and Table (ref), for the list of the brands. The average number of shopping trips per household ($T_i$) is 13.7, and the price $P_{ijt}$ of each brand $j \in {1,\ldots,J}$ is represented by a size-weighted price index constructed from prices at the UPC level. For the analysis, we used the first 10 months of data for the estimation and reserved the last two months for the prediction outside the sample. Further details on data preparation are provided in the Supplementary Material.

Conditional on the consideration set $\{\mathcal{C}_i\}$, in the most general version of the model, we enter the fixed effects and random effects in the MNL model $V_{ijt} =\delta_j + P_{ijt}(\beta + b_i)$, where $i \in \{1, \ldots, 1880\}$ indexes households, and $t \in \{1, \ldots, T_i\}$ indexes purchase occasions. In this model, $\delta_j$ represents the brand-specific fixed effect for brand $j$, with the normalization $\delta_J = 0$. The parameter $\beta$ is the common fixed effect, and $b_i \sim \mathcal{N}(0, D)$ is the random effect for household $i$. We consider four variants of the MNL, differentiated by the inclusion of random effects and/or consideration set heterogeneity, as detailed in models (1)–(4) of Table (ref). In addition, models (5) and (6) assume an independent consideration structure (i.e.,\ $K=1$). Each of these cases is estimated using the simulation method developed in Section (ref), by omitting the components not present in the full hierarchical model (MNL_RC).

Empirical Results

We obtained 20,000 MCMC draws for each of the six models in Matlab on a desktop with a 4.9GHz processor and 64GB RAM. The average of the inefficiency factors is around 7.38 with standard deviation 2.41, indicating that the MCMC output mixes well. Broadly speaking, the estimated parameters of the response model from the approaches (1)-(4) shown in Table (ref) are similar to those in the literature. For instance, when consideration set heterogeneity is incorporated, the magnitude of the slope parameter $\beta$ on price increases and the number of significant brand-specific terms $\delta_j$'s decreases (See Table (ref) for the list of estimated $\delta_j$'s under the MNL_RC model). These patterns are consistent with previous studies based on smaller models, including Stopher1980, SwaitBenAkiva1986, ChiangChibNarasimhan1998, and vanNieropPaap2010. However, without the scalable fitting methodology developed in this paper, it was unclear if those patterns would persist in a model of the scale we have estimated.

Moreover, when we control for consideration sets, the posterior mean of $D^{1/2}$ decreases, which aligns with the findings in ChiangChibNarasimhan1998, Morozov2021MrkSci that random effect heterogeneity is overestimated in models that omit consideration set heterogeneity.

Under the independent consideration assumption ($K=1$), i.e.,\ (5) and (6), the estimated parameters are similar to the proposed flexible approach i.e.,\ (3) and (4) except that the estimated $D^{1/2}$ under (6) is slightly larger than (2), which contradicts with the previous studies. In general, it is possible that the obtained estimates under $K=1$ are biased, as shown in simulation studies in Section (ref). We conduct the test for independent consideration, which is introduced and studied in Supplementary Material. Under both (3) and (4), the estimated posterior probability of the alternative hypothesis (dependent consideration) is very close to one, and we conclude that the considerations of cereal products in this particular market are dependent.

Table (ref) also shows the computational time per 1,000 MCMC draws. The extra burden of estimating latent consideration sets using our proposed approach is reasonable. For instance, when consideration sets are estimated along with random effects, the computational time roughly doubles (67 mins.\ for MNL_R and 124 mins.\ for MNL_RC). Not surprisingly, compared to the fully flexible estimator, the estimators that assume independent consideration take less computational time but only slightly.

Estimated parameters in the mixture model

We begin by reporting in Table (ref) the posterior mean and standard deviation (s.d.) of the 100 brand fixed effect parameters. The 95% posterior credibility intervals of most of these brand-specific intercepts exclude zero indicating that these brands are endowed with significant brand equity. We next investigate the clustering of households according to the proposed mixture model. The posterior mode of the number of nonempty clusters under the full-specification (MNL_RC) is six. The Supplementary Material shows further estimation results on the number of clusters and the DP concentration parameter $\alpha$.

To understand how households are clustered, we computed the posterior mean of the event that a given pair of households $(i,k)$ are clustered together i.e.\ $\{S_i = S_{k}\}$. This results in a $n \times n$ similarity matrix, which can be found in the Supplementary Material.

An examination of how households are clustered reveals interesting points. Take household A as an example whose actual choices consist of $\{4,37,62,64,73\}$. Define an estimator $\hat{\mathcal{C}}_i$ of the consideration set for household $i$ as the set of brands $j$ whose posterior probability that $C_{ij}=1$ is greater than 0.2658, the prior median of $ q_{hj}$. This results in the estimated set $\hat{\mathcal{C}}_{A}=\{4,26,37,59,60,61,62,64,68,73,79 , 101\}$. The upper panel of Table (ref) lists the three households with the highest posterior similarity to subject $A$. There are several observations.

First, the actual choices of the households tend to overlap within a cluster; each household purchased at least one of brands 5, 73, or 89. Second, the estimated consideration sets $\hat{\mathcal{C}}_{k}$ are similar between households in a cluster. For example, household A did not choose brands 26, 79, and 101, but other households did, and they are in $\hat{\mathcal{C}}_{A}$. Third, the stronger the purchase overlap, the higher the chance of being in the same cluster. The lower panel of Table (ref) shows the results for household B. In this cluster, brands 79 and 81 were purchased by all the four households, brands 26, 51, 68, and 101 were each purchased by three households, and we see higher similarity scores ($\geq 0.60$). In this way, our algorithm discovers the probabilistic grouping patterns in the choice data.

Price sensitivity of demand

To analyze household shopping behavior, we randomly select 100 units and report in Figure (ref)

the percentage decrease in aggregate demand when the price of a brand increases by 1% under the MNL_R and MNL_RC models. For all but one brand, this sensitivity is higher under consideration set heterogeneity, in conformity with previous findings that were derived in a small $J$ setting.

Predictive performance

We next assess the predictive performance of the proposed model using the last two months of data as an out-of-sample period. Let $\mathcal{O}\subset \{1,\ldots, n\}$ denote the set of subjects who made purchases in the out-of-sample period. This set contains 1079 subjects. For each $i \in \mathcal{O}$, we predict $\bm Y^f_i = \{ Y_{iT_i+s}: s=1,\ldots,h_i \}$, given the covariates $\bm w^f_i = \{ \bm w_{iT_i+s}: s=1,\ldots,h_i \}$, where $h_i$ denotes the forecast horizon for the subject $i$. Let $\bm y^f_i=\{ y_{iT_i+s}: s=1,\ldots,h_i \}$ be the actual set of responses for the subject $i \in \mathcal{O}$. Then, as a measure of predictive performance, we calculate the predictive likelihoods

where the response probability conditional on a consideration set is given in (ref). Figure (ref) gives the log-predictive likelihood for each household under the (MNL_R) and (MNL_RC) models. The higher predictive likelihood under the latter model shows that including consideration set heterogeneity tends to improve predictive performance. More details about this are given in the Supplementary Material.

Discussion

In this concluding section, we discuss the broader relevance of the work, especially to the modeling of excess zeros in high-dimensional sparse microbiome data sets. In a microbiome dataset with $n$ samples and $J$ taxa, let $u_{ij}$ denote the measured count for taxon $j$ in sample $i$, and $T_i = \sum_{j=1}^J u_{ij}$ represent the total count over taxa in the $i$th sample, where $i = 1, \ldots, n$ and $j = 1, \ldots, J$. Typically it is assumed that the vector $\bm{u}_i = (u_{i1}, \ldots, u_{iJ})'$ follows a multinomial distribution with index $T_i$ and a vector of probabilities $\bm{\rho}_i = (\rho_{i1}, \ldots, \rho_{iJ})'$, where $0 < \rho_{ij} < 1$ and $\sum_{j=1}^J \rho_{ij} = 1$. In our notation, $u_{ij}$ relates to $Y_{it}$ through $u_{ij}=\sum_{t=1}^{T_i}1(Y_{it}=j)$. To address the high dimensionality and sparsity of such datasets, Zeng2023zero propose the following hierarchical model translated in our terminology as

where $\bm C_i = (C_{i1},\ldots,C_{iJ})'$, and $1-C_{ij}$ are latent indicators for excess zeros, and the $q_{j}$ are the corresponding probabilities. The $\bm{f}_i$ are latent factors, and the $\bm{\beta}_j$ denote the loadings of the associated factors. Note that the excess zeros across taxa are independent in this modeling. In other words, the model above corresponds to the independent consideration model that we review in the introduction. In this context, complex dependency patterns across taxa in the excess zeros can be captured by our modeling. To do this, we would let $\bm C_i$ be correlated vectors and i.i.d.\ with the density for the infinite mixture of independent consideration models (3): \[ \Pr(\bm C_i = \bm c_i )= \sum_{h=1}^\infty \omega_h \prod_{j=1}^J \left\{ q_{hj}^{c_{ij}} \left( 1-q_{hj}\right)^{1-c_{ij}} \right\}. \] We can adapt the MCMC framework of this paper to estimate this model. Our approach for updating the $\bm C_i$'s and the mixture parameters can be used in conjunction with existing approaches for simulating the factor-related objects.

Another key issue in practice is variable selection when many subject-level covariates are available. This challenge can be addressed using shrinkage priors. A natural extension of our framework involves modeling consideration sets that change at one or two points in time due to learning from past choices. This would require incorporating the learning process into the model and modifying the theoretical analysis accordingly. We leave this promising direction for future work.