On the "Poisson Trick" and its Extensions for Fitting Multinomial Regression Models
Jarod Y. L. Lee Peter J. Green Louise M. Ryan
stat.ME stat.AP stat.CO




1 2]Jarod Y.L. Lee
1 3]Peter J. Green
1 2 4]Louise M. Ryan
[1]School of Mathematical and Physical Sciences, University of Technology Sydney, Australia.
[2]Australian Research Council Centre of Excellence for Mathematical & Statistical Frontiers, The University of Melbourne, Australia.
[3]School of Mathematics, University of Bristol, U.K.
[4]Department of Biostatistics, Harvard T.H. Chan School of Public Health, U.S.

On the “Poisson Trick” and its Extensions for Fitting Multinomial Regression Models
    [
    December 30, 2023
===================================================================================








	This article is concerned with the fitting of multinomial regression models using the so-called “Poisson Trick”. The work is motivated by <cit.> and <cit.> which have been criticized for being computationally inefficient and sometimes producing nonsense results. We first discuss the case of independent data and offer a parsimonious fitting strategy when all covariates are categorical. We then propose a new approach for modelling correlated responses based on an extension of the Gamma-Poisson model, where the likelihood can be expressed in closed-form. The parameters are estimated via an Expectation/Conditional Maximization (ECM) algorithm, which can be implemented using functions for fitting generalized linear models readily available in standard statistical software packages. Compared to existing methods, our approach avoids the need to approximate the intractable integrals and thus the inference is exact with respect to the approximating Gamma-Poisson model. The proposed method is illustrated via a reanalysis of the yogurt data discussed by <cit.>. 
	
2mm
Keywords:  Discrete choice model; Longitudinal data; Mixed logit model; Multinomial mixed model; Nominal polytomous data;  Unobserved heterogeneity.

 






§ INTRODUCTION


Data with correlated categorical responses arise frequently in applications. This may arise from units grouped into clusters (clustered data) or multiple measurements taken on the same unit (longitudinal data). For instance, we might expect the unemployment outcomes (employed, unemployed, not in labour force) of residents living in the same region to be correlated, due to similar job opportunities and socioeconomic levels. Ignoring the correlation structure and assuming that all observations are independent by fitting an ordinary multinomial regression model may result in biased estimates and inaccurate predictions. Multinomial mixed models can account for correlation by using group level random effects <cit.>. 
 

For multinomial mixed models, it is a common practice to assume a multivariate normal distribution for the random effects. The multivariate normal distribution is easy to interpret and is convenient when we want to build more complicated correlation structures into our model. However, the resulting likelihood involves multidimensional integrals that cannot be solved analytically. The computational effort to evaluate the likelihood increases with the number of groups and categories, making it not suitable for large scale applications. In fact, <cit.> showed that closed-form likelihoods for multinomial mixed models do not exist regardless of the random effect distribution, except for the special case when there are no covariates. 



Various methods have been proposed to circumvent the computational obstacle for fitting multinomial mixed models. Among them are quadrature <cit.>, Monte Carlo EM algorithm, pseudo-likelihood approach <cit.> and Markov Chain Monte Carlo methods <cit.>. <cit.> proposed a random effects estimation approach using a discrete probability distribution approximation. Simulation based methods such as method of simulated moments <cit.> and method of simulated maximum likelihood <cit.> are widely used in the  econometrics literature. Recently, <cit.> proposed a fast moment-based estimation method that scales well for large samples and  which arguably can be extended for fitting multinomial mixed models. <cit.> used the fact that the multinomial model is a multivariate binary model and exploited a procedure proposed by <cit.> for model fitting. Their approach has been criticized by <cit.> as they failed to realize that a multivariate link function is needed in the context of multinomial models. An alternative strategy using clustered bootstrap was subsequently proposed by <cit.>. Although some authors have considered the Dirichlet-multinomial model that results in a closed-form likelihood, it does not allow the incorporation of individual level covariates. 






<cit.> advocate using Poisson log-linear or non-linear mixed models, both with random effects, as surrogates to multinomial mixed models. Their method capitalizes on existing mixed models software packages for fitting generalized linear models with random effects. This allows multinomial mixed models to be fitted using an approximate likelihood from the Poisson surrogate models. Their results are based on extensions of the well known “Poisson Trick” <cit.> that relates multinomial models with Poisson models via a respecification of the model formulae. Although clever, their methods have been criticized for being computationally inefficient <cit.> and sometimes producing nonsense results <cit.>. This might be due to the intractable likelihoods of their models and the various approximation methods being used in different software packages. The considerable execution time is especially problematic, where it can take up to months to fit the model to a moderate-sized  dataset <cit.>! In this article, we propose a new approach based on an extension of the Gamma-Poisson model <cit.>, where the likelihood can be expressed in closed-form. Using the proposed estimation procedure, the parameters can be estimated via readily available packages for fitting generalized linear models. 


The remaining paper is organized as follows. Section <ref> reviews the “Poisson Trick” for multinomial regression models with independent responses, and suggests a parsimonious fitting strategy when all covariates are categorical. In Section <ref> we propose a new approach for approximating the likelihood of multinomial regression models with random effects. The empirical performance of the proposed model is demonstrated via a simulation study and a reanalysis of the yogurt brand choice dataset as discussed by <cit.> in Section <ref>. Finally, we conclude with a discussion and a summary of our findings in Section <ref>.












§ "POISSON TRICK" FOR INDEPENDENT MULTINOMIAL RESPONSES
 

This section describes the relationship between multinomial and Poisson regression models for independent responses, which we refer to as the “Poisson Trick”. The results are based on the well known fact that given the sum, Poisson counts are jointly multinomially distributed <cit.>. 



 §.§ Derivation
 
Let Y_j = (Y_jq)_q=1^Q be the Q × 1 response vector for observation j with the corresponding probability vector p_j = (p_jq)_q=1^Q, where q indexes the multinomial category. A common approach to satisfy the two characteristics of probability (i) 0 ≤ p_jq≤ 1 for all j and q; and (ii) ∑_q=1^Q p_jq = 1 is via
    p_jq = ζ_jq/ζ_j+,
 
where ζ_jq is a positive user-specified function of covariates x and fixed effects γ, and ζ_j+ = ∑_q=1^Q ζ_jq. Depending on the type of variable under consideration, x and γ can be indexed by various combinations of j and q <cit.>. Conditional on the multinomial sums Y_j+ = ∑_q=1^Q Y_jq, Y_js are independently multinomially distributed for each j, i.e.

    Y_j|Y_j+∼ℳ(Y_j+,p_j) .
 

In multinomial models, Y_j+ = y_j+ is treated as fixed. Suppose we instead treat Y_j+ as random and assume 

    Y_j+∼𝒫(δ_j ζ_j+),

independently for each j. This results in a multinomial-Poisson mixture with the following joint probability function for each j:

    P(Y_j = y_j ∩ Y_j+=y_j+)    =  P(Y_j+=y_j+) P(Y_j = y_j| Y_j+=y_j+) 
       = e^-δ_jζ_j+(δ_jζ_j+)^y_j+/y_j+!×y_j+!/∏_q y_jq!∏_q (ζ_jq/ζ_j+)^y_jq
       = ∏_q {e^-δ_j ζ_jq(δ_j  ζ_jq)^y_jq/y_jq!}  iff   Y_j+ = y_j+.
 

The marginal probability of Y_j can then be obtained by summing the joint probability over all possible values of Y_j+:


    P(Y_j=y_j)    = ∑_Y_j+=0^∞∏_q {e^-δ_j ζ_jq (δ_j ζ_jq)^y_jq/y_jq!}
       = ∏_q {e^-δ_j ζ_jq (δ_j  ζ_jq)^y_jq/y_jq!}.
  
Thus, allowing the multinomial sums to be random according to a Poisson distribution results in

    Y_jq∼𝒫(δ_j ζ_jq),

independently for each j and q. Summing over all observations, the log-likelihood is

    ∑_j ℓ^𝒫(δ_jζ_j+;Y_j+) + ∑_j ℓ^ℳ(ζ_j;Y_j|Y_j+) = ∑_j ∑_q ℓ^𝒫(δ_jζ_jq; Y_jq),

where ℓ^𝒫 and ℓ^ℳ denote the Poisson and multinomial log-likelihood functions respectively, and ζ_j = (ζ_jq)_q=1^Q. The second term on the left hand side is the model we would like to fit, and the term on the right hand side is the model we actually fit. 


To show that the Poisson surrogate model is an exact fit to the multinomial model, first note the log-likelihood corresponding to the multinomial model is

    ∑_j log(y_j+!) - ∑_j∑_q log(y_jq!) + ∑_j∑_q y_jqlogζ_jq - ∑_j y_j+logζ_j+,

and the log-likelihood of the Poisson surrogate model is

    - ∑_jδ_j ζ_j+ + ∑_j y_j+logδ_j + ∑_j∑_q y_jqlogζ_jq - ∑_j log(y_j+!).
 

Differentiating Equation <ref> with respect to δ_j and setting it to 0, we obtain δ̂_j = y_j+/ζ_j+. Plugging in the maximizing value of δ_j into Equation <ref>, we have

    -∑_j y_j+ + ∑_j y_j+log y_j+ - ∑_j y_j+logζ_j+ + ∑_j∑_q y_jqlogζ_jq.
 

Equation <ref> is identical to Equation <ref>, up to an additive constant. It follows that the maximum likelihood estimates, their asymptotic variances and tests for the fixed effects can be exactly recovered under the Poisson surrogate model <cit.>. That is, likelihood inference for ζ_jq is the same whether we regard Y_j+ as fixed (multinomial) or randomly sampled from independent Poissons. This result applies to the fixed effects model, with any parameterization of ζ_jq, including:

	
  * Exponential transformations of linear combinations of categorical variables and regression coefficients <cit.>,
	
  * Exponential transformations of linear combinations of continuous variables and regression coefficients,
	
  * Any monotonic transformations of linear combinations of covariates and regression coefficients,
	
  * Nonlinear functions of covariates and regression coefficients,
	
  * Nonparametric formulations.
 
The Poisson surrogate model eliminates ζ_j+ from the denominator of the multinomial probabilities. This makes sense intuitively, as we do not expect the multinomial sums to provide any useful information in estimating the fixed effects.  
Given that δ̂_j can also be obtained by setting the fitted values of the multinomial sums Ŷ_j+=E(Y_j+) equal to the observed counts y_j+ in the Poisson surrogate model, δ_j has the effect of recovering the multinomial sums. The key idea is to include a separate constant δ_j for each unique combination of covariates in the Poisson surrogate models.




 §.§ Specifying Model Formulae for Poisson Surrogate Models


For purposes of exposition, the model formulae in this section are written in terms of the  language <cit.>, although this article is not concerned with software packages per se. Multinomial models are fitted using the multinom() function within the nnet package <cit.>; Poisson models are fitted using the glm() function within the stats package.


For concreteness, consider the non-parallel baseline category logit models. The “baseline category logit” assumption refers to the following: treating category 1 as the baseline category with ζ_j1 = 1  ∀ j without loss of generality, we model log(p_jq/p_j1) = log(ζ_jq) as a linear function of covariates x and regression coefficients γ. This assumption is not necessary, but chosen so that the model formulae can be illustrated using functions within the stats package. The “non-parallel” assumption refers to covariate effects that vary across categories <cit.>, i.e. all elements of the γ vector are indexed by q. That is, if the set of logits are plotted against the covariate on the same graph, a set of straight lines with slopes that are in general different will be obtained. Later we shall discuss cases where we relax this assumption. 


Consider a hypothetical dataset with two predictors X_1 and X_2 (these can be categorical or continuous) and a multinomial outcome vector Y with Q=3 categories. In short format, each row of data represents an observation with a 3-dimensional outcome vector (Y_1,Y_2,Y_3). Poisson models treat the outcomes of each observation as independent and glm() requires data to be presented in long format. This requires an additional factor C that denotes the category memberships. Each row now comprises a scalar outcome, resulting in 3 rows of data per observation. The first few rows of data in both short and long format are shown in Table <ref>. 






Table <ref> shows the equivalant relationship between non-parallel multinomial models and the corresponding Poisson models, where the parameters satisfy the usual constraints for identifiability. The Poisson surrogate models possess several important features:

	
  * The model includes an indicator variable I (that corresponds to logδ_j in Section <ref>) for each observation, although this can be simplified when all covariates are categorical. 
	This is to ensure the exact recovery of the multinomial sums, as the fixed sums are treated as random in the Poisson models. As a result, we do not interpret the coefficients of I since they are just nuisance parameters. 
	
  * The category membership indicator C enters as a covariate in the Poisson models, where the coefficients correspond to the intercepts in the multinomial modelsv so that the counts are allowed to vary by category.
	
  * The model includes interaction terms between X and C (denoted by * in the model formula), where the coefficients correspond to the slopes in the multinomial models. This is due to the non-parallel assumption where each category has a separate slope, and also the fact that multinomial models treat the response counts jointly for each observation, whereas Poisson models treat each response count as a separate observation. It is important that these interaction terms are included even if they are not significant. For multinomial models where some (partial models) or all (parallel models) of the covariate effects do not vary across categories, the equivalent Poisson models can be obtained by modifying the interaction structure between X and C accordingly. For instance, in parallel models where all categories share the same covariate effects, there is no need to include the interation terms between X and C, since the slopes do not vary across categories. 
 

 

When writing the model formula, it is important to specify I and C as factors due to their categorical nature. This can be achieved via the factor() function in . 


 

Special Case: Categorical Covariates 

When all covariates are categorical, the model formulae in Table <ref> offer a more parsimonious option for fitting the Poisson models without having to estimate a separate parameter for each observation. 

 

As stated above, the key to achieving the 1-1 correspondence between multinomial and Poisson models (with the same link function) is to include a separate constant for each unique combination of covariates. For categorical covariates, this can be achieved by including the full interaction among the predictors in the Poisson model. When all the covariates are categorical, the interaction term has the precise effect of pooling groups of observations with identical covariates. Fitting such models is equivalent to fitting the observation index as a factor (Table <ref>), but the pooling results in a smaller effective data frame, and therefore smaller storage requirements and faster fitting speed, with no loss of information. Of course, if there are many factors, there may not be much saving, because it will be comparatively rare for different observations to have all the same factor level combinations. 








§ EXTENDING THE “POISSON TRICK” FOR CORRELATED MULTINOMIAL RESPONSES





 §.§ Derivation

Consider a set of observations which fall into a collection of I groups and let λ_i = (λ_iq)_q=1^Q be a vector-valued random effect for group i. Each observation belongs to only a single group. Extending the notation in Section <ref>, the Q × 1 response vector for observation j in group i is Y_ij = (Y_ijq)_q=1^Q, with the corresponding probability vector p_ij = (p_ijq)_q=1^Q, where p_ijq = λ_iqζ_ijq/∑_q=1^Q λ_iqζ_ijq. Conditional on the multinomial sums Y_ij+ = ∑_q=1^Q Y_ijq and the random effects λ_i, 
the counts are Multinomially distributed:

    Y_ij|Y_ij+,λ_i ∼ℳ(Y_ij+, p_ij).
 
In analogy to the results in Section <ref>, given the random effects, we treat Y_ij+ as random and assume

    Y_ij+|λ_i ∼𝒫(δ_ij∑_q=1^Q λ_iqζ_ijq),

independently for each i and j. This gives

    Y_ijq|λ_iq∼𝒫(δ_ijλ_iqζ_ijq),

independently for each j and q. The probability argument in Equation <ref> still holds, now conditional on the random effects:

    ∑_i ∑_j ℓ^𝒫(δ_ij∑_q=1^Q λ_iqζ_ijq;Y_ij+|λ_i) + ∑_i ∑_j ℓ^ℳ(ζ_ij;Y_ij|Y_ij+,λ_i) = ∑_i ∑_j ∑_q ℓ^𝒫(δ_ijλ_iqζ_ijq; Y_ijq|λ_iq),
 
where  ζ_ij = (ζ_ijq)_q=1^Q. If the random effects are observed, the conditional probability statement above imply a 1-1 exact correspondence between the multinomial and the Poisson surrogate models. However, due to the unobserved nature of the random effects, interest lies in the marginal distribution, obtained by integrating out the random effects. This results in an approximate relationship between the two models. It turns out that the marginal likelihood of the approximating Poisson surrogate model (right hand side of Equation <ref>) can be expressed in closed-form if we assume an independent Gamma model for the random effects, with E(λ_iq) = α_q/β_q and Var(λ_iq) = α_q/β_q^2, i.e. λ_iq∼𝒢(α_q,β_q) <cit.>. 


With this assumption for the distribution of the random effects, the marginal likelihood of the multinomial model that we would like to fit (second term on the left hand side of Equation <ref>) is given by

    L^ℳ   = ∏_i ∫⋯∫∏_j 	 {y_ij+!/∏_q y_ijq!∏_q (λ_iqζ_ijq/∑_q=1^Q λ_iqζ_ijq)^y_ijq}×∏_qβ_q^α_qλ_iq^α_q-1e^-β_qλ_iq/Γ(α_q) dλ_i1…dλ_iQ.
 
This does not generally exhibit a closed-form solution regardless of the random effect distribution, unless in the special cases of no covariate or with only group specific covariates <cit.>. Numerical or simulation methods can be used to approximate the likelihood, with computational efforts increasing with increasing number of groups and categories. On the other hand, the marginal likelihood of the Poisson surrogate model can be expressed in closed-form:

    L^P    = ∏_i {∫⋯∫∏_j ∏_q	 
    e^-δ_ijλ_iqζ_ijq(δ_ijλ_iqζ_ijq)^y_ijq/y_ijq!∏_qβ_q^α_qλ_iq^α_q-1e^-β_qλ_iq/Γ(α_q) dλ_i1…dλ_iQ}
       = ∏_i {∏_qΓ(α_q+y_i+q)β_q^α_q/Γ(α_q)(β_q+∑_jδ_ijζ_ijq)^α_q+y_i+q×∏_j∏_q(δ_ijζ_ijq)^y_ijq/y_ijq!}.
 
The Poisson surrogate model is an extension of the Gamma-Poisson model as proposed by <cit.> and <cit.> to allow the modelling of counts for multiple categories. 


As a consequence of Equation <ref>, we have 

    E(Y_ijq) = α_q/β_qδ_ijζ_ijq.
 
Refer to the appendix for details. This is the population-averaged expected value and is not suitable for prediction in general, as it does not take into account the cluster effect. However, it can be useful for out of sample prediction, when there are no samples present in a particular group. 


Special Case: Var(λ_iq) approaches 0 

When Var(λ_iq) approaches 0 for all q, the model reduces to the special case of no random effects as outlined in Section <ref>, and the exact correspondence between the multinomial and the Poisson models can be regained. 





 §.§ Identifiability
 
There is some lack of identifiability with the model formulation given by Equation <ref>, characterized by non-uniqueness of the maximum likelihood estimates. 
There is an identifiability issue between λ_iq and δ_ij, and also between λ_iq and the category intercepts. To fix this, we impose the constraint of α_q = 1/β_q so that  E(λ_iq) = 1. As a consequence, λ_iq∼𝒢(1/β_q,β_q) and Var(λ_iq) = β_q. Also, we only require a random effect for each logit, and thus a constraint for the random effects associated with the baseline category q=1 is needed. Denote u_iq = logλ_iq. Several authors such as <cit.> (pp.514) and <cit.> considered a multivariate normal distribution for the random effects, i.e. (u_iq)_q=2^Q ∼𝒩(0,Σ), where Σ is a Q-1 by Q-1 variance-covariance matrix. This is equivalent to saying that u_i1=0 for all i, or σ_11 = 0. The equivalent statement in our proposed model is to fix λ_i1 = 1 for all i. This is tantamount to saying Var(λ_i1) = β_1 approaches 0, and thus α_1 approaches ∞. 



 §.§ Prediction of Random Effects and Fitted Values

We focus on the best predictor (BP) for random effects prediction, i.e. the predictor that minimises the overall mean squared error of prediction.  <cit.> shows that the BP is given by the posterior expectation of the random effect. Under the proposed Poisson surrogate model, the BP is given by 

    BP(λ_iq) = λ̂_iq≡λ^⋆argmin E(λ_iq-λ^⋆)^2 := E(λ_iq|y),


which can be calculated via 

    λ̂_iq = ∫_-∞^∞λ_iq f(λ_iq) f(y|λ_iq)  dλ_iq/∫_-∞^∞ f(λ_iq) f(y|λ_iq)  dλ_iq .
 
Solving for the integral, the BP is

    λ̂_iq  = Y_i+q + 1/β_q∑_j  δ_ijζ_ijq + β_q,

where Y_i+q = ∑_j Y_ijq. λ̂_iq depends on the parameters δ_ij, γ and β_q, in which we replace by their estimators, leading to the empirical best predictor (EBP). The fitted values can then be defined as 

    Ŷ_ijq = δ_ijλ̂_iqζ_ijq,

where we replace δ_ij and ζ_ijq by their respective estimators δ̂_ij and ζ̂_ijq.




 §.§ Parameter Estimation
 
Consider the parameterization ζ_ijq = exp(η_ijq) which is widely adopted in practice, where η_ijq = x_ijq^Tγ, where x_ijq and γ are both vectors. The chosen index structure for x and γ encompasses a variety of possible scenarios: (i) category-specific predictors with generic coefficients x_ijq^Tγ, (ii) category-specific predictors with category-specific coefficients x_ijq^Tγ_q, and (iii) observation-specific predictors with category-specific coefficients x_ijγ_q. This can be achieved by creating the appropriate interaction terms between the predictor and the category indicator variable, thus modifying the model matrix. Note that observation-specific predictors must be paired with choice-specific coefficients. Otherwise they will disappear in the differentiation when we consider the log-odds. 


Denote θ = (γ,(β_q)_q=2^Q), where γ includes the incidental parameters log(δ_ij) for all i and j. Algorithm <ref> presents an Expectation/Conditional Maximization (ECM) algorithm <cit.> for parameter estimation of the Poisson surrogate model. Refer to the appendix for a detailed derivation.










§ YOGURT BRAND CHOICE DATASET
 
We consider the yogurt brand choice dataset previously analyzed by <cit.> and <cit.>. <cit.> approximated the likelihood of a multinomial logit model with Gaussian random effects using a discrete distribution. <cit.> approximated the multinomial logit model using the Poisson log-linear model and Poisson nonlinear model, both with Gaussian random effects. 


The dataset consists of purchases of yogurt by a panel of 100 households in Springfield, Missouri, and were originally provided by A. C. Nielsen. The data
were collected by optical scanners for about two years and correspond to 2,412 purchases. Variables collected include brand, price and presence of newspaper feature advertisements for each purchase made by households in the panel. Price and feature advertisements are choice-specific variables. We assume a parallel baseline logit model by assigning generic coefficients γ to these variables, as we do not expect the effect of price and feature advertisements on the probability of purchase to vary according to brands. The four brands of yogurt: Yoplait, Dannon, Weight Watchers, and Hiland account for market shares of 34%, 40%, 23%, and 3% respectively. Following <cit.>, we put Hiland as the reference brand. Table <ref> presents the yogurt data in both long and short format. The letters `f' and `p' represent the feature and price variables respectively, with the letter that follows denoting the brand. For instance, `fy' stands for `feature of Yoplait' and `pd' stands for `price of Dannon'. 




We fit the models proposed in Sections <ref> and <ref>, and compare our results to that of <cit.>, fitted using the SAS macro GLIMMIX and the SAS procedure NLMIXED. The results are presented in Table <ref>. The preference ordering of the brands are the same for all models, i.e. Yoplait is the most preferred brand, followed by Dannon, Weight Watchers and Hiland. The slope parameters estimates have the expected signs for all models. An increase in price is associated with a decrease in the probability of purchase. Feature advertisement tends to increase the chance of purchase. The household-to-household variation in the probability of purchase for Weight Watchers is much larger than the other brands, although none are significant (p > 0.05). 


In comparing the estimates between models, we note that the fixed effects model is likely to produce biased estimates as it did not take into account of the correlation induced by multiple purchases from the same household. The parameter estimates of NLMIXED and Gamma-Poisson are uniformly larger than that of GLIMMIX, except for the intercept associated with Weight Watchers (for NLMIXED) and for the slope associated with price (for Gamma-Poisson). The estimates of the standard errors from NLMIXED and Gamma-Poisson are also uniformly larger than that of GLIMMIX. These differences can be attributed to the different distributional assumptions of the random effects, and also the different approximations used in GLIMMIX and NLMIXED to estimate the intractable likelihood. In this regard, our model exhibit a closed-form likelihood that allows exact inference to be performed with respect to the approximating model. 


We tried to fit a simplified version of GLIMMIX using the glmer() function within the lme4 package in , with just a random effect per household (ignoring the choice effect). However, the model failed to converge within a few months, even though <cit.> claimed that the GLIMMIX model coverged in SAS.





 






§ CONCLUDING REMARKS
 
In this article, we presented methods for fitting various multinomial regression models via the so-called “Poisson Trick” and its extensions. The “Poisson Trick” for fitting fixed effects multinomial regression models is handy when the direct fitting of multinomial models is not supported, for instance the INLA package <cit.> in . For multinomial regression models with random effects, there exist a variety of experience for using the existing extensions proposed by <cit.>, from taking months to fit a moderate sized dataset <cit.>, producing nonsense results <cit.> to non-convergence in our experience of fitting the yogurt brand choice dataset. We proposed an extension of the “Poisson Trick” using Gamma (multiplicative) random effects. In contrast to the models by <cit.>, our model exhibits a closed-form likelihood and can be maximized using existing functions for fitting generalized linear models that are stable and heavily optimized, without having to approximate the integrals.

















§ APPENDIX










 §.§ Derivation of the Population-Averaged Expected Values in Equation <ref>

Equation <ref> is also equivalent to

    ∏_i {∏_q  [ Γ(α_q+y_i+q)/Γ(α_q) y_i+q!(  ∑_j δ_ijζ_ijq/β_q+∑_jδ_ijζ_ijq) ^y_i+q( β_q/β_q+∑_jδ_ijζ_ijq)^α_q]  ×∏_q  [ y_i+q!/∏_j y_ijq!∏_j (δ_ijζ_ijq)^y_ijq/(∑_jδ_ijζ_ijq)^y_i+q] }.
 
This results in two different interpretations for the extended Gamma-Poisson surrogate model: 

	
  * For each i and q, Y_ijq is independent negative multinomial 𝒩ℳ(α_q, δ_ijζ_ijq/β_q+ ∑_j δ_ijζ_ijq) (Equation <ref>). 
	
  * For each i and q, the category sums  Y_i+q are independent Negative Binomial 
 𝒩ℬ(α_q,∑_jδ_ijζ_ijq/β_q+∑_jδ_ijζ_ijq), and conditional on the Y_i+q, Y_ijq is independent multinomial ℳ(Y_i+q,  δ_ijζ_ijq/∑_jδ_ijζ_ijq) (Equation <ref>). 
 
Taking expectation of both Equations <ref> and <ref> with respect to Y_ijq gives rise to the population-averaged expected value given in Equation <ref>. The definitions of negative multinomial and negative binomial distributions are given in the following subsections.



  §.§.§ Negative Multinomial Distribution

This is the distribution on the n+1>2 non-negative integers outcomes { X_0,…,X_n }, with corresponding probability of occurence p = { p_0,…,p_n } and probability mass function

    Γ( ∑_i=0^n x_i ) p_0^x_0/Γ(x_0)∏_i=1^n p_i^x_i/x_i!,


for parameters x_0 > 0 and p=(p_i)_i=1^n, where p_i ∈ (0,1) for all i, ∑_i=0^n p_i = 1 and Γ(·) is the Gamma function. We write Y ∼𝒩ℳ(x_0,p). For positive integer x_0, the negative multinomial distribution can be recognized as the joint distribution of the n-tuple { X_1,…,X_n } when performing sampling until X_0 reaches the predetermined value x_0. The mean vector of negative multinomial distribution is given by x_0/p_0 p.



  §.§.§ Negative Binomial Distribution

This is the distribution on the non-negative integers outcome X, with corresponding probability of occurence p and probability mass function

    Γ(r+x)/x!Γ(r) (1-p)^r p^x,

for parameters r>0 and p∈(0,1). We write X∼NB(r,p). For positive integer r, the negative binomial distribution can be recognized as the distribution for the number of heads before the rth tail in biased coin-tossing, but it is a valid distribution for all r>0. In engineering, it is sometimes called the Pólya distribution in the case where r is not integer.





 §.§ Derivation of the Expectation/Conditional Maximixation (ECM) Algorithm in Section <ref>

Treating λ = λ_iq for all i and q=2 to Q as missing data and y = y_ijq for all i, j and q as observed data, the complete data is (y_ijq,λ). Denote θ = (γ,(β_q)_q=2^Q), where γ includes the incidental parameters log(δ_ij) for all i and j. The complete data log-likeliood ℓ(θ|y,λ) is

    -∑_i∑_q≠1λ_iq( ∑_j e^x_ijq^Tγ) + ∑_i∑_q≠1 y_i+qlogλ_iq + ∑_i∑_j∑_q x_ijq^Tγ y_ijq + ∑_i∑_j∑_q log(y_ijq!) +  
    ∑_i∑_q≠1 (1/β_q-1)logλ_iq - ∑_i∑_q≠1λ_iq/β_q - ∑_i∑_q≠1logβ_q/β_q-∑_i∑_q≠1logΓ(1/β_q).
 
The (t+1)th E-step involves finding the conditional expectation of the complete data log-likelihood with respect to to the conditional distribution of λ given y and the current estimated parameter θ^(t). Straightforward algebra establishes that

    λ_iq|y_ijq,θ^(t)∼𝒢( y_i+q+1/β_q^(t) , (∑_j e^x_ijqγ^(t)+1/β_q^(t))^-1),

independently for each i and q, where the gamma distribution is parameterized in terms of scale parameter. It follows that

    λ̂_iq^(t+1) = E(λ_iq|(y_ijq)_j,θ^(t)) = y_i+q+1/β_q^(t)/∑_j e^x_ijqγ^(t)+1/β_q^(t)
    χ̂_iq^(t+1) = E(log(λ_iq)|(y_ijq)_j,θ^(t)) = ψ(y_i+q+1/β_q^(t)) - log(∑_j e^x_ijqγ^(t)+1/β_q^(t)).
 
Thus, in the (t+1)th E-step, we replace λ_iq and χ_iq = log(λ_iq) in Equation <ref> with λ̂_iq^(t+1) and χ̂_iq^(t+1), giving Q(θ|θ^(t)). The (t+1)th CM-step then finds θ^(t+1) to maximize Q(θ|θ^(t)) via a sequence of conditional maximization steps, each of which maximizes the Q function over a subset of θ, with the rest fixed at its previous value. In our application, it is natural to partition θ into γ and β_q for each q=2 to Q. Differentiating Equation <ref> with respect to γ, we obtain

    -∑_i∑_j∑_q λ_iq x_ijq e^x_ijq^Tγ + ∑_i∑_j∑_q x_ijq y_ijq,


which is the score equation of the Poisson log-linear model <cit.> with an additional offset λ_iq. This allows us to leverage existing functions for fitting generalized linear models available in most statistical software packages for maximizing γ in the CM step. This is an important feature as γ often contains a huge amount of parameters in our applications, due to the inclusion the incidental parameter log(δ_ij) for every unique combination of covariates. Existing functions for fitting generalized linear models are typically stable and heavily optimized, even for a large number of parameters. Maximizing β_q in the CM step for each q is straightforward, as it only involves univariate optimization.





§ ACKNOWLEDGEMENTS

Lee's research is partially supported by the Australian Bureau of Statistics. The authors are grateful to Zhen Chen for providing the yogurt data.








agsm 

 

