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.
115,531 characters · 8 sections · 158 citation commands
Verifying the existence of maximum likelihood estimates for generalized linear models
\pagenumbering{gobble}
\and
\and
}
\thispagestyle{empty} JEL Classification Codes: { C13, C18, C23, C25}\\ Keywords:{ Nonlinear models, GLM, Separation, Pseudo-maximum likelihood, Panel data, Gravity models}
\pagenumbering{arabic}
\setcounter{page}{1} Estimators based on count data models are widely used in applied economic research (cameron_regression_2013; winkelmann_econometric_2013). In particular, Poisson regression has exploded in popularity since the publication of santos_silva_log_2006.\footnote{Poisson pseudo-maximum likelihood (PML) estimators have emerged as a workhorse approach for studying health outcomes (manning_estimating_2001), patent citations (figueiredo_industry_2015), trade (orourke_aer2019), migration (bertoli_rational_2019), commuting (brinkman2019freeway), auctions (bajari_rand2003), finance (cohn2021count), and many other economic applications with nonnegative dependent variables and/or with data assumed to be generated by a constant-elasticity model.} Given this widespread and longstanding popularity, it is genuinely surprising that economists have only become aware relatively recently that count data models are not guaranteed to have maximum likelihood (ML) solutions. More precisely, santos_silva_existence_2010 show that the first-order conditions that maximize the likelihood of Poisson models might not have a solution if regressors are perfectly collinear over the subsample where the dependent variable is nonzero. Beyond this observation, however, santos_silva_existence_2010 caution that “it is not possible to provide a sharp criterion determining the existence” of Poisson ML estimates. Moreover, although nonexistence is a well known issue in binary outcome models, it seemingly remains unknown if similar issues could arise in other nonbinary outcome models besides Poisson, and the connections between these various cases remain seemingly unknown as well.
In this paper, we resolve several key aspects of this ambiguity. {We document that} nonexistence of ML estimates is a potential problem for a broad class of generalized linear models (GLMs), including Poisson, binary outcome models such as probit and logit, as well as several other models. {So as not to overstate our own theoretical contribution, the former result is taken from verbeek1989compactification, and similar results for related classes of models can be found in aickin1979existence, geyer1990likelihood, clarkson_computing_1991, geyer2009likelihood, and fienberg2012maximum. As we discuss below, these results are often not given enough attention. }We also clarify that this problem continues to be salient for pseudo-maximum likelihood (PML) estimators of these models and, furthermore, that some common PML estimators are affected by nonexistence in ways that cannot be remedied without changing the estimator. For cases in which simpler remedies are possible, we discuss computational methods for detecting and resolving nonexistence and propose a novel algorithm that works well even in settings that require a complex array of high-dimensional covariates, such as panel data models with multiple levels of fixed effects.
We derive our main results in part by drawing on a largely uncredited contribution by verbeek1989compactification, who established necessary and sufficient conditions governing the existence of ML estimates for a broad class of GLMs.\footnote{As of this writing, both printings of verbeek1989compactification,verbeek_compactification_1992 together have only eleven unique citations listed on Google Scholar. Another notable earlier contribution by aickin1979existence, discussed below, appears to be similarly uncredited. }
Using Verbeek's earlier results as our starting point, we show that for many GLMs, even when the ML estimates can nominally be said to “not exist”, at least some of the linear parameters can usually be consistently estimated. We also add new results for PML estimation approaches that have only become popular {in more recent years} and that turn out not to share these useful properties. For example, the log-link gamma PML estimator sometimes recommended in fields such as international trade and health care economics has very different conditions governing nonexistence than Poisson and suffers from more dire consequences when it occurs.
In addition to discussing how to detect such a problem, we also provide guidance on what can be done about it. At the moment, this is another area in need of clarity. Even for binary response models, where nonexistence is well known as the so-called “separation” problem, textbooks that mention the topic generally stop short of suggesting remedies (zorn_solution_2005,eck2018computationally). The binary outcome model literature has filled this gap primarily by presenting a choice between two main ways of solving the problem, each with its limitations. On the one hand, the most common approach is to drop a regressor from the model (zorn_solution_2005,allison2008convergence,rainey2016dealing). This is also the approach that has been discussed most often in the context of models for nonbinary outcomes (santos_silva_existence_2010,larch_currency_2017). On the other hand, dropping a regressor has implications for the estimation and identification of the other parameters, and often it is not obvious which regressor is the “right” one to drop. Thus, a leading alternative recommended for binary outcome settings is to impose a penalty on the likelihood function, often interpreted as assuming the parameters have been drawn from a known prior distribution (heinze2002solution,gelman2008weakly). These methods can be adapted to settings with nonbinary outcomes as well (firth1993bias,kosmidis2020mean), and they have the advantage that they can produce finite estimates for all of the model parameters even when separation occurs. {However, because they modify the objective being maximized, the estimates they yield are not ML estimates and thus are not directly comparable to those from the original unpenalized GLMs.} Furthermore, they are not currently compatible with models that include high-dimensional fixed effects, which are widely used in the international trade literature (head_gravity_2014,yotov_advanced_2016) and are becoming increasingly popular in applied work in general.\footnote{{ The term high-dimensional fixed effects refers to the inclusion of multiple sets of fixed effects (e.g., individual, firm, time), where at least one set contains so many categories that estimating the model with dummy variables directly is computationally impractical.} The popularity of this type of model is only likely to increase in the near future thanks to a series of computational innovations that have made models with multiple levels of fixed effects more feasible to compute (see figueiredo_industry_2015,larch_currency_2017,stammann2017fast,berge2018efficient,ppmlhdfe) as well as a growing literature on bias corrections for incidental parameter bias (see arellano_understanding_2007,fernandez-val_individual_2016).}
Our own {preferred} remedy, which {in practice} only involves withholding the separated observations from the estimation sample, is generally very simple to implement {through the new algorithms we introduce} and does not have any of these limitations.\footnote{Our recommendation to withhold separated observations from the estimation is ostensibly similar to allison2008convergence's suggestion to “do nothing”, as doing nothing could result in approximately valid estimates and inferences for at least some of the model parameters. However, in general, doing nothing could result in lack of numerical convergence or---in the worst case---convergence to incorrect values. In the words of geyer2019slides, “no one knows how much applied statistics is garbage because of this.” Also, some software packages drop separated observations by default (e.g., Stata's probit command), but they generally are not adept at detecting these observations, nor do they usually provide theoretical justification for this practice in their documentation. Our companion \href{https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/README.md}{website} offers examples and discussion; see \url{github.com/sergiocorreia/ppmlhdfe/blob/master/guides}.} {The theoretical justification for withholding these observations} comes from an insight advanced independently by aickin1979existence, verbeek1989compactification, geyer1990likelihood, and clarkson_computing_1991: a model suffering from separation can often be nested within a “compactified” model where the conditional mean of each observation is allowed to go to its boundary values. The (pseudo-)likelihood function always has a maximum somewhere in the compactified parameter space; thus, we can transform the problem of nonexistence to one of possible corner solutions. More importantly, observations with a conditional mean at the boundary in the more compactified model are effectively perfectly predicted observations. These observations offer no information about the parameters with interior solutions and, as we will show, can be quickly detected even for very complex models. Removing these observations then results in a standard (non-compactified) version of the model that is assured to produce the same model fit as the compactified version, as well as the same point estimates and inferences of the parameters with interior solutions. {As such, our approach is equivalent in practice to fitting what geyer2009likelihood calls the “limiting conditional model” that conditions on the separated observations. }We also show that {the estimates for the estimable parameters} are consistent and that correct inference requires only careful attention to which of the regressors are involved in separation. The resulting output, on the whole, is no different than what one would observe with a perfectly collinear regressor, and the problems of interpretation and inference turn out to be very similar as well.
{Although separation becomes equivalent to perfect collinearity after excluding the separated observations, the two concepts differ in important aspects. In the case of separation, all regressors are important to the fit of the model, including those without finite estimates. Furthermore, while the separated observations are withheld from the estimation step, the model yields predicted values for them that are consistent with fitting the model over the full sample. For models estimated via ML, it is also possible to obtain one-sided confidence bounds for the parameters whose estimates diverge to infinity. For discrete-response GLMs with canonical links, one can even obtain meaningful inferences on the predicted values of the separated observations (see eck2018computationally).}
{ Since the possible non-existence of a finite MLE is a well known problem in the context of binary outcome models, it may be surprising that equally fundamental results for Poisson regression and other nonbinary GLMs remain undercited by comparison. This imbalance reflects both history and practice: the problem was easier to visualize geometrically for binary logit and probit models and more frequently encountered in applied research, leading to its widespread acknowledgment. By contrast, the analogous results for Poisson and other nonbinary GLMs have mostly been circulated in the statistical theory literature. This literature begins with haberman_log-linear_1973,haberman_analysis_1974's derivation of a necessary and sufficient condition for the existence of estimates for log-linear frequency table models. Notably, it was known at the time that this condition was difficult to verify for higher-dimensional tables (see albert_existence_1984), a still-unsettled problem we indirectly solve in this paper. Soon thereafter, wedderburn_existence_1976 independently derived a sufficient but not necessary condition for the existence of estimates across a wide class of GLMs that included Poisson and Gamma.\footnote{His result can be shown to be equivalent to santos_silva_existence_2010's later result for the Poisson model.} The first statement of a necessary and sufficient condition for the existence of Poisson regression estimates was given by aickin1979existence, who derived results for models in the linear discrete exponential family.}
{ silvapulle_existence_1981 and albert_existence_1984 are then credited with popularizing the concepts of “separation” and “overlap” for binary outcome models. silvapulle_existence_1981 gave these terms algebraic meaning, while albert_existence_1984 provided an influential geometric interpretation that distinguished between “complete” and “quasi-complete” separation. A few years later, silvapulle_existence_1986 showed how the conditions for separation studied in silvapulle_existence_1981 can be formulated as a linear programming problem, thus marking an important step towards detecting the issue in practice. Finally, verbeek1989compactification, geyer1990likelihood, and clarkson_computing_1991 independently extended these ideas to broader classes of models. Verbeek (1989) was the first to unify the conceptualizations of the problem that were being used in the binary outcome literature and the more general GLM setting, while geyer1990likelihood and clarkson_computing_1991 developed similar insights for the linear exponential family and for “models with a linear part” (stirling1984iteratively), respectively.\footnote{{The “models with a linear part” concept studied in clarkson_computing_1991 is a broader category than GLMs, as it includes any other models where the parameters enter the likelihood function via a linear index function. The Tobit model for censored data is an example of a non-GLM model conforming to this framework.}} Despite these advances, awareness of the possible nonexistence of estimates for nonbinary GLMs did not spread widely in applied work until it was highlighted in santos_silva_existence_2010. }
We add to this earlier literature in three main ways. First, by considering an expanded set of estimation approaches, we offer a more detailed treatment of how the separation problem varies across GLMs used {in economics research} and estimators thereof. For example, a significantly stricter set of conditions governs the existence of estimates for gamma PML and inverse Gaussian PML than for Poisson, logit, and probit---a result that raises concerns about applications of the former estimators to settings where zero outcomes are common, such as health care cost analysis and international trade. Importantly, these are precisely the settings where these estimators have come into common usage; see, e.g., manning_estimating_2001,egger_glm_2015.\footnote{manning_estimating_2001 leave aside the issue of zero outcomes in their paper, but indicate that gamma PML is generally a good model for health care cost data and also remark that “there is ostensibly nothing in the above analysis that would preclude applications to data where realizations of $y$ are either positive or zero, as is common in many health economics applications.” Our own findings indicate that zeroes do pose a distinct problem for gamma PML estimation that must be carefully taken into account.} Second, we clarify that at least some of the linear parameters can be consistently estimated in the presence of separation as well as how to obtain valid asymptotic inferences---though, again, it is important to note these results do not extend to all of the estimators we consider.\footnote{gourieroux_pseudo_1984 (1984, Appx 1.1) and fahrmeir1985consistency (Sec. 2.2) both assume in their proofs of consistency that the solutions for the linear parameters are interior. We present a proof that relies on a suitable reparameterization of the separated model such that the results of gourieroux_pseudo_1984 apply directly.}
Finally, we introduce a simple-but-powerful method for detecting separation in models with a large number of fixed effects, a conceptually nontrivial task that would ordinarily require solving a high-dimensional linear programming problem. Because our algorithm relies on repeated iteration of a least-squares-with-equality-constraints regression, it can take advantage of the recent innovations of correia_linear_2017, who shows how to solve high-dimensional least-squares problems in nearly linear time. To our knowledge, the only other method that has been suggested for detecting separation in large ML settings is that of eck2018computationally. {While their primary focus is on obtaining inference {for the separated observations} when the ML estimate lies in the Barndorff-Nielsen completion (barndorff1978information}), their method also detects separation by computing the null eigenvectors of the Fisher information matrix, which reveals the estimable and non-estimable components. In contrast, our methods avoid large matrix operations altogether, which should make them substantially more scalable. {In addition, their approach detects separation ex post after first running the estimation algorithm on the full data set, while ours aims to detect separation ex ante. {As such, our contribution concerns the task of detecting the separated observations rather than the downstream problem of how to do inference on them like in eck2018computationally.} Due to its scalability, our approach can provide a solution for models that fall outside the scope of currently available implementations for existing remedies, including not only detection methods such as eck2018computationally and Kosmidis2021detectseparation but also penalized likelihood methods such as gelman2008weakly and kosmidis2020mean.}
The rest of the paper proceeds as follows. Section 2 formally establishes the problem of separation in GLMs, including its sufficient and necessary conditions. Section 3 discusses how to address separation in setups with and without fixed effects. Section 4 provides an empirical example. Section 5 concludes. Further details are available in the Appendix, including additional proofs and results of interest. We have also created a \href{https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/README.md}{website} dedicated to the separation problem, which provides numerous examples illustrating the methods and principles described in this paper. {These examples include demonstrations for logit and multinomial logit taken from the literature as well as 17 examples for Poisson and Poisson PML that we ourselves have curated.}
The class of GLM-based estimators we consider is defined by the maximization of the following log (pseudo-)likelihood objective function, corresponding to distributions of the exponential family:
For brevity, we will generally refer to this objective function as the “likelihood”, though we will use “pseudo-likelihood” when strictly discussing PML estimators. The individual term $\ell_i$ will be the “likelihood contribution” or “pseudo-likelihood contribution”. $y_{i}\ge0$ is an outcome variable, $x_{i}$ is a set of $M$ regressors ($x_{1},x_{2},\ldots,x_{M}$), and $\beta\in\mathbb{R}^{M}$ is an $M\times1$ vector of parameters to be estimated. The function $\alpha_{i}(\varphi)>0$ is usually of the form $w_{i}/\varphi$, where $w_{i}$ is a known weight, and $\varphi$ is a potentially unknown scale or dispersion parameter.\footnote{$\varphi$ is not associated with the problem of separation and will henceforth be treated as known. The results we document apply to models with unknown scaling factors without loss of generality. Table (ref) gives examples. For more information on this class of models, see McCullaghNelderGLM Section 2.2.2.} $\theta_{i}=\theta(x_{i}\beta;\nu)$ is the canonical location parameter, which links the linear predictor of a given observation $x_{i}\beta$ to its likelihood contribution $l_{i}$ and its conditional mean $\mu_{i}\equiv E[y_{i}|x_{i}]=b^{\prime}(\theta_{i})$. Note that $\theta(x_{i}\beta;\nu)$ is continuous, strictly increasing, and twice differentiable in $x_{i}\beta$ and that $b(\theta_{i})$ is continuous, increasing, and convex in $\theta_{i}$. Notably, these last few restrictions together ensure that the quantities $\theta_{i}$, $x_{i}\beta$, and $\mu_{i}$ are each increasing with respect to one another and that $l\left(\beta\right)$ is continuous in $\beta$. We further assume that $\lim_{x_{i}\beta\rightarrow-\infty}\mu_{i}=0$ to rule out the simple linear model, which always has a solution.\footnote{Also note that the linear predictor term $x_{i}\beta$ is often denoted as $\eta_{i}$. We keep it as $x_{i}\beta$ to economize on notation.} $\nu$ is an additional dispersion parameter that allows us to also consider negative binomial models (see Table (ref)). Lastly, $c(\cdot)$ is a known real-valued function that depends on the specific GLM.
The first-order condition for the $m$-th individual parameter, $\beta_{m}$, follows from the GLM score function:
Examples of models conforming to this framework notably include binary outcome models, count models such as Poisson and negative binomial, and a variety of other closely related, non-GLM models such as conditional logit, multinomial logit, and the Cox proportional-hazards model. In addition, as the score vectors of many of these models can also be used to construct PML estimators for continuous data, this framework also applies to PML estimators such as Poisson, gamma, Gaussian, inverse Gaussian, and Bernoulli PML without loss of generality. Note that our interest in PML estimators represents an important deviation from verbeek1989compactification because PML estimation does not impose any restrictions on $c(y_{i},\varphi)$. As such, we can consider potential nonexistence problems in models where $y_{i}=0$ values would otherwise be inadmissible, such as log-link gamma PML and other PML estimators with similar score functions.\footnote{For more on the wide applicability of PML, see gourieroux_pseudo_1984, manning_estimating_2001, and santos_silva_log_2006.}
On top of these general restrictions, we use two further assumptions to derive a necessary and sufficient condition for existence that holds across most of these estimators. First, we assume that the matrix of regressors \ensuremath{X=x_1, x_2, \ldots, x_M}\xspace is of full column rank. This rank assumption allows us to set aside the more widely understood case of perfectly collinear regressors, although in Section (ref), we will find it useful to draw a comparison between nonexistence and perfect collinearity. Second, we assume for the moment that the individual likelihood contributions $l_{i}(\beta)$ have a finite upper bound. Later on, we will consider two estimators for which $\ell_i$ is not guaranteed to have a finite upper bound, gamma PML and inverse Gaussian PML. We show that the relevant criteria governing existence are not the same as when this assumption is met.
To extend and generalize the earlier result from santos_silva_existence_2010 for Poisson models, we are now ready to prove the following proposition:
The proof of this proposition follows verbeek1989compactification, while also drawing on an earlier proof by silvapulle_existence_1981 specifically for binary outcome models.\footnote{In verbeek1989compactification, the relevant theorems are Theorem 6, which establishes conditions under which the likelihood function has a local maximum that lies on the boundary of the parameter space, and Theorem 4, which establishes that any local maximum on the boundary is a global maximum if the likelihood function is concave. Note that we have relaxed the concavity assumption since it is straightforward to show the weaker result that if there is a local maximum at the boundary, the global maximum can only occur at the boundary.} In addition, the necessity of the condition on the boundedness of $l_{i}(\cdot)$ function is due to clarkson_computing_1991; note that Proposition (ref) later in the paper explores the implications of relaxing this assumption.
The general idea is that we want to show that if a vector $\gamma^{*}$ satisfying (ref)-(ref) exists, then the likelihood function $l(\beta)$ will always be increasing if we search for a maximum in the direction associated with $\gamma^{*}$.\footnote{This is the same concept as what geyer2009likelihood calls the “generic direction of recession”. } Otherwise, if no such $\gamma^{*}$ exists, then searching in any direction from any starting point in $\mathbb{R}^{M}$ under the noted conditions will cause $l(\beta)$ to eventually decrease, such that the function must reach a maximum for some finite $\beta^{MLE}\in\mathbb{R}^{M}$.
To proceed, let $\gamma=(\gamma_{1},\gamma_{2},\ldots,\gamma_{M})\in\mathbb{R}^{M}$ be an arbitrary nonzero vector of the same dimension as $\beta$ and let $k>0$ be a positive scalar. Now consider the function $l(\beta+k\gamma)$, which allows us to consider how the likelihood changes as we search in the same direction as $\gamma$ from some initial point $\beta$. Differentiating $l(\beta+k\gamma)$ with respect to $k$, we obtain
Suppose there is a $\gamma^{*}$ such that $z_{i}=x_{i}\gamma^{*}$ satisfies (ref)-(ref). In this case, setting $\gamma=\gamma^{*}$ the above expression becomes
with the inequality following because $b^{\prime}$ and $\theta^{\prime}$ are both positive and because $b^{\prime}=\mu<\overline{y}$. Notice also that the inequality is strict because we must have at least one observation for which $z_{i}\neq0$; otherwise, our full rank assumption would be violated, and we would be in the case of perfect collinearity. Because this expression is always positive, $l(\beta+k\gamma^{*})>l(\beta)$ for any $k>0$ and for any $\beta\in\mathbb{R}^{M}$. Thus, there is no finite solution $\beta^{MLE}\in\mathbb{R}^{M}$ that maximizes $l\left(\cdot\right)$, and estimates are said not to exist. Intuitively, the objective function will always be increasing as either $x_{i}\beta\rightarrow-\infty$ for at least one observation where $y_{i}=0$ or $x_{i}\beta\rightarrow\infty$ for at least one observation where $y_{i}=\overline{y}$.
Alternatively, suppose that, for any $\gamma$, we always have that $x_{i}\gamma\neq0$ for at least one interior observation ($0<y_{i}<\overline{y}$). Importantly, this ensures that $\lim_{k\rightarrow\infty}l_{i}(\beta+k\gamma)=-\infty$ for at least one observation. Since $l_{i}(\cdot)$ is continuous in $\beta$ and (by assumption) has a finite upper bound, we can therefore always identify a finite scalar $\overline{k}$ such that $k>\overline{k}$ implies that $l(\beta+k\gamma)=\sum_{i}l_{i}(\beta+k\gamma)<\sum_{i}l_{i}(\beta)=l(\beta)$, for any $\beta,\gamma\in\mathbb{R}^{M}$. In other words, searching for an ML or PML solution $\beta^{MLE}$ in any direction from any starting point in $\mathbb{R}^{M}$ space will always eventually yield a decrease in $l\left(\cdot\right)$. Because $l\left(\cdot\right)$ is continuous, this guarantees the existence of a finite $\beta^{MLE}\in\mathbb{R}^{M}$ maximizing $l\left(\cdot\right)$.
Next, note that, for any $y_{i}=0$ observation such that $x_{i}\gamma>0$ , $l_{i}(\beta+k\gamma)$ is monotonic in $k$ with $\lim_{k\rightarrow\infty}l_{i}(\beta+k\gamma)=-\infty$. Similarly, note that $\mu<\overline{y}$ ensures the same is true for any $y_{i}=\overline{y}$ observation such that $x_{i}\gamma<0$.\footnote{To be clear, if $\overline{y}=\infty$, we never have that $y_{i}=\overline{y};$ only conditions (ref) and (ref) are salient. On the other end of the spectrum, models for “fractional” data such as Bernoulli PML (cf., papke_econometric_1996,silva2014estimating) allow the dependent variable to vary continuously over $[0,1]$. For these models, all three conditions stated in Proposition (ref) are relevant.} Thus, we can again always find a sufficient $\overline{k}$ such that $k>\overline{k}$ implies $l(\beta+k\gamma)<l(\beta)$ so long as we always have that either $x_{i}\gamma>0$ for at least one observation where $y_{i}=0$ or $x_{i}\gamma<0$ for at least one observation where $y_{i}=\overline{y}$.\footnote{Readers should be wary of the weight carried by the word “always” here. It could be the case, for example, that $x_{i}\gamma=0$ for all $y_{i}>0$ with $x_{i}\gamma\geq0$ for all $y_{i}=0$. This is still a case where estimates do not exist, since $\gamma^{*}=-\gamma$ would satisfy the needed conditions.} Finally, note that we do not consider the case where there exists a vector $\gamma$ such that $x_{i}\gamma=0$ for all $i$, as this is the case where $X$ is not of full rank. Therefore, the only possible scenario in which estimates do not exist is the one where we can find a linear combination of regressors $z_{i}=x_{i}\gamma^{*}$ satisfying (ref)-(ref). $\blacksquare$
To tie in some standard terminology from the binary outcome literature (cf., albert_existence_1984), we will say that when estimates maximizing (ref) do not exist, the linear combination of regressors defined by $z_{i}=x_{i}\gamma^{*}$ “separates” the observations for which $z_{i}\gtrless0$ from the rest of the sample. For the sake of providing a more unified perspective, we will henceforth refer to the nonexistence with the term “separation”. A particular point of interest for us is how to also adapt the related terms “complete separation” and “quasi-complete separation” to this more general context. For binary outcome models, separation is usually considered “complete” if either $z_{i}<0$ for all $y_{i}=0$ or $z_{i}>0$ for all $y_{i}=\overline{y}=1$, since in these cases the value of $z_{i}$ perfectly predicts whether $y_{i}$ is $0$ or $1$. Otherwise, we have only “quasi-complete separation”, where only some $y_{i}$ outcomes are perfectly predicted. Outside of binary outcome models, however, as long as $y_{i}$ takes on at least two positive values, it will never be the case that $z_{i}=0$ perfectly predicts all positive $y_{i}$, regardless of whether $z_{i}<0$ perfectly predicts all $y_{i}=0$ outcomes or only some of them. Thus, for lack of an analogous vocabulary for discussing separation in the nonbinary outcome case, we would suggest that separation occurring in these models should generally be regarded as “quasi-complete”.
In addition, for those readers more familiar with santos_silva_existence_2010's results for Poisson models specifically, another term that is useful for us to clarify for the nonbinary outcome context is “overlap”. In santos_silva_existence_2010, Poisson estimates are shown to exist so long as there are no regressors that are perfectly collinear over the subsample where $y_{i}>0$. In our way of phrasing the issue, this criterion equates to saying there exists no linear combination of regressors satisfying equation (ref). However, as santos_silva_existence_2010 are careful to note, this criterion is only sufficient, rather than necessary and sufficient. As the remaining elements of the preceding proof show, even if such a linear combination exists, separation is still avoided if $z_{i}$ takes on both positive and negative values when $y_{i}=0$, such that its maximum and minimum values over $y_{i}=0$ “overlap” the $z_{i}=0$ values it takes on when $y_{i}>0$.\footnote{The question of when overlap occurs is precisely the point left ambiguous in santos_silva_existence_2010. See page 311 of their paper.}
Interestingly, we do not know of a widely accepted label for what we have called the “linear combination of regressors that separates the data” (i.e., $z_{i}=x_{i}\gamma^{*}$). Clearly, $z$ plays a central role in the analysis of separation, and the literature could use a concise name for it. We propose the term “certificate of separation.” The idea is that we can easily “certify” whether any such $z$ separates the data by verifying that i) its values satisfy (ref)-(ref), and ii) that the $R^{2}$ of a regression of $z$ against the regressors $x$ is equal to one.\footnote{Our proposed name borrows from optimization, where phrases such as “certificate of feasibility”, “certificate of convexity”, “certificate of nonnegativity”, and so on are used with a similar purpose.} Note that there can be multiple $z$'s certifying separation of different observations, and that adding up two or more $z$'s preserves their properties.\footnote{If $z^{*}$ and $z^{**}$ are valid certificates of separation with associated coefficients $\gamma^{*}$ and $\gamma^{**}$, then $z^{***}=z^{*}+z^{{**}}$ is also a valid certificate of separation, as (ref)-(ref) hold trivially and $z^{***}$ is a linear combination of $x$ with associated coefficient $\gamma^{{***}}=\gamma^{*}+\gamma^{**}$.} Thus, we refer to an “overall certificate of separation” $\overline{z}$ that can be used to identify all separated observations. For any $\gamma^{*}$ associated with a certificate of separation, we will tend to use the term “separating vector” (although another name for it is the “direction of recession”; see geyer2009likelihood). We will use $\overline{\gamma}$ to denote the separating vector associated with $\overline{z}$.
Results for gamma PML and inverse Gaussian PML. One stipulation that sticks out in Proposition (ref) is our requirement that the individual likelihood contribution $l_{i}(\cdot)$ have a finite upper bound. To our knowledge, the implications of relaxing this assumption have not been touched upon in the prior literature. Rewinding some of the last few details behind the above proof, the specific role played by this restriction is that it ensures that if $\lim_{k\rightarrow\infty}l_{i}(\beta+k\gamma)=-\infty$ for any $i$, the overall objective function $l(\beta+k\gamma)=\sum_{i}l_{i}(\beta+k\gamma)$ also heads toward $-\infty$ for large $k$. However, this might not hold if $l_{i}(\cdot)$ is not bounded from above. In this case, even if the data exhibit “overlap” (as defined above), this alone will not be sufficient to ensure that $l(\cdot)$ has a maximum. Instead, stronger conditions may be needed.
For illustration, the two models we will consider where $l_{i}(\cdot)$ does not necessarily have a finite upper bound are gamma PML and inverse Gaussian PML.\footnote{Note that ML estimation of either a gamma distribution or an inverse Gaussian distribution will not admit $y_{i}=0$ values. Thus, we consider PML versions of these estimators only. In general, what gamma PML and inverse Gaussian PML have in common is that their score functions place a relatively larger weight on observations with a smaller conditional mean. Similar results will apply to other estimators with comparable score functions. } As shown in Table (ref), the form of the pseudo-likelihood function for gamma PML regression is
and the form for inverse Gaussian PML is
In either case, notice that the associated $b_{i}$ function from equation (ref), $x_{i}\beta$ for gamma and $-\exp(-x_{i}\beta)$ for inverse Gaussian, has a lower bound of $-\infty$, as $\lim_{x_{i}\beta\rightarrow-\infty}b_{i}=-\infty$. Thus, in either case, while $l_{i}(\cdot)$ continues to have a finite upper bound for observations where $y_{i}>0$, if $y_{i}=0$, we have that $\lim_{x_{i}\beta\rightarrow-\infty}l_{i}(\cdot)=\infty$. With this in mind, the following Proposition collects results that apply to either of these estimators:
Part (a) of Proposition (ref) follows from again considering the function $l(\beta+k\gamma)$, this time specifically for gamma PML. Using (ref), it is straightforward to show that $\lim_{k\rightarrow\infty}l(\beta+k\gamma)=-\infty$ if $x_{i}\gamma<0$ for at least one observation with $y_{i}>0$. By a continuity argument similar to the one used above, this implies that $l(\beta+k\gamma)$ must eventually become decreasing in $k$ for sufficiently large $k$.
Next, consider what happens if there exists a linear combination of regressors $z_{i}=x_{i}\gamma$, which is always $\ge0$ when $y_{i}>0$. In this case, because $\lim_{k\rightarrow\infty}\sum_{z_{i}\neq0}-\alpha y_{i}\exp(-x_{i}\beta-kz_{i})=0$, we have that \[ \lim_{k\rightarrow\infty}l(\beta+k\gamma)=\lim_{k\rightarrow\infty}\sum_{i}-\alpha\left(x_{i}\beta-kz_{i}\right)+\sum_{z_{i}=0}-\alpha y_{i}\exp(-x_{i}\beta). \] There are four possibilities for the above limit. If $\sum_{i}z_{i}<0$, the gamma pseudo-likelihood function is always increasing in the direction associated with $\gamma$, such that finite estimates do not exist. Alternatively, if $\sum_{i}z_{i}>0$, the limit equals $-\infty$ and we are again assured that this function must eventually decrease with $k$, such that estimates will exist.
The remaining two possibilities occur when $\sum_{i}z_{i}=0$. In this case, the effect of an increase in $k$ on the likelihood function is always given by \[ \frac{dl(\beta+k\gamma)}{dk}=\sum_{z_{i}>0}\alpha y_{i}\exp(-x_{i}\beta-kz_{i})z_{i}\ge0. \] Inspecting the above expression, $dl(\beta+k\gamma)/dk>0$ with strict inequality if $z_{i}>0$ for at least one observation with $y_{i}>0$, ensuring again that finite estimates do not exist. The final possibility is if $z_{i}=0$ for all $y_{i}>0$ observations, in which case $l(\beta+k\gamma)=l(\beta)$ for any $k>0$. In other words, regardless of which initial $\beta$ we consider, the likelihood will always be weakly higher when we increment $\beta$ by some positive multiple of $\gamma$, implying either that a finite solution for $\beta$ maximizing $l(\cdot)$ does not exist (if $\sum_{i}z_{i}=0$ with $z_{i}>0$ for at least one $y_{i}>0$) or that any finite solution will be nonunique (if $z_{i}=0$ for all $y_{i}>0$). Thus, taking all of these results together, gamma PML estimation of $\beta$ will not have a finite solution if there exists a linear combination of regressors satisfying (ref) and (ref) and may not necessarily have a unique solution even if these conditions are not met.
For proving part (b), which pertains instead to inverse Gaussian PML, it is again convenient to work with the derivative of the $l(\beta+k\gamma)$ function with respect to $k$. Continuing to let $z_{i}=x_{i}\gamma$, and after dividing up terms appropriately, this derivative can be expressed as
Let us start with the conditions highlighted in part (b), where $z_{i}\ge0$ for all $y_{i}>0$ and where $z_{i}<0$ for at least one observation where $y_{i}=0$. We can see that the second and third terms in (ref) will go to $0$ in the limit where $k$ becomes infinitely large. The first term, meanwhile, heads to infinity. Thus, the pseudo-likelihood function increases asymptotically for large $k$, and it is clear there is no finite solution for $\beta$.
However, we still need to verify what happens if we cannot find a linear combination $z_{i}$ satisfying both of the conditions stated in part (b). This part requires slightly more work. If $z_{i}\ge0$ for all $i$, for example, all three terms in (ref) go to zero for $k\rightarrow\infty$\textemdash a result that is not in itself all that informative. Likewise, if we consider what happens when $z_{i}$ may be less than zero for $y_{i}>0$, the first and third terms could potentially head toward $+\infty$, while the second term heads toward $-\infty$. In all of these seemingly ambiguous scenarios, we can use L'H\^{o}pital's rule to clarify that $dl(\beta+k\gamma)/dk<0$ for sufficiently large $k$, indicating that the pseudo-likelihood function will always eventually decrease in the direction associated with $\gamma.$ $\blacksquare$
To our knowledge, we are the first to study the general circumstances under which estimates from gamma PML and inverse Gaussian PML exist.\footnote{Even wedderburn_existence_1976, in his original derivation of a sufficient condition for the existence of GLM estimates, specifically avoids commenting on what conditions would be needed for gamma estimates to exist if the dependent variable is allowed to be zero.} That these estimators have not been specifically looked at in this context is perhaps not all that surprising, since these models have not traditionally been used with zeroes and since the increase in popularity of PML estimation in applied work has only occurred relatively recently. Indeed, thanks to contributions such as manning_estimating_2001, santos_silva_log_2006, and head_gravity_2014, the main context in which researchers will likely be familiar with gamma PML is in settings where zeroes are common, such as data for international trade flows and health care costs. Inverse Gaussian PML is also sometimes considered for these types of applications (see egger_glm_2015) but is significantly less popular, likely because it is more difficult to work with numerically.
In this light, the results contained in Proposition (ref) can be read in one of two ways. On the one hand, we confirm that gamma PML and inverse Gaussian PML can, in principle, be used with datasets that include observed zeroes, even though their ML equivalents cannot. Since the ability to admit zeroes on the dependent variable is one of the reasons researchers have recently become curious about these estimators, this confirmation seems useful.\footnote{The other main reason is that the traditional practice of applying a log transformation to the dependent variable and estimating a linear model is now widely known to introduce a bias whenever the error term is heteroskedastic.} On the other hand, we can see from a comparison of Propositions (ref) and (ref) that the criteria required for gamma PML and inverse Gaussian PML to have finite solutions are considerably more strict than the equivalent criteria required for most other standard GLM estimators. Furthermore, as we will see in the next section, these fundamental differences also imply that gamma PML and inverse Gaussian PML lack some appealing properties that enable us to more easily remedy situations where estimates do not exist for other models. For these reasons, we recommend researchers to exercise extra caution when using either of these two estimators with datasets that include zeroes in the dependent variable.
This section describes recommendations for dealing with separation in practice, including in high-dimensional environments with many fixed effects and other nuisance parameters. Before digging into these details, it is important to make two general points. First, as we have shown, the implications of separation differ depending on the estimator; thus, the appropriate remedy should similarly depend on the estimator being used. Second, the appeal of our own preferred alternative\textemdash {withholding the separated observations from the estimation sample beforehand}\textemdash is likely to depend on one's comfort level with allowing the linear predictor $x_{i}\beta$ to attain what would ordinarily be an inadmissible value. One method we caution against is simply removing from the model one of the regressors involved in the separation, as this affects the identification and estimation of all remaining parameters, with the effect differing depending on which regressor is dropped.
In subsection (ref), we will first show that when $x_{i}\beta$ is allowed to attain $\pm\infty$, separated observations often do not affect the score function for $\beta$ under fairly general circumstances. This insight provides a theoretical justification for the practice of {withholding} separated observations from the estimation, at which point the separation problem becomes one of perfect collinearity within the remaining estimation sample. This insight is particularly useful for models with many fixed effects, as perfect collinearity amongst the fixed effects is generally not a problem for identifying the coefficients of the non-fixed effect covariates.\footnote{{Indeed, one of the computational advantages of modern software packages for estimating models with high-dimensional fixed effects (e.g., correia_linear_2017) that they do not explicitly estimate unique coefficients for each fixed effect parameter. This usefully means the researcher does not need to concern themselves with which fixed effects may be collinear when implementing the estimation.} } Once these results are established, subsections (ref) and (ref) then focus on detecting and addressing separation, including in high-dimensional environments.
We now turn to discussing how identification of at least some of the model parameters can be achieved when separation occurs. We start with the concept utilized in aickin1979existence, verbeek1989compactification, geyer1990likelihood, and clarkson_computing_1991 of a “compactified” (or “extended”) GLM where the parameter space is extended to admit its boundary values. We can phrase this compactification in one of several equivalent ways. For example, we could express the domain for $\beta$ as $[-\infty,+\infty]^{M}$, the compact closure of $\mathbb{R}^{M}$.\footnote{As discussed in verbeek1989compactification, one way to justify the inclusion of infinitely large values in the admissible parameter space is to observe that we could just as easily perform the maximization over a homeomorphic space where the parameters of interest are instead bounded by a finite interval (e.g., $[-1,1]^{M}$ instead of $[-\infty,\infty]^{M}$). A version of this concept is also described in haberman_analysis_1974. It is also sometimes referred to as the “Barndorff-Nielsen completion” (barndorff1978information).} However, it is also often convenient to work with the linear predictor $x_{i}\beta$, which in turn also may vary over $[-\infty,+\infty]$ for each $i$. In particular, note how the conditional mean $\mu_{i}$ behaves as $x_{i}\beta$ attains either of its two limits: when $x_{i}\beta\rightarrow-\infty$, we have that $\mu_{i}\rightarrow0$, whereas when $x_{i}\beta\rightarrow\infty$ (a situation that is only relevant for binary response models and fractional data models), we have that $\mu_{i}\rightarrow\overline{y}.$ { It is straightforward to show that estimates for $\beta$ maximizing the likelihood always exist when we compactify the model in this way.}
With this adjustment to the parameter space in mind, consider what happens to the score function $s(\beta)$ and information matrix $\mathbf{F}(\beta):=\mathbb{E}[\partial s(\beta)/\partial\beta]$ in the limit as $k\rightarrow\infty$ in the case of separation outlined above. In other words, consider
and
where we take $\gamma^{*}$ to be a vector satisfying the applicable conditions for nonexistence. At this point, it will also be useful to state the following lemma:
The utility of this lemma (which we prove in our Appendix) is that, together with (ref) and (ref), it delivers the following proposition:
Part (a) follows from our proof of Proposition (ref) (and is also a central result from verbeek1989compactification). After allowing $\beta$ to take on either $-\infty$ or $+\infty$, we rule out the cases where estimates would otherwise be said not to exist. Parts (b) and (c) are analogous to the insights contained in clarkson_computing_1991's Theorem 2. After invoking Lemma (ref), the score function in (ref) can be rewritten as
The key insight presented in (ref) is that the contribution of any observation with $x_{i}\gamma^{*}\neq0$ always drops out of the overall score function under these circumstances.
As a result, it must be the case that any $\beta$ that maximizes $l(\beta)$ in the compactified model must also maximize $\sum_{x_{i}\gamma^{*}=0}l_{i}(\beta)$ (i.e., the likelihood associated with the observations not separated by $x_{i}\gamma^{*}$). Otherwise, we would have that $\sum_{x_{i}\gamma^{*}=0}s_{i}\left(\beta\right)\neq0$ and $\sum_{x_{i}\gamma^{*}\neq0}s_{i}\left(\beta\right)=0$, implying that the joint likelihood of the non-separated observations can be increased without affecting that of the separated observations.
Parts (b) and (c) then follow because if $\beta=\beta^{*}$ maximizes the likelihood of the non-separated observations $\sum_{x_{i}\gamma^{*}=0}l_{i}(\beta)$, then any coefficient vector of the form $\beta^{*}+k\gamma^{*}$ maximizes it as well. That is, the estimates of $\beta$ will be different with the separated observations than without them. But the quantities $x_{i}\beta$, $\theta_{i}$, and $\mu_{i}$ will not be affected, as stated in part (b), since $x_{i}(\beta^{*}+k\gamma^{*})=x_{i}\beta^{*}$ over the subsample where $x_{i}\gamma^{*}=0$ (and since $\theta_{i}$ and $\mu_{i}$ are functions of $x_{i}\beta)$. Consequently, for any $m$ such that $\gamma_{m}^{*}=0$, the individual parameter estimate $\beta_{m}^{*}+k\gamma_{m}^{*}=\beta_{m}^{*}$ is clearly the same in either case, as stated in part (c).
{Part (d) then follows from similar arguments but relies on a more detailed explanation. In short, valid inference must account for the fact that some regressors become collinear when estimation is restricted to the non-separated observations. In our full proof of part (d), provided in our Appendix, we show that withholding the separated observations from the estimation is equivalent to estimating the finite components of a re-parameterized model that explicitly introduces the certificates of separation as regressors with infinite coefficients. A valid information matrix can then be formed using only the regressors that have finite estimates in this re-parameterized model, which importantly preserves all model coefficients that have finite estimates in the original model. Consequently, standard inference procedures can be used for these coefficients.}$\blacksquare$
For researchers encountering separation problems, the key takeaways from Proposition (ref) are likely to be parts (c) and (d): even if one or more of the elements of the MLE for $\beta$ “does not exist” (i.e., is $\pm\infty$), it is still often the case that $\beta$ has some finite elements that are identified by the model's first-order conditions and that can be consistently estimated. Specifically, as long as separation is “quasi-complete” and there are at least some observations with $x_{i}\overline{\gamma}=0$, coefficients for regressors that do not play a role in the separation can be consistently estimated by first {withholding} any separated observations, and then performing the estimation over the subsample where $x_{i}\overline{\gamma}=0$. Meanwhile, for the parameters that are estimated to be infinite, one can often still estimate finite combinations of these parameters and conduct standard inference on them, {as we show in our proof of part (d) in the Appendix}.\footnote{For example, in a Poisson model where $\mu=\exp\left[\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+\beta_{3}x_{3}\right]$ and $z_{i}=x_{i1}+x_{i2}$ is a linear combination of $x_{1i}$ and $x_{2i}$ that equals $0$ for all $y_{i}=0$ and is $<0$ for some $y_{i}=0$, then $\exp\left[\beta_{0}+\beta_{1}z_{i}+\left(\beta_{2}-\beta_{1}\right)x_{2}+\beta_{3}x_{3}\right]$ is a re-parameterization of $\mu$ that presents the same information about $\beta_{3}$. Here, we know that {the ML estimates for ${\beta}_{1}$ and ${\beta}_{2}$ will both be $+\infty$.} The {ML estimate for the combined parameter $\beta_{2}-\beta_{1}$ is finite}, however, and, {importantly,} the re-parameterized model allows us to take into account its covariance with $\widehat{\beta}_{3}$ in drawing inferences. Moreover, it can still be possible to construct one-sided confidence intervals for the individual parameters $\beta_1$ and $\beta_2$ in cases like these; we provide an example of this type of inference on our \href{https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/README.md}{website}.}
The practical implications of these insights vary based on the model and estimator. For binary choice models, if the data exhibit complete separation instead of only quasi-complete separation, then meaningful estimation is impossible with or without the separated observations. Furthermore, Proposition (ref) is of no use for estimators with potentially unbounded (pseudo-)likelihood functions such as gamma PML, as in these cases the compactified model will have infinitely many solutions when there is separation of any kind. However, as we have discussed, the degree of separation for many other commonly used GLMs can only be quasi-complete. A Poisson model, for example, can always be estimated by first identifying and {withholding the} separated observations from the estimation sample. For these situations, Proposition (ref) lends significant theoretical justification to this approach, especially when the researcher's focus is only on a particular subset of regressors (as is often the case with fixed effect models, for example).
To flesh out some additional intuition behind these results, it is helpful to draw a connection between separation and the better understood result of perfect collinearity between regressors. Under perfect collinearity, there is at least one redundant regressor which, given the other regressors, conveys no additional information about the observed outcomes.\footnote{More precisely, this occurs when $x_{i}\gamma=0$ over the entire sample, for some nonzero vector $\gamma$.} Therefore, the estimated effect of the redundant regressor could theoretically take any value without affecting the score function or the estimates of variables it is not collinear with. Separation is similar in that, because the regressors implicated in $x_{i}\overline{\gamma}$ are only identifiable from the observations where $x_{i}\overline{\gamma}\neq0$, {there must again be at least one regressor that provides no information for estimating the coefficients of the regressors that are not involved in separation.} The two issues are still fundamentally distinct, since separation involves estimates of the problematic regressors becoming infinite rather than indeterminate. In either case, however, it is important that a researcher note that the choice of which regressor to drop from the estimation performed by the estimation algorithm is often arbitrary and that the computed coefficients of some of the remaining regressors (i.e., those that are involved in either separation or perfect collinearity) may need to be interpreted as being relative to an omitted regressor or omitted regressors, as shown in (ref) {in our Appendix}. Indeed, we would generally advise that researchers should be very cautious when a regressor is shown to be dropped by the statistical software they are using, regardless of the underlying cause.
{As another way of illustrating the difference between perfect collinearity and separation, consider the consequences of simply dropping one of the problematic regressors but continuing to use the full sample. When the model exhibits separation, unlike in the case of perfect collinearity, dropping a regressor will have meaningful implications for both the fit of the model and the estimates of all of the model coefficients. Our suggested approach, by contrast, delivers the same model fit one would obtain using maximum likelihood to estimate the full model over the full sample. Since the separated observations are perfectly predicted by the model, their fitted values are obtained as a by-product of the initial step that detects which observations are separated. The fitted values for the remaining observations are then obtained simply by estimating the model on the subsample of non-separated observations. This estimation step also yields the maximum likelihood estimates for all coefficients that have finite estimates in the original model.}\footnote{{ As discussed in the Appendix, it is often possible to recover the signs of the coefficients whose estimates diverge under separation, and it is always possible to estimate certain finite combinations of these coefficients. In practice, recovering the signs of the infinite coefficient estimates requires first identifying the combinations of regressors that separate the data and then carefully re-writing the linear part of the model as in (ref) in our Appendix. Our \href{https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/README.md}{website} includes examples of how to use our ppmlhdfe Stata command to implement these steps. Also see footnote (ref) for simple descriptive example.}}
{As a final remark, another similarity between separation and perfect collinearity is that separation is neither strictly a “small sample” issue nor a “large sample” issue. It may be resolved by obtaining a larger sample if the underlying reason is that there is not enough variation in one or more of the regressors in the current sample. However, it may also occur in large samples either because of fundamental co-dependence between $y_{i}$ and some of the regressors (e.g., a trade embargo may always predict zero exports) or because the number of regressors increases with the sample size (as is typically the case with panel data models and other fixed effects models commonly used in economics research). This motivates our interest in clarifying the large-sample properties of the estimates that are obtained after withholding the separated observations.}
The discussion thus far has been strictly theoretical, but the practical aspects of the separation problem are also interesting. To date, most discussion in the existing literature about how to detect separation has focused on binary outcome models, where the only relevant conditions governing separation are (ref) and (ref). However, for applications with nonbinary outcomes, there will usually be many observations with $0<y_{i}<\overline{y}$, such that the third condition stated in (ref) becomes key. In some cases, this condition can greatly simplify the task of detection. For instance, santos_silva_existence_2010 show that if $X$ is of full column rank over $0<y_{i}<\overline{y}$, then equation (ref) cannot be satisfied and there is no separation. Likewise, if the rank of $X$ over $0<y_{i}<\overline{y}$ is $M-1$, such that there is only one $\gamma^{*}$ that satisfies equation (ref), it is generally easy to compute values for $z_{i}=x_{i}\gamma^{*}$ over the rest of the sample and check whether or not they satisfy the other conditions for separation.
However, detecting separation becomes much more complicated if there are multiple linear combinations of regressors that satisfy equation (ref) (i.e., if $rank(X)<M-1$ over $0<y_{i}<\overline{y}$). Table (ref) gives a simple example of a dataset that presents this issue. In this instance, a check for perfectly collinear regressors over $y_{i}>0$ would quickly reveal that both $z_{1i}=x_{3i}-x_{4i}$ and $z_{2i}=x_{2i}-x_{4i}$ are always $0$ over $y_{i}>0$. The second- and third-to-last columns of Table (ref) then show that both $z_{1i}$ and $z_{2i}$ exhibit overlap over $y_{i}=0$, suggesting that estimates should exist. However, just by virtue of there being two such linear combinations of regressors satisfying (ref), then any other linear combination $z_{3i}$ of the form $z_{3i}=\alpha z_{1i}+(1-\alpha)z_{2i}$ also satisfies (ref). Thus, there are actually an infinite number of linear combinations of regressors one would need to check for overlap in this manner in order to determine existence. In this particular example, it is still possible to determine without too much effort that $z_{3i}=0.5z_{1i}+0.5z_{2i}$ separates the first observation. But for more general cases, a more rigorous approach is needed to take into account the many different ways the data could be separated.
In light of these complexities, silvapulle_existence_1986 and clarkson_computing_1991 have suggested using linear programming methods to detect separation. A suitable illustration of these approaches can be expressed using the following constrained maximization problem:
where $\mathds1_{x_{i}\gamma^{S}<0}$ and $\mathds1_{x_{i}\gamma^{S}>0}$ respectively denote indicator functions for observations with $x_{i}\gamma^{S}<0$ and $x_{i}\gamma^{S}>0$. If a nonzero vector, $\gamma^{S}$, can be found that solves the problem defined by (ref), then the linear combination $x_{i}\gamma^{S}$ clearly satisfies the conditions for separation described in Proposition (ref). {Furthermore, since $\gamma^{S}$ must maximize the number of separated observations, it follows that $\gamma^{S}=\overline{\gamma}$.} A simplex solver or a variety of other similar linear programming methods may be used to solve for $\gamma^{S}$; see konis2007linear for a thorough discussion.
A common weakness of linear programming methods in this context is that they suffer from the curse of dimensionality. Notice that the number of constraints associated with (ref) is equal to the number of observations, $N$, and the number of $\gamma$-parameters that need to be solved for is equal to the number of regressors, $M$. While there are standard operations that may be used to reduce the size of the problem to one with only $N-M$ constraints (cf., konis2007linear, p. 64), an obvious problem nonetheless arises if either $M$ or $N-M$ is a large number, as is increasingly the case in applied economics research.\footnote{As noted in the introduction, this popularity is largely driven by the wide adoption of fixed effects Poisson PML (FE-PPML) estimation for estimating gravity models. For example, figueiredo_industry_2015 estimate a gravity model for patent citations with $N\thickapprox26$ million and $M\thickapprox27,000$, and larch_currency_2017 estimate a similar model for international trade flows with $N\thickapprox880,000$ and $M\approx55,000$. However, high-dimensional fixed effects estimation is also likely to become more attractive for other GLM estimators aside from PPML as well; see stammann2016estimating, stammann2017fast, and fernandez-val_individual_2016 for some relevant innovations that have appeared in the past few years.} In these cases, the standard approach just described necessitates solving a high-dimensional linear programming problem, which may be difficult to solve even using the most computationally efficient linear programming solvers currently available.\footnote{Computationally efficient linear programming solvers typically involve inverting an $M\times M$ basis matrix (cf., hall2011high), a step we would prefer to avoid.} The following discussion, therefore, turns to the question of how to equip researchers to deal with the separation problem in models with many fixed effects and other nuisance parameters.
To introduce a notion of high dimensionality, we will now suppose the set of regressors can be partitioned into two distinct components: a set of $P$ non-fixed effects regressors $w_{i}=w_{1i},\ldots,w_{Pi}$, which we will treat as countable in number, and a set of $Q$ indicator variables $d_{i}=d_{1i},\ldots,d_{Qi}$, where $Q$ is allowed to be a large number. The total number of regressors $M=P+Q$ is therefore also large, and the combined matrix of fixed effect and non-fixed effects regressors can be expressed as $X=\{w_{i},d_{i}\}$. Note that this partition does not depend on the indexing of the fixed effects, but they could easily be subdivided into multiple levels (e.g., “two-way” or “three-way” fixed effects specifications) depending on the application.\footnote{In addition, note that the high-dimensional portion of the regressor set need not consist of only indicator variables; the methods we describe can also be applied to models where $d_{i}$ contains linear time trends, fixed effects interacted with non-fixed effect variables, and so on without loss of generality.} The number of observations, $N$, is assumed to be greater than $M$, with $N-M$ also generally treated as a large number.
Before describing our preferred method for solving this problem, we first briefly discuss the shortcomings of other feasible methods that might otherwise seem appealing.
One strategy is to reduce the dimensionality of the above linear programming problem using the Frisch-Waugh-Lovell theorem, extending an earlier strategy proposed by larch_currency_2017. By projecting out all other regressors---including fixed effects---from each regressor over the sample of positive observations, this reduces the number of parameters to be solved from $M$ to $P$, making the problem much more tractable. This initial projection step can be performed quickly even for very large $Q$ using the methods of {correia_linear_2017. As we discuss further in the Appendix, this approach is effective in many settings but cannot detect separation patterns that involve only fixed effects, since the fixed effects are purged at the outset. While the impact on the estimates of the non-fixed effect parameters from the latter type of separation may be benign, it may still cause numerical issues that slow or prevent convergence of the researcher's estimation algorithm.}
As another alternative, we could simply attempt to compute estimates without any precautions and consider any observation for which the conditional mean appears to be converging numerically to either $0$ or $\overline{y}$ to be separated.\footnote{ppmlhdfe, stammann2017fast, and berge2018efficient each describe algorithms that can accommodate high-dimensional models in a computationally efficient way. The iterative output from these algorithms can in principle be used to detect observations whose computed $\mu$ values are converging to inadmissible values. As discussed in the Appendix, we have made such a method available as an option for our {\tt ppmlhdfe} command.} This strategy has the advantage of being model-agnostic and simple to implement, but it is generally unreliable. As noted by clarkson_computing_1991, it is not guaranteed to detect separation correctly. Moreover, leaving these observations in the regression sample, even temporarily, may again lead to non-convergence and is likely to slow down convergence even in the best of cases. Implementing this method is especially challenging when the true distribution of $\mu$ is very skewed, as it becomes very difficult to numerically distinguish true instances of $\mu=0$ from mere small values of $\mu$.
Our algorithm, which is based on an application of weighted least squares, does not suffer from these types of issues. It can be applied to a very general set of estimation settings, is guaranteed to detect separation, and is both simple to understand and fast. Moreover, it can be implemented in any standard statistical package (without the need for a linear programming solver).
We now turn to describing how the algorithm works for the estimation of Poisson models and similar models with only a lower bound. We will then explain how it may be readily applied to binomial and multinomial models without loss of generality. To proceed, let $u_{i}$ be an artificial regressand such that $u_{i}\le0$ when $y_{i}=0$ and $u_{i}=0$ when $y_{i}>0$. Also, let $\omega_{i}$ be a set of regression weights, given by
with $K$ an arbitrary positive integer. The purpose behind these definitions is that we can choose a sufficiently large $K$ such that a weighted regression of $u_{i}$ on $x_{i}$ can be used to detect if the equality constraint in (ref) can be satisfied by the data. {This result is an application of what is sometimes called the “weighting method” (stewart1997weighting).}\footnote{{It is also similar to penalty methods and barrier methods, which have been used for decades as alternatives to simplex-based algorithms for solving linear programming problems (forsgren2002interior). The value-added of our approach is that it also takes advantage of recent computational innovations that are specific to the estimation of least-squares regressions and that readily accommodate models with arbitrarily high dimensionality.}} {We clarify how this technique works using the following lemma:}
To prove this statement, note first that the residual sum of squares (RSS) minimized by this regression will be at most $u^{\prime}u$. It is then useful to let $K$ equal the smallest integer that is ${>u^{\prime}u/\epsilon^{2}}$. If the weighted least-squares residual $e_{i}$ is greater than $\epsilon$ in absolute magnitude, then that observation will contribute more than ${K\epsilon^{2}}$ to the RSS and the RSS will be at least $K\epsilon^{2}$. If $K>{u^{\prime}u/\epsilon^{2}}$ then RSS$>u^{\prime}u$, which is a contradiction.$\blacksquare$
Because we can force the predicted values of $u_{i}$ from this regression to zero for observations with $y_{i}>0$, the coefficients computed from this regression therefore satisfy (ref). The only remaining step is to choose $u_{i}$ so that all separated observations have predicted values less than zero and all non-separated observations have predicted values equal to zero. We achieve this goal via the following algorithm:
The unweighted $R^{2}$ of the last regression iteration is always equal to $1.0$ when it converges (i.e., $u_{i}=\widehat{u}_{i}$ for all $i$). The following proposition establishes the convergence properties of this algorithm and its effectiveness at detecting separation:
We provide proof of this proposition in our Appendix. The main observation for our current purposes is that none of the above steps are significantly encumbered by the size of the data and/or the complexity of the model. Thanks to the recent innovations of correia_linear_2017, weighted linear regressions with many fixed effects can be computed in almost-linear time (as can more general high-dimensional models using time trends or individual-specific continuous regressors).\footnote{As discussed in guimaraes_simple_2010, this is because we can use the Frisch-Waugh-Lovell theorem to first “partial out” the fixed effects, $d_{i}$, from either side of the problem via a within-transformation operation and then regress the within-transformed residuals of $u_{i}$ on those of the non-fixed effect regressors, $w_{i}$, to obtain $e_{i}$. correia_linear_2017 then shows how to solve the within-transformation sub-problem in nearly linear time.} The above method can therefore be applied to virtually any estimation setting for which (ref) and (ref) are necessary and sufficient conditions for existence, even when the model features many levels of fixed effects and other high-dimensional parameters. Notably, this includes frequency table models\textemdash the original object of interest in haberman_analysis_1974\textemdash which themselves may be thought of as multi-way fixed effects models without non-fixed effect regressors.
The above algorithm still needs a name. Its defining features are that it iteratively uses weighted least squares in combination with a “linear rectifier” function\footnote{We borrow this term from the machine learning literature, where $min(y,0)$ and $max(y,0)$ are known as linear rectifiers or ReLUs (Rectified Linear Units). Despite their simplicity, ReLUs have played a significant role in increasing the accuracy and popularity of deep neural networks glorot2011deep.} to ensure $u_{i}$ eventually converges to the overall certificate of separation that identifies all separated observations. Thus, we have settled on the name “iterative rectifier” (or IR for short).\footnote{“Iteratively Rectified Weighted Least Squares” would have introduced acronym ambiguity with “Iteratively Reweighted Weighted Least Squares”.}
{Finally, it is important to clarify that our iterative rectifier algorithm readily extends to a broader class of models such as binary outcome models and censored models. For the binary outcomes---and, more generally, for multinomial discrete-choice models---the extension amounts to a simple re-parametrization. In particular, a logit model can be rewritten as a Poisson model (see albert_existence_1984, albert_existence_1984).\footnote{albert_existence_1984 conjectured that this type of equivalence between logit and Poisson models could be used to simplify the problem of detecting separation in frequency-table models. The notes we provide in our Appendix include a proof of albert_existence_1984's conjecture.} For larger problems involving binary outcome models, Appendix (ref) explains---and proves---how to use this transformation to write down equivalent Poisson models that are separated if and only if the original binary outcome models are separated. An analogous argument applies for fractional response models where $y_{i}$ can vary continuously over $[0,\overline{y}]$. Lastly, for censored models, as hinted by clarkson_computing_1991 and shown by koll2021, a Type I Tobit model left-censored at zero has the same separation conditions as a Poisson model; and thus our algorithm can be applied directly to this setup.\footnote{See also the \href{https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/nonexistence_examples.md\#tobit-type-i-tobit-model}{tobit example} in our companion website.} } {.67em}
As an illustrative example, we work from the application of baier2019widely. In their main analysis, a “two-stage” method is used to study heterogeneity in the effects of free trade agreements (FTA). In the first stage, a high-dimensional vector of coefficients for each FTA and each pair of countries is estimated using Poisson PML with additional fixed effects. In the second stage, the estimates from the first stage are regressed on a low-dimensional set of covariates in order to examine sources of heterogeneity. The FTA coefficients that are estimated in the first stage differ by the direction of trade, so that, for example, NAFTA-US-Canada and NAFTA-Canada-US are coded as separate indicator variables, each of which switches from 0 to 1 when the NAFTA trade agreement goes into effect in 1994. There are 910 such indicator variables in their data set, which also includes trade flows between 69 countries over the years 1986 to 2006.
The issue we use for illustration arises in the first stage of baier2019widely's procedure. Drawing from recommended practices in the empirical FTA literature, a model close to the one that they estimate is
The dependent variable $y_{ijt}$ is bilateral trade flows, triply indexed for origin country $i$, destination country $j$, and time $t$. The parameters $\xi_{it}$, $\zeta_{jt}$, and $\varphi_{ij}$ therefore are respectively fixed effects for origin-time, destination-time, and origin-destination (or “pair”). These fixed effects alternatively may be thought of as the coefficients of the dummy variables spanning these dimensions, an equivalence we invoke below. The “globalization” coefficient $b_{t}$ measures how much international trade grows relative to each country's domestic sales in each year. It is identifiable in spite of the fixed effects because the data set includes “internal trade” observations capturing sales made by domestic producers in their own markets. The heterogeneous FTA coefficients being estimated are given by $\delta_{A:(i,j)}$, where $A$ denotes a unique FTA and $(i,j)$ serves as an index for each of the origin-destination pairs involved in that FTA. $FTA_{ijt}$ is an indicator equal to 1 when countries $i$ and $j$ are subject to an FTA.
As baier2019widely note, it is not possible to identify a coefficient for the effect of Romania's 1993 FTA with the European Free Trade Area (EFTA) countries on Iceland-Romania trade. No exports from Iceland to Romania were recorded from the beginning of the sample until the FTA begins in 1993. Therefore, all Iceland-Romania observations before 1993 are separated by the following linear combination: \[ -1\times\left(D_{ISL-ROM}-FTA_{ijt}\times D_{ISL-ROM}\right), \] where $D_{ISL-ROM}$ is an indicator equal to 1 for the Iceland-Romania pair. Because the model includes pair fixed effects, $D_{ISL-ROM}$ is effectively one of the regressors. Likewise, $FTA_{ijt}\times D_{ISL-ROM}$ can be regarded as a regressor whose coefficient corresponds to the $\delta_{A:(i,j)}$ parameter for the Iceland-Romania pair upon the signing of the EFTA-Romania FTA. It can be easily verified that this combination induces separation. It is equal to --1 for the $y=0$ observations for Iceland-Romania before the year 1993, satisfying (ref), and is equal to 0 otherwise, satisfying (ref).
Figure (ref) displays the FTA coefficient estimates we obtain when we estimate (ref) using our ppmlhdfe command without any checks for separation. The true coefficient estimate for the Iceland-Romania pair should be infinity. However, though the estimated value for Iceland-Romania is indeed the largest estimate, it does not otherwise stand out as especially problematic given the other extreme values that are found and the overall shape of the distribution. Without checking for separation beforehand, a researcher could easily mistake the reported value for the Iceland-Romania FTA estimate to be legitimate, biasing the subsequent analysis.
This example is well chosen as a use case for our methods because of the high degree of complexity that the model in (ref) embodies. Counting all of the parameters that need to be estimated, there are on the order of 2,800 pair fixed effects, 2,200 exporter-time and importer-time fixed effects, 20 globalization coefficients, and 910 heterogeneous FTA coefficients. When considering possible approaches for detecting separation, the information matrix for this model is not straightforward to obtain or decompose, making it difficult to implement the methods of eck2018computationally. Existing algorithms based on linear programming, such as konis2007linear, are not viable either due to the high dimensionality.\footnote{Moreover, one can envision much larger examples along these lines. Neither the size of the data nor the number of fixed effect parameters are especially large compared with those in, e.g., larch_currency_2017 or french2024effects. } When we apply our iterative rectifier algorithm, which is implemented in ppmlhdfe through the sep(ir) option, the 7 Iceland-Romania observations preceding their FTA are correctly identified as being separated, as are 42 other observations that are perfectly predicted by the pair fixed effects.\footnote{The latter 42 observations that are separated because they are associated with pairs that never trade. Due to the pair fixed effects in the model, all of these observations are perfectly predicted zeroes. Because these cases each only involve a single fixed effect, they are simple to find and can also be detected beforehand using a separate check. For example, one can use the option sep(fe ir) to instruct ppmlhdfe to first check for observations perfectly predicted by any of the fixed effects and then apply the iterative rectifier algorithm to find the remaining 7 separated observations.} The computation takes only one iteration.
To facilitate a comparison of different methods for detecting separation, we next reduce the original data set of baier2019widely to a much smaller one that retains the same structure and the same separation issue involving Iceland-Romania trade. By randomly removing observations, while keeping Iceland-Romania, the reduced version of the data has only 1,176 observations instead of 58,989 and 14 FTA coefficients to be estimated instead of 910. Reducing the data in this way is helpful in part because having a smaller number of coefficients to report allows us to verify numerically that dropping the separated observations does not affect any of the coefficient estimates or their standard errors when done correctly. We also remove beforehand any pairs that never have positive trade. Table (ref) shows results for different approaches and options for tackling separation applied to this reduced data set. Columns 1-3 again show results obtained without any separation checks, only in these cases we experiment with varying the Poisson deviance criterion used to determine if the Poisson PML estimation algorithm has converged. In all cases, an estimate is erroneously reported for the FTA coefficient for Iceland-Romania, but it is interesting to observe how both the computed estimate and its implied statistical significance depend arbitrarily on the chosen tolerance. Column 4 then shows the results for when we apply our iterative rectifier algorithm beforehand. As with the full data set, it performs correctly, dropping the 7 separated observations. Furthermore, none of the other coefficient estimates are affected by dropping these observations as compared to column 1.
The remaining columns of Table (ref) evaluate the performance of the methods of santos_silva_existence_2010. By default, their ppml command checks for separation using two steps. First, it checks for collinearity among all of the regressors over the subsample of positive observations, which is equivalent to checking our condition (ref). Second, if it identifies a regressor as collinear in the first step, it then checks whether the mean value of that regressor over the positive subsample falls between the maximum and minimum values it takes over the sample of $y=0$ observations.\footnote{This checks a version of the relevant “overlap” condition in (ref) but differs because it focuses on whether an individual regressor exhibits overlap rather than on linear combinations of regressors. Importantly, finding in the first step that there is perfect collinearity between regressors over $y>0$ often does not indicate that any specific regressor is responsible for perfect collinearity.} If not, that regressor is dropped. This approach does identify the Iceland-Romania FTA indicator as problematic but is not equipped to detect the role played by the 7 Iceland-Romania observations preceding the FTA and thus does not drop them. Consequently, the results in column 5 for the other coefficient estimates and their standard errors are not the same as in the previous columns. They are still similar in this case due to the high degree of saturation in the model, but it is easy to appreciate how larger differences could arise in more general settings. Finally, column 6 employs the “strict” separation check option for ppml that only checks for collinearity over the $y>0$ sample and does not perform the second step. In this case, ppml now wrongly drops 34 observations that it takes to be perfectly predicted by the dummy variables we use to encode the fixed effects.\footnote{ppml drops observations when the excluded regressor is a dummy variable and when the less common value of the dummy variable occurs when $y=0$. Again, this criterion is not equivalent to our conditions (ref) and (ref) and thus is not able to identify the 7 observations that should be dropped in this case.} Again, we see numerical differences in the estimates as well as their standard errors as a result.
Since our ppmlhdfe Stata package includes several other relatively robust methods for detecting separation, we include in our appendix an expanded version of this example with additional results. We also include on our \href{https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/README.md}{accompanying website} many additional examples demonstrating the implementation of our methods, including for logit models and multinomial logit models. For example, for logit settings, we replicate well known examples from agresti2012categorical,agresti2015foundations, heinze2002solution, and Kosmidis2021detectseparation. For multinomial logit, we replicate the “alligators” example from kosmidis2017multinomial. For Poisson and Poisson PML, we construct 17 of our own examples that researchers working on separation detection methods can use as additional test cases. In addition, we replicate the contingency table example from geyer2009likelihood, including identifying the “direction of recession” in that example. We make both versions of our example based on baier2019widely available as well.
In this paper, we have provided an updated treatment of the concept of separation in the estimation of GLMs. While the result that all GLMs with bounded individual likelihoods suffer from separation under similar circumstances has been shown before by several authors, these results arguably have not received sufficient attention. Now that estimation techniques have progressed to the point where nonlinear models are regularly estimated via (pseudo-)maximum likelihood with many fixed effects, there is considerable ambiguity over whether the estimates produced by these models are likely to exist, what it means when they do not exist, and what can be done to ensure that the model can be successfully estimated.
We have brought more clarity to each of these topics by building on the earlier work of verbeek1989compactification and clarkson_computing_1991, which we have extended to incorporate estimators that have not been previously examined and that have their own more idiosyncratic criteria governing existence. An important takeaway from this analysis is that some, but not all, GLM estimators can still deliver uniquely identified, consistent estimates of at least some of the model parameters even if other parameter estimates are technically infinite.
We have also introduced a new method to detect separation in models with multiple levels of high-dimensional fixed effects, a task that would otherwise require solving an impractical or even infeasible high-dimensional linear programming problem. As GLM estimation with high-dimensional fixed effects increasingly becomes faster and more appealing to researchers, the need for methods that can detect and deal with separation in these models represents an important gap that we aim to fill.
{At the same time, though our methods represent concrete improvements over existing practices for addressing separation in high-dimensional fixed effects settings, they are not a panacea. Fundamentally, detecting separation relies on numerically precise calculations that permit fine distinctions between observations whose predicted values truly lie on the boundary versus those whose predictions differ from it by only a small amount. Verifying these distinctions can be especially difficult for large data sets where the outcome variable is extremely skewed (e.g., sectoral trade data). In these settings, observations that are merely close to the boundary can lead to non-convergence of the estimation algorithm even if the data are not technically separated by the model. For these cases, additional safeguards such as step-halving (marschner2011glm2) may help to prevent erroneous infinite-valued steps before reaching convergence. If possible, it may also be useful to examine predicted values from a preliminary model to identify observations that are likely to lie near the boundary and to assess the consequences of selectively trimming some of these observations as a diagnostic exercise. Formalizing such a procedure, including the consequences for inference, may be a productive avenue for future research.}
\setstretch{1.34} \setstretch{1.35}