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.
86,712 characters · 18 sections · 52 citation commands
Iterative Distributed Multinomial Regression
Discrete choice models, including logit and multinomial-logit (MNL) models, are widely used in applied social science research. With the growing availability of diverse data types and the integration of econometric models with textual and image data, researchers often encounter applications of MNL models with a massive number of choices; see the applications discussed in Section (ref). In these cases, even if the number of parameters for each choice is small, the total number of parameters, which increases linearly with the number of choices, will be large due to large choice sets. As a result, maximum likelihood estimation (MLE) can become computationally intractable due to the high cost of solving an optimization problem for a large number of parameters.\footnote{MLE of the MNL model refers to the conditional MLE given the covariate and total counts.}
Researchers in diverse disciplines have explored various approaches and proposed numerous methods to numerically solve the MLE; see the numerical algorithms discussed in Section (ref). However, there is a notable lack of theoretical results ensuring the consistency or asymptotic efficiency of the estimators derived from these methods. This gap in the literature between numerical computation and statistical inference motivates the present study.
Our proposed estimator utilizes the multinomial-Poisson (MP) transformation, which reformulates the multinomial likelihood into a Poisson likelihood by incorporating individual fixed effects into the MNL model. When all the covariates are categorical, baker1994multinomial establishes that the multinomial likelihood and the Poisson likelihood produce identical estimates of the parameters in the MNL model and advocates the computational advantage of maximizing the Poisson likelihood. However, the MP transformation has also been employed in MNL models with continuous or mixed discrete and continuous covariates. For instance, gentzkow2019measuring note that they “... approximate the likelihood of {[}their{]} multinomial logit model with the likelihood of a Poisson model...”.
As the first contribution of this paper, we establish an equivalence result for all types of covariates, continuous, discrete, or mixed, justifying the application of the MP transformation in these cases. To accomplish this, we re-interpret the Poisson likelihood from the MP transformation as a conditional quasi-log-likelihood function given the covariates. Maximizing this function provides a quasi-maximum likelihood estimator (QMLE) for the MNL model. The equivalence result is established by demonstrating that the resulting QMLE is identical to the MLE of the MNL model.
While the QMLE is computationally more efficient than the MLE, it remains costly when applied to MNL models with large choice sets. To address this computational challenge, taddy2015distributed exploits an important feature of the quasi-log-likelihood function: for any given fixed effects, the function is additively separable in the parameters across different choices of the MNL model. taddy2015distributed proposes to estimate parameters for each choice separately at a specific value of the fixed effects and calls the resulting estimator distributed computing estimator. As noted in taddy2015distributed, however, his distributed computing estimator is inconsistent except in a few very special cases.
To regain consistency and asymptotic efficiency, we adopt the idea of iterative backfitting algorithms studied in pastorello2003iterative, dominitz2005some, and fan2015maximization to compute the QMLE of the parameters in the MNL model and fixed effects iteratively.\footnote{Similar iterative algorithms have also been developed for dynamic discrete games of incomplete information where directly computing the MLE using the nested fixed-point algorithm proves computationally infeasible, see aguirregabiria2002swapping,aguirregabiria2007sequential for a nested pseudo-likelihood (NPL) algorithm. kasahara2012sequential further analyze the conditions necessary for the convergence of the NPL algorithm and derive its convergence rate. } During each iteration, we first solve for the parameters of interest through distributed computing, given approximate values of the individual fixed effects. Then, we update estimates of the individual fixed effects based on their expressions derived from maximizing the quasi-log-likelihood function using the previous estimates of the model parameters. We call our estimator iterative distributed computing (IDC) estimator. This is the second and main contribution of the paper. The IDC algorithm is fast because of distributed computing, even when the choice set is large. We consider three IDC estimators based on three different initial estimators: a consistent estimator based on pairwise binomial logistic regression, taddy2015distributed's distributed computing estimator, and a maximum likelihood estimator assuming that the distribution of total counts is Poisson. The latter two estimators are inconsistent in general. All three initial estimators are fast to compute because they allow for distributed computing.
As the third contribution, we establish theoretical results on the consistency and asymptotic efficiency of all three IDC estimators. When the initial estimator is consistent, we show that the IDC estimator with any finite number of iterations is always consistent and is asymptotically efficient under an information dominance condition when the number of iterations diverges with respect to the sample size $n$ at $\log\left(n\right)$ rate. With inconsistent initial estimators, the IDC estimators are consistent and asymptotically efficient under a stronger contraction mapping condition when the number of iterations diverges with respect to $n$ at a polynomial rate.
When the number of choices is large, conducting inference via plug-in estimation of the variance matrix becomes infeasible. This is because the Fisher information matrix has a very large dimension, causing the computation of the inverse of its estimator to be both time-consuming and unreliable. The fourth contribution of the paper is that we propose a parametric bootstrap inference procedure and show its consistency. Because the IDC estimator is fast to compute, our inference procedure is computationally feasible.
Lastly, we conduct extensive simulations to study the finite sample performance of our estimator and inference procedure. We are particularly interested in the computational time of the IDC estimator and its accuracy compared to the maximum likelihood estimator. The simulation results show that the IDC estimator is very fast to compute, with a running time approximately linear in the number of choices. Compared to the maximum likelihood estimator, our estimator has a very similar mean squared error in all the different model settings but is much faster to compute when the number of choices is large. We also study the finite sample behavior of the proposed bootstrap inference procedure. The results suggest that the procedure achieves the correct size and is consistent.
\paragraph{Applications}
The proposed IDC estimator can be applied to study various economics and computer science topics such as text analysis, dimensionality reduction, spatial choice models, image classification (russakovsky2015imagenet), and video recommendation (davidson2010youtube).
Text analysis: The integration of text data into econometric models is increasingly prominent in economics. For example, baker2006investor analyze investor sentiment’s effect on stock returns, while chen2021sentiment explore how hedge funds capitalize on sentiment changes. Modeling text data often involves treating word counts as a multinomial distribution, as taddy2015distributed demonstrates using Yelp reviews to predict outcomes based on user and business attributes. gentzkow2019measuring use the distributed computing estimator in taddy2015distributed for the multinomial regression to measure polarization in Congressional speeches. NBERwTextSelection extend this approach, using Hurdle Distributed Multinomial Regression to backcast, nowcast, and forecast macroeconomic variables from newspaper text.
Dimensionality Reduction: Our estimator aids in dimensionality reduction for inverse multinomial regression, as taddy2013multinomial discusses. Instead of inferring sentiment from text, taddy2013multinomial’s approach estimates word distribution given sentiment. He introduces a score based on word frequencies and regression parameters, which is useful in forward-regression models.
Spatial Choice Models: High-dimensional choices also appear in spatial models. bucholztaxi2019 models taxi drivers’ location choices with a dynamic spatial search, reducing dimensionality via discretization. Similarly, pellegrinifotheringham2002 apply hierarchical discrete choice models to immigration, while bettman1979memory addresses brand choices in limited-option settings, proposing hierarchical selection for high-dimensional cases.
\paragraph{Numerical Algorithms}
To address the computational difficulty of solving the MLE of the MNL model, researchers have proposed several numerical methods to find approximate solutions. bohning1988monotonicity and bohning1992multinomial propose replacing the Hessian matrix in the Newton-Raphson iteration with its easy-to-compute global lower bound and show that the approximate solution converges with the number of iterations. Because the convergence rate depends crucially on the difference between the Hessian matrix and its lower bound, the algorithm can be slow to run for certain model parameters. Additionally, based on our simulation exercise, if the choice probabilities vary significantly across different choices with some being close to zero, the algorithm becomes unstable. In comparison, our IDC algorithm is stable in all the simulation settings. boyd2011distributed introduce an alternating direction method of multipliers, which reformulates the original optimization problem by introducing redundant linear constraints. gopal2013distributed propose a log concavity method, which replaces the log partition function of the multinomial logit with a parallelizable upper bound. recht2011hogwild, raman2016ds, and fagan2018unbiased study a stochastic gradient descent method, which uses random training samples to calculate the gradient at each iteration. Although these methods can be computationally efficient, to the best of the authors' knowledge, no consistency or asymptotic efficiency result has been shown in these works.
Penalization methods have also been introduced to the MNL regression and some of the numerical methods discussed above are adopted in solving the penalized MNL regression, see e.g., friedman2010regularization, simon2013blockwise, and nibbering2022multiclass, The proposed IDC procedure in this paper can be combined with the aforementioned algorithms to further improve the performance of penalized MNL regressions.
\paragraph{Organization of the rest of this paper}
The remainder of this paper is organized as follows. In Section (ref) we present a comprehensive overview of the multinomial logistic regression model and the MP transformation. Section (ref) introduces our iterative distributed computing estimator along with some initial values. In Section (ref) we provide the asymptotic theory of the iterative distributed computing estimator. Section (ref) contains the simulation results. Finally, with Section (ref) we conclude. Appendix (ref) collects the notations and equations used in the paper. All the technical proofs are provided in Appendix (ref). The codes for implementing the estimation and inference procedures are available \href{https://github.com/yigitokar/IDMR}{here}.
\paragraph{Notations}
Throughout the paper, we use index $i\in\left\{ 1,\ldots,n\right\} $ for individual, $j\in\left\{ 1,\ldots,p\right\} $ for covariate, and $k\in\left\{ 1,\ldots,d\right\} $ for unique choice. Boldfaced symbols such as $\boldsymbol{C}$ and $\boldsymbol{V}$ are used to denote vectors; while elements of the vectors are denoted by plain symbol such as $C_{k}$ and $V_{j}$. Denote $\sim$ as “equality up to a constant”, such that $f\left(\theta\right)\sim g\left(\theta\right)$ is equivalent to $f\left(\theta\right)=g\left(\theta\right)+h$, where $h$ is a constant relative to $\theta$.
Let $\boldsymbol{C}_{i}\in\mathbb{R}^{d}$ denote the random vector of counts on $d$ different choices for individual $i=1,...,n$, summing up to $M_{i}=\sum_{k=1}^{d}C_{ik}$. We use the random vector $\boldsymbol{V}_{i}\in\mathbb{R}^{p}$ to denote the covariate vector that includes a constant.
Consider a correctly specified multinomial-logit (MNL) model. The conditional probability mass function is given by the following:
where for $k=1,\ldots d$, we let $\eta_{ik}^{\ast}\equiv\boldsymbol{V}_{i}^{\prime}\boldsymbol{\theta}_{k}^{\ast}$ with unknown parameters $\boldsymbol{\theta}_{k}^{\ast}\equiv\left(\theta_{k1}^{\ast},\ldots,\theta_{kp}^{\ast}\right)^{\prime}$ and $\Lambda_{i}^{\ast}\equiv\sum_{k=1}^{d}e^{\eta_{ik}^{\ast}}$. For the identification, we set $\boldsymbol{\theta}_{d}^{\ast}=\mathbf{0}$. Let $\boldsymbol{\theta}^{\ast}\equiv\left(\boldsymbol{\theta}_{1}^{\ast\prime},\ldots,\boldsymbol{\theta}_{d}^{\ast\prime}\right)^{\prime}\in\Theta$ denote the parameter vector of interest. Throughout the paper, we use the superscript $\ast$ to indicate the true value of the unknown parameter. Denote the parameter space of $\boldsymbol{\theta}_{k}^{\ast}$ for $k=1,\ldots,d$ as $\varTheta_{k}$. We have that $\varTheta_{d}=\left\{ \mathbf{0}\right\} $ and $\Theta\equiv\prod_{k=1}^{d}\varTheta_{k}$.
In this paper, we focus on the case where $d$ is large (but fixed) such that directly solving for the maximum likelihood estimator is computationally costly. Applications include text corpora, where $\boldsymbol{C}_{i}$ represents the counts of $d$ different words/phrases in a text of $M_{i}$ words; browser logs, where $\boldsymbol{C}_{i}$ indicates the number of times a website among $d$ total websites is visited by an individual; and location choices, where among $M_{i}$ number of locations traveled by the driver, $\boldsymbol{C}_{i}$ contains the number of times each location, among $d$ different ones, is visited.
Given a random sample of size $n$, let $\boldsymbol{\eta}_{i}\equiv\left(\eta_{i1},\ldots,\eta_{id}\right)^{\prime}\equiv\left(\boldsymbol{V}_{i}^{\prime}\boldsymbol{\theta}_{1},\ldots,\boldsymbol{V}_{i}^{\prime}\boldsymbol{\theta}_{d}\right)$ and $\Lambda_{i}=\sum_{k=1}^{d}e^{\eta_{ik}}$ for $i=1,...,n$. Ignoring terms that are independent of the parameter $\boldsymbol{\theta}$, the conditional log-likelihood function given the covariate $\boldsymbol{V}$ and total count $M$ takes the following form:
Let $L_{C\mid V,M}\left(\boldsymbol{\theta}\right)$ denote the probability limit of $\frac{1}{n}l_{C\mid V,M}\left(\boldsymbol{\theta}\right)$. Denote $B\left(\boldsymbol{\theta},\varepsilon\right)$ as an open ball in $\Theta$ centered at $\boldsymbol{\theta}$ with radius $\varepsilon$. We make the following assumption throughout the paper.
Assumption (ref) (i) implies that $\boldsymbol{\theta}^{\ast}$ is identified as $\boldsymbol{\theta}^{\ast}=\arg\max_{\boldsymbol{\theta}\in\Theta}L_{C\mid V,M}\left(\boldsymbol{\theta}\right)$. Define the following objective function:
Based on ((ref)), the conditional maximum likelihood estimator of $\boldsymbol{\theta}^{\ast}$ is:\footnote{The definition of $\widetilde{\boldsymbol{\theta}}$ implicitly assumes that the solution to the minimization problem is unique. This can be shown to hold with probability approaching one by the identification of the model. See mcfadden1973conditional. The same result holds for all the estimators defined in the paper. We ignore such mathematical subtlety for the remainder of the paper to simplify the discussion. }
Solving the above optimization problem analytically is impossible. In addition, due to the potentially large dimension $d$, numerical algorithms such as the Newton-Raphson method are difficult to implement either because they usually involve computing the inverse of the Hessian matrix, which is of dimension $pd\times pd$, during each iteration. In this paper, we propose an estimator that is both computationally attractive and asymptotically efficient.
In this section, we present the multinomial-Poisson (MP) transformation, based on which we develop our estimator. We reinterpret the Poisson likelihood as a quasi-likelihood conditional on the covariates.
Let $\boldsymbol{1}_{d}\equiv\left(1,\ldots,1\right)^{\prime}\in\mathbb{R}^{d}$. With a slight abuse of notation, define
where $\boldsymbol{\mu}\equiv\left(\mu_{1},\ldots,\mu_{n}\right)\in\mathbb{R}^{n}$. The following lemma shows that the two functions $l_{C\mid V,M}\left(\boldsymbol{\theta}\right)$ and $l_{C\mid V,M}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$ are the same for any $\boldsymbol{\theta}\in\Theta$ and $\boldsymbol{\mu}\in\mathbb{R}^{n}$. In other words, argument $\boldsymbol{\mu}$ in $l_{C\mid V,M}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$ does not affect the value of the function. The proof of the lemma is straightforward by realizing that $\boldsymbol{C}_{i}^{\prime}\boldsymbol{1}_{d}=M_{i}$ by definition.
Define the following two functions:
It is not difficult to see that $f\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$ takes the form of a log-likelihood function of $n$ conditional Poisson distributions with means $\sum_{k=1}^{d}e^{\left(\eta_{ik}+\mu_{i}\right)}$, $i=1,\ldots,n$ (after ignoring terms that are independent of $\boldsymbol{\theta}$ and $\boldsymbol{\mu}$). This in turn renders $ql_{C\mid V}$ a conditional quasi-log-likelihood function of which $C_{ik}$ given $V_{i}$ is drawn independently from a Poisson distribution with mean $e^{\left(\eta_{ik}+\mu_{i}\right)}$, $k=1,...,d$. This property underlies the naming of the MP transformation.
Based on the conditional quasi-log-likelihood function $ql_{C\mid V}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$, we can compute a conditional quasi MLE (QMLE) of $\boldsymbol{\theta}^{\ast}$:
baker1994multinomial shows that $\widehat{\boldsymbol{\theta}}=\widetilde{\boldsymbol{\theta}}$ if the covariate vector $\boldsymbol{V}$ contains only categorical random variables. The following lemma demonstrates that $\widehat{\boldsymbol{\theta}}=\widetilde{\boldsymbol{\theta}}$ holds irrespective of the type of the covariate vector, thereby generalizing the result of baker1994multinomial.
It is important to note that Lemma (ref) does not depend on the assumption that $f\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$ is the correct log-likelihood function of $M_{i}$, or equivalently that $\Pr\left(M_{i}\mid\boldsymbol{V}_{i}\right)=\mathrm{Po}\left(\sum_{k=1}^{d}e^{\left(\eta_{ik}+\mu_{i}\right)}\right)$ or $\Pr\left(C_{ik}\mid\boldsymbol{V}_{i}\right)=\mathrm{Po}\left(e^{\left(\eta_{ik}+\mu_{i}\right)}\right)$, where $\mathrm{Po}\left(\cdot\right)$ denotes the Poisson distribution. No assumption about the conditional distribution of $M_{i}$ given $\boldsymbol{V}_{i}$ is required for any of the results in the paper to hold. As we show in the following sections, the introduction of $f\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$ is merely a trick to achieve distributed computing.
Lemma (ref) shows that, instead of minimizing $Q_{n}^{\ast}\left(\boldsymbol{\theta}\right)$, we can minimize $Q_{n}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$ to obtain the QMLE of $\boldsymbol{\theta}^{\ast}$. However, computing $\widehat{\boldsymbol{\theta}}$ remains computationally intensive, as solving ((ref)) is impractical when $d$ is large. Nevertheless, the additive structure of $Q_{n}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$ enables the problem to be solved distributively.
As noted in taddy2015distributed, although minimizing $Q_{n}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$ with respect to $\boldsymbol{\theta}$ and $\boldsymbol{\mu}$ jointly is computationally infeasible, given any value of $\boldsymbol{\mu}$, solving $\arg\min_{\boldsymbol{\theta}}Q_{n}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$ is much easier because the optimization can be done separately for each $\boldsymbol{\theta}_{k}$ and be computed across machines. To see this, we rewrite $Q_{n}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$ as:
where for $k=1,\ldots,d$, $Q_{kn}\left(\boldsymbol{\theta}_{k},\boldsymbol{\mu}\right)\equiv\sum_{i=1}^{n}\left(e^{\left(V_{i}^{\prime}\boldsymbol{\theta}_{k}+\mu_{i}\right)}-C_{ik}\left(V_{i}^{\prime}\boldsymbol{\theta}_{k}+\mu_{i}\right)\right)$. In consequence, it holds that for any $\boldsymbol{\mu}$,
Based on ((ref)), solving $\arg\min_{\boldsymbol{\theta}\in\Theta}Q_{n}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$ for any given $\boldsymbol{\mu}$ is equivalent to solving $d$ optimizations: $\arg\min_{\boldsymbol{\theta}_{k}\in\varTheta_{k}}Q_{kn}\left(\boldsymbol{\theta}_{k},\boldsymbol{\mu}\right)$ for each $k=1,\ldots,d$, where each optimization is a Poisson regression.\footnote{Since $\varTheta_{d}=\left\{ \mathbf{0}\right\} $, solving for $\arg\min_{\boldsymbol{\theta}_{d}\in\varTheta_{d}}Q_{dn}\left(\boldsymbol{\theta}_{d},\boldsymbol{\mu}\right)$ is trivial.} Since $\boldsymbol{\theta}_{k}$ has only $p$ dimensions, $\arg\min_{\boldsymbol{\theta}_{k}\in\varTheta_{k}}Q_{kn}\left(\boldsymbol{\theta}_{k},\boldsymbol{\mu}\right)$ is easy to compute. In addition, the optimizations for $k=1,\ldots,d$ can be computed across machines allowing for distributed computing.
By Equation ((ref)), it is not difficult to see that $\widehat{\boldsymbol{\theta}}=\arg\min_{\boldsymbol{\theta}\in\Theta}Q_{n}\left(\boldsymbol{\theta},\widehat{\boldsymbol{\mu}}\right)$. Because $\widehat{\boldsymbol{\theta}}$ is equivalent to the MLE $\widetilde{\boldsymbol{\theta}}$ by Lemma (ref), it has the desired properties such as being both consistent and asymptotically efficient. As a result, we would hope to obtain $\widehat{\boldsymbol{\mu}}$ first and then compute $\widehat{\boldsymbol{\theta}}$ by distributed computing. However, the value of $\widehat{\boldsymbol{\mu}}$ depends on $\widehat{\boldsymbol{\theta}}$, which itself is difficult to calculate. On the other hand, given any value of $\boldsymbol{\theta}$, solving $\arg\min_{\boldsymbol{\mu}}Q_{n}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$ is also fast, and the solution even has a closed form. Denote the solution to $\arg\min_{\boldsymbol{\mu}}Q_{n}\left(\boldsymbol{\theta},\boldsymbol{\mu}\right)$ as $\overline{\boldsymbol{\mu}}_{n}\left(\boldsymbol{\theta}\right)$. A simple calculation would show that
Let $\widehat{\boldsymbol{\mu}}_{T}=\left(\log\left(M_{1}\right),\ldots,\log\left(M_{n}\right)\right)^{\prime}$. Instead of solving for $\widehat{\boldsymbol{\mu}}$ using ((ref)), taddy2015distributed proposes an estimator $\widehat{\boldsymbol{\theta}}_{T}=\arg\min_{\boldsymbol{\theta}\in\Theta}Q_{n}\left(\boldsymbol{\theta},\widehat{\boldsymbol{\mu}}_{T}\right)$ and calls it the distributed computing estimator. Such an estimator is fast to compute but fails to be consistent except in the special cases discussed in taddy2015distributed.
We propose an iterative distributed computing (IDC) estimator, such that during each iteration we solve ((ref)) with $\boldsymbol{\mu}$ updated from the previous step estimate of $\boldsymbol{\theta}$ via ((ref)). Our IDC estimator is defined by the following steps.
Step $\mathbf{0}$. Compute an initial estimator of $\boldsymbol{\theta}^{\ast}$, denoted as $\widehat{\boldsymbol{\theta}}^{\left(0\right)}$.
Step $\mathbf{1},\ldots,\mathbf{S}$. For step $s$, where $s=1,\ldots,S$, we first update $\boldsymbol{\mu}$ using estimator $\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}$ from the previous step via $\overline{\boldsymbol{\mu}}_{n}\left(\cdot\right)$. Then we update $\boldsymbol{\theta}$ given the value of $\boldsymbol{\mu}$:
The iterative estimator with $S$ iterations is defined as $\widehat{\boldsymbol{\theta}}^{I}\equiv\widehat{\boldsymbol{\theta}}^{\left(S\right)}$. For any $\boldsymbol{\theta}$, the value of $\overline{\boldsymbol{\mu}}_{n}\left(\boldsymbol{\theta}\right)$ can be directly computed from ((ref)). In each step, we compute $\arg\min_{\boldsymbol{\theta}_{k}}Q_{kn}\left(\boldsymbol{\theta}_{k},\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}\right)\right)$ for $k=1,\ldots,d$ on $d$ parallel computers. This amounts to running $d$ Poisson regressions with $p$ parameters, and the computational burden for each step is low. The algorithm is described in Algorithm (ref).
Unlike many existing algorithms that numerically solve for the MLE, such as gradient descent or stochastic gradient descent, the IDC estimator does not inherently involve any tuning parameter. This is advantageous because the performance of the classical (stochastic) gradient descent is generally sensitive to the learning rate.
Similar to all iterative optimization procedures, the initial estimator plays a critical role. In finite samples, a good initial guess of $\boldsymbol{\theta}^{\ast}$ can improve the performance of the IDC estimator. Asymptotically, a consistent initial estimator can lead to consistent and asymptotically efficient iterative estimators under weaker assumptions than an inconsistent initial estimator. In this section, we propose three initial estimators: a consistent initial estimator of $\boldsymbol{\theta}^{\ast}$ based on binomial MLE, taddy2015distributed's estimator, and the MLE based on the Poisson assumption of $M_{i}$. The latter two are inconsistent without any assumption on the distribution of $M_{i}$.
\paragraph{A Consistent Initial Estimator }
Let $N_{ik}\equiv C_{ik}+C_{id}$. The following lemma results from the MNL model defined in ((ref)).
Lemma (ref) shows that we can consistently estimate $\boldsymbol{\theta}_{k}^{\ast}$ based on a binomial logistic regression with the log-likelihood function given by:
Let $\check{\boldsymbol{\theta}}_{k}=\arg\min_{\boldsymbol{\theta}_{k}\in\varTheta_{k}}-l_{C_{k},C_{d}\mid V,N_{k}}\left(\boldsymbol{\theta}_{k}\right)$ for $k=1,\ldots,d-1$ and $\check{\boldsymbol{\theta}}=\left(\boldsymbol{\check{\theta}}_{1}^{\prime},\ldots,\boldsymbol{\check{\theta}}_{d-1}^{\prime},\boldsymbol{\check{\theta}}_{d}^{\prime}\right)^{\prime}$ with $\boldsymbol{\check{\theta}}_{d}^{\prime}=\mathbf{0}$. The consistency of $\check{\boldsymbol{\theta}}$ follows from standard arguments in the maximum likelihood estimation.
Compared to $\widetilde{\boldsymbol{\theta}}$, the conditional probability used in constructing the above binomial logistic log-likelihood function does not use all the available information. Therefore, $\check{\boldsymbol{\theta}}$ is less efficient than $\widetilde{\boldsymbol{\theta}}$. However, each component of $\check{\boldsymbol{\theta}}$, $\check{\boldsymbol{\theta}}_{k}$, can be calculated independently, allowing for parallel computing. The substantially short running time of $\check{\boldsymbol{\theta}}$ makes it a great candidate for the initial value $\widehat{\boldsymbol{\theta}}^{\left(0\right)}$.
\paragraph{Inconsistent Initial Estimators }
Even though taddy2015distributed's estimator fails to be consistent in general, it could serve as a candidate for the initial value in our Algorithm (ref). Another option is to replace $\widehat{\boldsymbol{\mu}}_{T}$ with a zero vector to obtain another estimator denoted as $\widehat{\boldsymbol{\theta}}_{P}=\arg\min_{\boldsymbol{\theta}\in\Theta}Q_{n}\left(\boldsymbol{\theta},\mathbf{0}\right)$. It can also be computed across machines for each $k=1,...,d$. Like taddy2015distributed's estimator, $\widehat{\boldsymbol{\theta}}_{P}$ could also serve as a candidate for the initial value in our algorithm. Moreover, under an extra condition that $M_{i}$ follows a Poisson distribution, $\widehat{\boldsymbol{\theta}}_{P}$ is the maximum likelihood estimator of $\boldsymbol{\theta}^{\ast}$.
Unlike $\check{\boldsymbol{\theta}}$, neither $\widehat{\boldsymbol{\theta}}_{T}$ nor $\widehat{\boldsymbol{\theta}}_{P}$ is a consistent estimator of $\boldsymbol{\theta}^{\ast}$ without any additional assumption.
In some applications, researchers may have prior knowledge of some linear equality constraints among parameters. Taking the constraints into consideration during the estimation would further improve asymptotic efficiency. In this section, we discuss how to modify our IDC estimator introduced in Section (ref) to incorporate equality constraints. The initial estimators introduced in Section (ref), although they do not account for the equality constraints, can be utilized to obtain the initial $\widehat{\boldsymbol{\theta}}^{\left(0\right)}$.
We consider two different types of constraints: constraints on parameters for the same choice and different choices. The procedures for different types of constraints differ in the optimization problems during each iteration. For each type, we use an example to illustrate our procedure.
For the first type, the constraint is on components of individual $\boldsymbol{\theta}_{k}^{\ast}$. Take the constraint $\theta_{k1}^{\ast}=\theta_{k2}^{\ast}$ for $k=1,\ldots d$ as an example. When computing $\widehat{\boldsymbol{\theta}}^{\left(s\right)}$ in Step $s$, we solve the constrained optimization problem: \[ \arg\min_{\boldsymbol{\theta}_{k}\in\varTheta_{k},\theta_{k1}=\theta_{k2}}Q_{kn}\left(\boldsymbol{\theta}_{k},\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}\right)\right) \] for each $k$. Because the constraint is on each $\boldsymbol{\theta}_{k}^{\ast}$, the original distributed computing scheme remains.
The second type of constraints involves components of $\boldsymbol{\theta}_{k}^{\ast}$ across different choices. For example, researchers may impose restrictions like $\theta_{11}^{\ast}=\cdots=\theta_{q1}^{\ast}$, where $q<d$. For $k=1,\ldots,q$, let $\boldsymbol{\theta}_{k,-1}$ be the subvector of $\boldsymbol{\theta}_{k}$ that excludes its first element. From Steps $1$ to $S$, we first update $\boldsymbol{\mu}$ using estimator $\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}$ from the previous step. Then, we compute $\widehat{\boldsymbol{\theta}}^{\left(s\right)}$ from the following optimization problems:
Optimization problems ((ref)) can be solved using parallel computers. And ((ref)) is an optimization problem with only one argument.\footnote{Let a general linear equality constraint of the second type be written as $R\boldsymbol{\theta}^{\ast}=r$, where $R$ and $r$ are known with dimensions $l_{R}\times pd$ and $l_{R}\times1$ respectively. The matrix $R$ is assumed to have full row rank so that there is no redundant constraint. We can always rearrange and decompose $R$ as $\left[R_{c},\mathbf{0}\right]$, where $R_{c}$ has dimension $l_{R}\times q$ and has no zero column. The number of arguments in the optimization problem ((ref)) for $R\boldsymbol{\theta}^{\ast}=r$ is $q-l_{R}\geq0$. The case where $q=l_{R}$ corresponds all $q$ number of elements in $\boldsymbol{\theta}^{\ast}$ having prespecified values. The optimization problem ((ref)) is no longer needed in this case.} In consequence, each step incurs a low computational burden. The IDC estimator with such a constraint is also fast to compute.
The aforementioned two procedures can be generalized to accommodate any linear equality constraint. In particular, the two procedures can be combined in a straightforward way when constraints contain both types.
In this section, we establish the consistency and asymptotic normality of our IDC estimator introduced in Section (ref).\footnote{Asymptotic properties of the constrained estimators discussed in Section (ref) can be established in the same way. Moreover, the procedures can be used to construct a likelihood ratio test for testing the null hypothesis on linear equality constraint among parameters.} Technical proofs are collected in Appendix (ref). We first impose the following two assumptions.
Assumption (ref) is standard. Assumption (ref) (ii) requires that the moment generating function of $\boldsymbol{V}$ exists within $\Theta$. Almost all commonly seen distributions satisfy Assumption (ref) (ii) and (iii).
In order to analyze the asymptotic properties of the IDC estimator, we need to study the way in which the iterative estimator updates itself in each step. In Step $s$ of the iteration, we first compute $\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}\right)$ using $\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}$ from the previous step and then calculate $\arg\min_{\boldsymbol{\theta}\in\Theta}Q_{n}\left(\boldsymbol{\theta},\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}\right)\right)$. To explicitly distinguish the argument in $\overline{\boldsymbol{\mu}}_{n}\left(\cdot\right)$ from the argument in $Q_{n}\left(\cdot,\boldsymbol{\mu}\right)$ for any given $\boldsymbol{\mu}$, we introduce $\boldsymbol{\vartheta}$ and use it as the argument in $\overline{\boldsymbol{\mu}}_{n}\left(\cdot\right)$. Define $Q_{n}^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)\equiv Q_{n}\left(\boldsymbol{\theta},\overline{\boldsymbol{\mu}}_{n}\left(\boldsymbol{\vartheta}\right)\right)$. We have that \[ Q_{n}^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)\equiv\sum_{i=1}^{n}\sum_{k=1}^{d}\left(\frac{M_{i}e^{\boldsymbol{V}_{i}^{\prime}\boldsymbol{\theta}_{k}}}{\sum_{k=1}^{d}e^{\boldsymbol{V}_{i}^{\prime}\boldsymbol{\vartheta}_{k}}}-C_{ik}\boldsymbol{V}_{i}^{\prime}\boldsymbol{\theta}_{k}-C_{ik}\log\left(\frac{M_{i}}{\sum_{k=1}^{d}e^{\boldsymbol{V}_{i}^{\prime}\boldsymbol{\vartheta}_{k}}}\right)\right). \] Further define function $Q^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)$ as the probability limit of $\frac{1}{n}Q_{n}^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)$:
In the following lemma, we provide some properties of $Q^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)$, which are crucial for obtaining the consistency and asymptotic normality of our IDC estimator.
Essentially, $\overline{\boldsymbol{\theta}}\left(\boldsymbol{\vartheta}\right)$ summarizes the operation in each step with $\boldsymbol{\vartheta}$ being the input and $\overline{\boldsymbol{\theta}}\left(\boldsymbol{\vartheta}\right)$ being the output when the sample size goes to infinity. By part (ii) of Lemma (ref), $\overline{\boldsymbol{\theta}}\left(\cdot\right)$ is well-defined. Part (v) of Lemma (ref) plays the most important role. Heuristically, for any given $\boldsymbol{\vartheta}$, the value of function $\overline{\boldsymbol{\theta}}\left(\boldsymbol{\vartheta}\right)$ is obtained by solving $\frac{\partial Q^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)}{\partial\boldsymbol{\theta}}=\mathbf{0}$. At the same time, function $\frac{\partial Q^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)}{\partial\boldsymbol{\theta}}$ relates to the first order derivative of $L_{C\mid V,M}\left(\boldsymbol{\theta}\right)$, the population objective function defined in Section (ref). By the identification assumption and the convexity of $-L_{C\mid V,M}\left(\boldsymbol{\theta}\right)$, only the true value $\boldsymbol{\theta}^{\ast}$ satisfies that $\frac{\partial}{\partial\boldsymbol{\theta}}L_{C\mid V,M}\left(\boldsymbol{\theta}\right)=\mathbf{0}$, which implies that $\frac{\partial Q^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)}{\partial\boldsymbol{\theta}}\mid_{\boldsymbol{\vartheta}=\boldsymbol{\theta}}=\mathbf{0}$ holds only at $\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)=\left(\boldsymbol{\theta}^{\ast},\boldsymbol{\theta}^{\ast}\right)$.
If a consistent initial estimator is used, such as $\check{\boldsymbol{\theta}}$, then Lemma (ref) is sufficient for the consistency of the IDC estimator as stated below.
On the other hand, if the initial estimator is consistent only under extra assumptions, such as $\widehat{\boldsymbol{\theta}}_{T}$ and $\widehat{\boldsymbol{\theta}}_{P}$, or even inconsistent, then we need a contraction mapping assumption on $\overline{\boldsymbol{\theta}}\left(\cdot\right)$.
Admittedly, Assumption (ref) is a high-level assumption. Based on the evidence from the simulation, the assumption holds for various values of $\boldsymbol{\theta}^{\ast}$ and distributions of $\boldsymbol{V}$ and $M$. Assumption (ref) relates to the contraction mapping assumption (Assumption 6) in pastorello2003iterative but is weaker, the reason being that the true $\boldsymbol{\theta}^{\ast}$ is the unique fixed point, see Lemma (ref) (v). Specifically, Assumption (ref) only requires that the distance between $\overline{\boldsymbol{\theta}}\left(\boldsymbol{\vartheta}\right)$ and $\boldsymbol{\vartheta}$ get smaller after both being mapped by $\overline{\boldsymbol{\theta}}\left(\cdot\right)$. Instead, the contraction mapping assumption (Assumption 6) in pastorello2003iterative requires that the distance between two arbitrary $\boldsymbol{\vartheta}^{1}$ and $\boldsymbol{\vartheta}^{2}$ get smaller after both being mapped by $\overline{\boldsymbol{\theta}}\left(\cdot\right)$.
Under Assumptions (ref) and (ref), the MLE $\widetilde{\boldsymbol{\theta}}$ is asymptotically normally distributed with asymptotic variance given by the Fisher information matrix \[ \mathcal{I}\left(\boldsymbol{\theta}^{\ast}\right)\equiv\frac{\partial^{2}}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{\prime}}Q^{\ast}\left(\boldsymbol{\theta}^{\ast}\right), \] where $Q^{\ast}\left(\boldsymbol{\theta}^{\ast}\right)\equiv p\lim_{n\rightarrow\infty}\frac{1}{n}Q_{n}^{\ast}\left(\boldsymbol{\theta}^{\ast}\right)$. In this section, we show that our IDC estimator has the same asymptotic distribution as $\widetilde{\boldsymbol{\theta}}$, based on which we introduce a valid bootstrap inference procedure.
The conditions required for proving the asymptotic distribution result depend on the initial estimator $\widehat{\boldsymbol{\theta}}^{\left(0\right)}$. If a consistent initial estimator is used, then the following assumption is sufficient. For any matrix $A$, denote $\left\Vert A\right\Vert $ as its spectral norm.
The detailed expressions of $\frac{\partial^{2}Q^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{\prime}}$ and $\frac{\partial^{2}Q^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)}{\partial\boldsymbol{\theta}\partial\boldsymbol{\vartheta}^{\prime}}$ can be found in Appendix (ref). Assumption (ref) is often called the information dominance condition and is tantamount to the local contraction mapping condition. It is weaker than Assumption (ref). Because we have an initial consistent estimator of $\boldsymbol{\theta}^{\ast}$, Assumption (ref) can be verified. Additionally, because the matrix $\frac{\partial^{2}Q^{\dagger}\left(\boldsymbol{\theta},\boldsymbol{\vartheta}\right)}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{\prime}}$ is block diagonal with each block having dimensions $p\times p$, computing its inverse is feasible.
The following theorem shows that when $S$ is sufficiently large, the IDC estimator $\widehat{\boldsymbol{\theta}}^{\left(S\right)}$ is equal to $\widetilde{\boldsymbol{\theta}}$ up to a term of order smaller than $n^{-1/2}$.
Theorem (ref) shows that we do not lose efficiency when employing the proposed IDC estimator as long as $S$ is large enough. A direct implication of the theorem is that $\widehat{\boldsymbol{\theta}}^{\left(S\right)}$ has the same asymptotic distribution as $\widetilde{\boldsymbol{\theta}}$ for sufficiently large $S$.
To conduct inference on $\boldsymbol{\theta}^{\ast}$ based upon $\widehat{\boldsymbol{\theta}}^{\left(S\right)}$, we need to consistently estimate the Fisher information matrix and compute its inverse. Because the dimension of $\mathcal{I}\left(\boldsymbol{\theta}^{\ast}\right)$ is $dp\times dp$, calculating the inverse of its estimator is not only time-consuming but also unreliable when $d$ is large. As a result, we proceed by applying the following parametric bootstrap, which is feasible thanks to the fact that the IDC estimator is fast to compute.
Given $\left\{ \left(\boldsymbol{V}_{i},M_{i}\right)\right\} _{i=1}^{n}$ of the original sample, we draw the bootstrap sample $\boldsymbol{C}_{in}^{\ddagger}$ for $i=1,\ldots,n$ from the multinomial logistic regression model with the conditional probability mass function given by \[ \mathrm{MNL}\left(\boldsymbol{C}_{in}^{\ddagger};\widehat{\boldsymbol{\eta}}_{i},M_{i}\right)\textrm{, where }\widehat{\eta}_{ik}=\boldsymbol{V}_{i}^{\prime}\widehat{\boldsymbol{\theta}}_{k}^{\left(S\right)}\textrm{ for }k=1,\ldots d. \] The bootstrap version of the iterative estimator $\widehat{\boldsymbol{\theta}}^{\ddagger\left(S\right)}$ is obtained by applying the algorithm introduced in Section (ref) with bootstrap sample $\left\{ \left(\boldsymbol{C}_{in}^{\ddagger},\boldsymbol{V}_{i},M_{i}\right)\right\} _{i=1}^{n}$.
Assume that we start with a consistent initial estimator. Based on Theorem (ref) (i), we have that for $S\geq\log\left(n\right)$,
Define the score function for $\boldsymbol{\theta}$ as \[ \dot{\mathbf{l}}\left(\boldsymbol{\theta}\mid\boldsymbol{c},\boldsymbol{v},m\right)\equiv\frac{d}{d\boldsymbol{\theta}}\log\mathrm{MNL}\left(\boldsymbol{c};\boldsymbol{\eta},m\right)\textrm{, where }\eta_{k}\equiv\boldsymbol{v}^{\prime}\boldsymbol{\theta}_{k}. \] The iterative estimator $\widehat{\boldsymbol{\theta}}^{\left(S\right)}$ is asymptotically linear with influence function $\widetilde{\mathbf{l}}\left(\boldsymbol{\theta}^{\ast}\mid\boldsymbol{c},\boldsymbol{v},m\right)$: \[ \sqrt{n}\left(\widehat{\boldsymbol{\theta}}^{\left(S\right)}-\boldsymbol{\theta}^{\ast}\right)=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\widetilde{\mathbf{l}}\left(\boldsymbol{\theta}^{\ast}\mid\boldsymbol{C}_{i},\boldsymbol{V}_{i},M_{i}\right)+o_{p}\left(1\right), \] where $\widetilde{\mathbf{l}}\left(\boldsymbol{\theta}^{\ast}\mid\boldsymbol{c},\boldsymbol{v},m\right)\equiv\mathcal{I}^{-1}\left(\boldsymbol{\theta}^{\ast}\right)\dot{\mathbf{l}}\left(\boldsymbol{\theta}^{\ast}\mid\boldsymbol{c},\boldsymbol{v},m\right)$. Applying the same derivation, we can show that the bootstrap version of the estimator is also asymptotically linear with the influence function evaluated at $\widehat{\boldsymbol{\theta}}^{\left(S\right)}$: \[ \sqrt{n}\left(\widehat{\boldsymbol{\theta}}^{\ddagger\left(S\right)}-\widehat{\boldsymbol{\theta}}^{\left(S\right)}\right)=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\widetilde{\mathbf{l}}\left(\widehat{\boldsymbol{\theta}}^{\left(S\right)}\mid\boldsymbol{C}_{in}^{\ddagger},\boldsymbol{V}_{i},M_{i}\right)+o_{p}\left(1\right). \] The Lindeberg-Feller central limit theorem proves the bootstrap consistency. The proof for the case of an inconsistent initial estimator is analogous. Let $\stackrel{d^{\ddagger}}{\rightarrow}$ denote the convergence in bootstrap distribution.
In this section, we evaluate the performance of our IDC estimator from various perspectives. We present the finite sample performance of the IDC estimator by looking at the effects of separately increasing $d$ and $n$ on mean squared error (MSE) and running time. We include the maximum likelihood estimator to show that the IDC estimator performs similarly to the MLE in terms of MSE and is always feasible even in cases where MLE is intractable. Lastly, we study the finite sample size and power of our bootstrap inference procedure.
In what follows, we present results on the finite sample performance of the IDC estimator in terms of MSE and running time in four tables. The reported running times in the tables are obtained from a cluster of 25 AWS EC2 instances with 12 vCPUs and 16GB memory. Such a configuration can be formed on commonly used cloud computing platforms within minutes. We employ this basic configuration to illustrate that our IDC estimator achieves superior performance compared to existing estimators, even when computational resources are suboptimal for distributed computing. All of the results presented in this section are repeated five hundred times and averaged.
We first study the finite sample performance of the IDC estimator with consistent $\check{\boldsymbol{\theta}}$ initialization, specifically what happens to MSE ($n$ increasing, $d$ fixed) and running time ($d$ increasing, $n$ fixed). We consider the following data generating process (DGP):
DGP-A {[}MNL{]}: We set $p=5$. The covariate vector $\boldsymbol{V}$ follows the standard normal distribution; $M$ follows the discrete uniform distribution on $\left[20,30\right]$; and the values of $\mathbf{\boldsymbol{\theta}}^{*}$ are obtained by random draws from the standard normal distribution.
From Table (ref), we observe that the MSE of the IDC estimator decreases as the sample size increases. When the number of iterations $S$ increases, we see an improvement in the MSE. However, the improvement is marginal, suggesting that the IDC estimator with $\check{\boldsymbol{\theta}}$ initialization stabilizes with only a few iterations. We also see from Table (ref) that the running time is approximately a linear function of $d$. When the number of cores available does not exceed the number of choices, the additional computational cost of increasing $d$ is very small. After the cores are fully occupied by the number of processes, the running time becomes approximately linear. The nominal value of running time depends on the hardware specifications.
To compare the performance of our IDC estimator with the MLE $\widetilde{\boldsymbol{\theta}}$, we simulate the MSE and running time of $\widetilde{\boldsymbol{\theta}}$ from the same DGP and present the result in Table (ref).\footnote{We also write code to try estimators in bohning1988monotonicity, bohning1992multinomial, and simon2013blockwise for comparison. However, our simulation result suggests that their performance depends crucially on the number of iterations.} It can be seen that the MSE of the IDC estimator with $\check{\boldsymbol{\theta}}$ initialization is very close to that of the MLE even when the number of iterations is only $10$. Note that the main advantage of the parallel estimator is best observed for high enough $d$ because for low $d$, the communication between parallel processes is unnecessary and hence parallel computing increases the running time unnecessarily. The superior performance of the IDC estimator is apparent when $d$ is large. For instance, when $d=150$ and $n=2000$, the IDC estimator with $S=10$ achieves a similar MSE as the MLE with only about one-tenth of the running time. For higher-dimensional cases, such as when $d$ exceeds $150$, computing the MLE becomes computationally intensive and may not be practical for many applications. In comparison, the IDC estimator with $S=10$ demonstrates more efficient computation times, requiring approximately 5, 10, and 20 minutes for sample sizes $n=500$, $1000$, and $2000$, respectively. Moreover, the running time for the IDC estimator can be further decreased if more compute instances are used. For example, using 96 instances, the running time of the IDC estimator for $d=150$ can be further reduced to $34$, $51$, and $188$ seconds for $n=500$, $1000$, $2000$ respectively even for $S=40$. Compared to the corresponding running time of MLE, the running time of the IDC estimator using 96 instances is more than 10, 30, and 50 times shorter.
In Table (ref), we present the MSE and running time of IDC estimators with $\widehat{\boldsymbol{\theta}}_{T}$ and $\widehat{\boldsymbol{\theta}}_{P}$ as initial estimators, respectively, for two different numbers of iterations $S$. Comparing MSEs of three IDC estimators with different initial values: $\check{\boldsymbol{\theta}}$ (Table (ref)) and $\widehat{\boldsymbol{\theta}}_{T}$ (Table (ref)) or $\widehat{\boldsymbol{\theta}}_{P}$ (Table (ref)), we observe that the IDC estimator with the consistent initial estimator reduces the MSEs for the same number of iterations.
Table (ref) presents MSEs of different estimators for three $\left(d,n\right)$ pairs each. We set the largest $d$ be $50$ so that MLE can be computed in a reasonable time. $\widehat{\boldsymbol{\theta}}_{PB}^{I}$, $\widehat{\boldsymbol{\theta}}_{T}^{I}$, and $\widehat{\boldsymbol{\theta}}_{P}^{I}$ denote the IDC estimators with $\check{\boldsymbol{\theta}}$, $\widehat{\boldsymbol{\theta}}_{T}$, and $\widehat{\boldsymbol{\theta}}_{P}$ as the initial estimators respectively. Besides DGP-A, we consider two additional DGPs to study the performance of the IDC estimator under different data settings. In all DGPs, we let $p=5$.
DGP-B {[}Poisson{]}: The random variable $\boldsymbol{V}$ follows a standard normal distribution; $C_{ik}$ follows a Poisson distribution with mean $e^{\eta_{ik}^{*}}$; and $M$ is obtained by summing up realizations of the Poisson draws for different choices.
DGP-C {[}Mixture{]}: We let $\boldsymbol{V}$ follow a mixture of Gaussian distributions with means $0$ and $4$ with standard deviations $1$ for both distributions. $M$ is also set to follow a mixture of Gaussian distributions with means $10$ and $60$ and rounded to the closest integer. The standard deviations are $1$ and $5$ respectively. We have made these modifications so that some choices are rarely selected and ensure the robustness of our estimator in those cases.
Based on the simulation result, our IDC algorithm is successfully executed for all DGPs and exhibits stability. In contrast, we encounter errors for DGP-C when computing the estimators in bohning1988monotonicity, bohning1992multinomial, and simon2013blockwise. We see from Table (ref) that $\widehat{\boldsymbol{\theta}}_{PB}^{I}$ performs close to $\widetilde{\boldsymbol{\theta}}$ for all DGPs and $\left(d,n\right)$ pairs. In DGP-B, the initial estimator $\widehat{\boldsymbol{\theta}}_{P}$ is the maximum likelihood estimator. As a result, $\widehat{\boldsymbol{\theta}}_{P}^{I}$ starts with not only a consistent but asymptotically efficient initial estimator. Even in this case, $\widehat{\boldsymbol{\theta}}_{PB}^{I}$ has comparable MSEs.
In summary, the IDC estimators with all three initial estimators have finite sample performance similar to the MLE for the DGPs studied in this section. They are much faster to compute than the MLE for large $d$ and are feasible even when the MLE might be intractable. Moreover, if the IDC estimator starts with the consistent initial estimator $\check{\boldsymbol{\theta}}$, its finite sample performance will be further improved and is almost the same as the MLE.
In this section, we illustrate the bootstrap inference procedure introduced in Section (ref). We investigate the finite sample performance of the procedure including the size and power. All the results are based on one thousand Monte Carlo repetitions, where the number of bootstrap repetitions is five hundred.
We consider the null hypothesis that some element of $\boldsymbol{\theta}^{\ast}$ equals to a specific value. Data are generated from DGP-A introduced in the previous section. Let the null and the alternative hypotheses be that $H_{0}:\theta_{11}^{\ast}=0$ and $H_{1}:\theta_{11}^{\ast}\neq0$. The test statistic is computed as $\left|\frac{\widehat{\theta}_{11}^{I}}{\widehat{se}_{b}\left(\widehat{\theta}_{11}^{I}\right)/\sqrt{n}}\right|$, where $\widehat{se}_{b}\left(\widehat{\theta}_{11}^{I}\right)$ is the bootstrap estimate of the standard error of $\widehat{\theta}_{11}^{I}$. The number of iterations is 10 when computing the IDC estimator. We set the nominal size as $5\%$ and use the $97.5\%$ quantile of the standard normal distribution as the critical value.
In Table (ref), we report the finite sample rejection probabilities of our test for different values of $\theta_{11}^{\ast}$. Values in the first row of the table indicate the deviation of $\theta_{11}^{\ast}$ from the null hypothesis. When the deviation is zero, the null hypothesis is true. It can be seen from the table that the finite sample rejection rates get closer to the nominal size when the sample size increases. And when the true value $\theta_{11}^{\ast}$ deviates more from the null hypothesis, the rejection probabilities increase. The same pattern appears for both $d=20$ and $d=50$. The finite sample performance of the test when $d=50$ is not as good as that when $d=20$. This is predictable because there are many more unknown parameters in the model when $d=50$ than when $d=20$. We expect the results to improve as the sample size increases for any fixed $d$.
In this paper, we propose an iterative distributed computing estimator for the multinomial logistic model that is fast to compute even when the number of choices is large. When the number of iterations goes to infinity, we show that our estimator is both consistent and asymptotically efficient. Based on the simulation study, the computational time of our estimator increases linearly with the number of choices. Moreover, our estimator has comparable finite sample performance to MLE when the latter is computationally feasible.
Extensions abound. First, our IDC estimator can be combined with several existing algorithms to accommodate more complex settings. For example, when minimizing $Q_{kn}\left(\boldsymbol{\theta}_{k},\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}\right)\right)$ for $k=1,\ldots,d$, we can replace the gradient of the objective function with its stochastic approximation calculated from a randomly selected subset of the data. Such an algorithm is an online algorithm and might reduce the running time especially when $n$ is large. We can also employ a one-step Newton-Raphson to compute $\arg\min_{\boldsymbol{\theta}_{k}\in\varTheta_{k}}Q_{kn}\left(\boldsymbol{\theta}_{k},\overline{\boldsymbol{\mu}}_{n}\left(\widehat{\boldsymbol{\theta}}^{\left(s-1\right)}\right)\right)$ in each iteration or even replace the Hessian matrix with its dominant. These modifications to the IDC estimator have the potential to further enhance computational efficiency, depending on the application. However, their theoretical properties require further investigation. Second, the asymptotic properties of the IDC estimator in this paper are established for a large but fixed number of choices. Asymptotic theory allowing for the number of choices to diverge with the sample size is yet to be established. Third, in cases where the number of covariates is also large, $\ell_{1}$ or $\ell_{2}$ regularization could be adopted. In a companion paper, we develop an asymptotic theory for a regularized iterative distributed computing estimator.