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.
71,161 characters · 16 sections · 68 citation commands
Robust Estimation and Inference for Categorical Data
Research in the biomedical, psychological, and social sciences, among others, is often concerned with modeling categorical variables. However, certain errors might be present in the specification of the postulated model or the data themselves. Examples include inattentive responding arias2020,huang2015ier,meade2012 and bot responses ilagan2023 in online Likert-type questionnaire items, misspecification of latent variable models for rating data foldnes2022,foldnes2020polycor, zero-inflated count data lachin2014,lambert1992, and specification errors in models for grouped personal income data victoriafeser1997. These and other studies demonstrate that if unaccounted for, such errors, henceforth collectively called contaminations, can cause model estimates to be severely biased. It is therefore often of interest to obtain contamination-robust model estimates.
Conventional robust estimators huber2009,hampel1986 are generally devised for problems with unbounded sample space and may therefore not be appropriate for models of categorical data, whose sample space may not be bounded and numerical. Hence, the study of robustness for categorical data necessitates non-standard approaches. It turns out that the non-standard nature of this problems yields unexpected results with respect to balancing efficiency and robustness, which is outlined in the following.
We propose a general class of robust estimators for the estimation of models for categorical data, called $C$-estimators (“$C$” for categorical). We show that all estimators in this class, save for one notable exception, are fully efficient at the postulated model, yet they may possess certain robustness properties. This is in contrast to elementary robustness theory according to which there is a fundamental tradeoff between efficiency and robustness huber2009,hampel1986. Another surprising result arises for the aforementioned notable exception: This estimator yields improved first-order robustness but has an unusual asymptotics: At the true model, the estimator does not converge in distribution, but under contamination, it is asymptotically Gaussian, suggesting that inference is easier when the data are contaminated.
As a second contribution, we use the developed limit theory to propose a novel test to formalize and identify categorical “outliers”. In essence, the test assesses if the frequency of a given categorical class can be modeled by the postulated model. If this null hypothesis is rejected, the class is considered “outlying”. We thus not only propose a way to robustly fit models for categorical data, but also a diagnostic tool to identify potential sources of contamination therein.
Motivated by the recent empirical interest in inattentive responding and bot responses arias2020,ilagan2023, we apply the theory developed in this paper to models for questionnaire responses in our companion paper welz2024polycor. Using simulated and empirical data, we find that $C$-estimators are very effective at mitigating bias and identifying potentially inattentive respondents.
This paper is organized as follows. Section (ref) defines the class of $C$-estimators, and Section (ref) describes the relationship between this class and existing estimators. Section (ref) derives the asymptotic properties of the class, and Section (ref) proposes a diagnostic test for detecting categorical outliers. Section (ref) demonstrates the practical usefulness of the developed theory by means of a simulation study on zero-inflated Poisson data. Section (ref) discusses our results and concludes.
This section introduces notation and describes the defining properties of $C$-estimators.
Let $\boldsymbol{Z} = (Z_1,Z_2,\dots, Z_k)^\top$ be a $k$-dimensional categorical, but not necessarily ordinal, random vector taking values in an $m$-dimensional countable sample space $\boldsymbol{\mathcal{Z}}$. Specifically, if $m$ is finite, the sample space is defined as $\boldsymbol{\mathcal{Z}} = \{\boldsymbol{z}_1, \boldsymbol{z}_2, \dots, \boldsymbol{z}_m\}$, and if $m=\infty$, it is defined as the countably infinite set $\boldsymbol{\mathcal{Z}} = \{\boldsymbol{z}_n : n\in\mathbb{N}\}$.
Suppose one postulates a statistical model for vector $\boldsymbol{Z}$ that is subject to a deterministic unobserved $d$-dimensional parameter vector $\bm{\theta}\in\bm{\Theta}\subset\mathbb{R}^d$. Models of categorical variables, denoted $\{\boldsymbol{p} \left( \bm{\theta} \right) = (p_{\boldsymbol{z}}\left(\bm{\theta}\right))_{\boldsymbol{z}\in\boldsymbol{\mathcal{Z}}} : \bm{\theta}\in\bm{\Theta}\}$ here, are characterized by assigning to each discrete outcome $\boldsymbol{z} = (z_1,z_2,\dots,z_k)^\top\in\boldsymbol{\mathcal{Z}}$ a nonnegative probability \[ p_{\boldsymbol{z}}\left(\bm{\theta}\right) = \mathbb{P}_{\bm{\theta}} \left[ \boldsymbol{\mathcal{Z}} = \boldsymbol{z} \right] = \mathbb{P}_{\bm{\theta}} \left[ Z_1 = z_1, \dots, Z_k = z_k \right] \quad\textnormal{such that}\quad \sum_{\boldsymbol{z}\in\boldsymbol{\mathcal{Z}}}p_{\boldsymbol{z}}\left(\bm{\theta}\right) = 1, \] which depends on the parameter $\bm{\theta}\in\bm{\Theta}$. The function $\boldsymbol{z}\mapstop_{\boldsymbol{z}}\left(\bm{\theta}\right)$ is known as probability mass function (PMF).
We now state some minimal assumptions to ensure that the model is well-behaved and allows for statistical inference.
In essence, Assumption (ref) requires the model's PMF to be sufficiently smooth, the model parameters to be point-identifiable, and the model probabilities to be strictly positive. The latter assumption is not strictly required, but makes the notation easier and avoids having to deal with unconstructive technicalities.
In this paper, we focus on situations in which the postulated model is potentially misspecified due to contamination. We adopt the contamination model of ruckstuhl2001, which assumes that observed data of $\boldsymbol{Z}$ is distributed according to a perturbed version of the model. This perturbed model is characterized by a probability function
where the unobserved but fixed $\varepsilon\in[0.5)$ is called the contamination fraction, and $h(\boldsymbol{z})$ is an unspecified arbitrary probability function on $\boldsymbol{\mathcal{Z}}$, called the contamination PMF. Furthermore, $\bm{\theta}_*\in\bm{\Theta}$ in (ref) can be thought of as a true parameter value under which the postulated model generates uncontaminated data, in contrast to $h(\cdot)$, which introduces contamination. Given a nonnegative contamination fraction $\varepsilon>0$, we say that the postulated model is misspecified if $f_{\varepsilon} (\boldsymbol{z}) \neq p_{\boldsymbol{z}}\left(\bm{\theta}\right)$ for any $\boldsymbol{z}\in\boldsymbol{\mathcal{Z}}$ and $\bm{\theta}\in\bm{\Theta}$. Note that the contamination model in (ref) can be seen as an extension of the classical contamination model of huber1964 to categorical data.
Suppose we observe an $N$-sized sample $(\boldsymbol{Z}_i)_{i=1}^N$ of independent copies of $\boldsymbol{Z}$ being distributed according to the possibly contaminated PMF $f_{\varepsilon}$ in (ref). Using this sample, the statistical problem is to estimate the true parameter $\bm{\theta}_*$ with only little bias when the postulated model is possibly misspecified. Denote by \[ \widehat{f}_{N}(\boldsymbol{z}) = \frac{1}{N} \sum_{i=1}^N \mathds{1}\left\{ \boldsymbol{Z}_i = \boldsymbol{z} \right\}, \qquad \boldsymbol{z}\in\boldsymbol{\mathcal{Z}}, \] the empirical probability function, which is a strongly consistent estimator of $f_{\varepsilon} (\boldsymbol{z})$ as $N\to\infty$ vandervaart1998.
$C$-estimation is based on Pearson residuals lindsay1994. For a given parameter $\bm{\theta}\in\bm{\Theta}$, the sample Pearson residual of cell $\boldsymbol{z}\in\boldsymbol{\mathcal{Z}}$, defined as \[ \delta_{\boldsymbol{z},N}\left(\bm{\theta}\right) = \frac{\widehat{f}_{N}(\boldsymbol{z})}{p_{\boldsymbol{z}}\left(\bm{\theta}\right)} - 1, \] measures the discrepancy between the sample probability, $\widehat{f}_{N}(\boldsymbol{z})$, and the model probability, $p_{\boldsymbol{z}}\left(\bm{\theta}\right)$.\footnote{lindsay1994 coins $\delta_{\boldsymbol{z},N}\left(\bm{\theta}\right)$ as “Pearson” residual because Pearson's $\chi^2$-distance can be expressed as a model-weighted sum of squared residuals, namely $\sum_{\boldsymbol{z}\in\boldsymbol{\mathcal{Z}}}p_{\boldsymbol{z}}\left(\bm{\theta}\right)\delta^2_{\boldsymbol{z},N}(\bm{\theta})$.} By construction, Pearson residuals assume values in the interval $[-1,+\infty)$. Values close to 0 indicate good model fit, whereas values away from 0 indicate poor model fit. Indeed, if the model is correctly specified $(\varepsilon = 0)$, then $\widehat{f}_{N}(\boldsymbol{z})\stackrel{\mathrm{a.s.}}{\longrightarrow} f_0(\boldsymbol{z}) = p_{\boldsymbol{z}}\left(\bm{\theta}_*\right)$, so there exists a parameter vector in $\bm{\Theta}$ at which the Pearson residual converges to 0 (in probability), namely the true $\bm{\theta}_*$. Conversely, if the model is misspecified $(\varepsilon > 0)$, then there exists no value in $\bm{\Theta}$ at which the Pearson residual converges (in probability) to 0.
In the same fashion as minimum disparity estimation lindsay1994, the key idea behind $C$-estimation is to downweigh cells whose Pearson residuals are far away from 0, either towards $-1$ or $+\infty$. This is achieved through a prespecified discrepancy function $\rho:[-1:\infty]\to\mathbb{R}$ to map individual Pearson residuals. Applying this mapping to all Pearson residuals gives rise to the empirical risk
with minimizer
The empirical risk (ref) is a special case of the the $f$-divergence of csiszar1963, measuring the divergence between the empirical PMF and the model PMF through the discrepancy function $\rho$. Assumption (ref) collects minimal assumptions on $\rho$.
While discrepancy functions may take negative values, the empirical risk in (ref) as well as its population counterpart are always nonnegative. This follows from Assumption (ref) in combination with equation (35) in csiszar1963.
The fact that Assumption (ref) permits discrepancy functions that are not everywhere twice or thrice differentiable is particularly appealing for robust estimation. Indeed, some of the most popular robust loss functions for continuous variables, such as Huber loss or Tukey's bisquare loss, are defined in piecewise fashion so that the influence of large function arguments (in absolute terms) can be controlled. Such functions may not be twice (or thrice) differentiable at their threshold points.
The discrepancy function $\rho$ gives rise to a specific function, the residual adjustment function lindsay1994, which determines the robustness and efficiency properties of the empirical risk minimizer $\widehat{\bm{\theta}}_N$ in (ref). Denoting $\psi = \rho'$, the RAF $A:[-1,\infty)\to\infty$ is defined as \[ A(x) = (x + 1) \psi(x) - \rho(x) \] with first derivative \[ A'(x) = (x + 1) \psi'(x). \] The RAF naturally arises from the (negative) loss gradient. Indeed, the empirical risk minimizer $\widehat{\bm{\theta}}_N$ in (ref) can equivalently be defined as solution to the estimating equation
where the gradient operator $\nabla_{\bm{\theta}}$ denotes differentiation with respect to $\bm{\theta}$. To ensure that the RAF exhibits meaningful behavior, we impose two minimal assumptions, listed in Assumption (ref).
The assumption that the RAF and its derivative respectively equal 0 and 1 at the origin is without loss of generality for the following reason. Recall that the RAF stems from the loss gradient in (ref). Note that since $\sum_{z\in\boldsymbol{\mathcal{Z}}}p_{\boldsymbol{z}}\left(\bm{\theta}\right) = 1$, we have for any $\bm{\theta}\in\bm{\Theta}$ that $\sum_{\boldsymbol{z}\in\boldsymbol{\mathcal{Z}}}\nabla_{\bm{\theta}}p_{\boldsymbol{z}}\left(\bm{\theta}\right) = \boldsymbol{0}$. Hence, we can replace $A(x)$ by the affine transformation $\widetilde{A}(x) = aA(x) + b$ with $a \neq 0, b\in\mathbb{R}$ being constants, and the solution to the estimating equation would remain unaffected because $\sum_{\boldsymbol{z}\in\boldsymbol{\mathcal{Z}}}\nabla_{\bm{\theta}}p_{\boldsymbol{z}}\left(\bm{\theta}\right) \widetilde{A}(\delta_{\boldsymbol{z},N}\left(\bm{\theta}\right)) \propto \sum_{\boldsymbol{z}\in\boldsymbol{\mathcal{Z}}}\nabla_{\bm{\theta}}p_{\boldsymbol{z}}\left(\bm{\theta}\right) A(\delta_{\boldsymbol{z},N}\left(\bm{\theta}\right))$. Furthermore, because of the possible non-differentiability at certain points (cf. Assumption (ref)), the RAF derivative $A'(0)$ may not exist. The assumption that the RAF is weakly increasing is analogous with a well-known assumption in $M$-estimation, namely that the derivative of the objective function is weakly increasing.
We are now ready to formally define $C$-estimators.
Before we study the theoretical properties of $C$-estimators in Section (ref), we first relate them to existing approaches in the following section.
There exist a number of estimators that have been proposed for the robust estimation of models for categorical data, which are devised for either general or specific models. Specifically, these are $M$-estimators, minimum disparity estimators, general maximum likelihood estimators for grouped data, and $E$-estimators. In the following, we outline how $C$-estimators relates to each of them.
Minimum disparity estimators markatou1997,lindsay1994,simpson1987mde are fully contained in the class of $C$-estimators. MDEs require that the discrepancy function $\rho$ in (ref) is thrice continuously differentiable on $[1, \infty)$. Examples of MDEs are the minimum Hellinger distance estimator, the negative exponential estimator, as well as the the maximum likelihood estimator (MLE). Table (ref) lists their associated discrepancy functions as well as RAFs, and Figure (ref) provides visualizations. lindsay1994 shows that all MDEs have the same influence function and are therefore fully efficient at the true model. Yet, some MDEs possess better robustness properties than the MLE markatou1997,lindsay1994,he1993,simpson1987mde; in particular, lindsay1994 derives a certain breakdown result. However, MDEs are not first-order robust due to having the same influence function as the MLE.
$E$-estimators ruckstuhl2001 are a class of robust estimators of the binomial model. Assuming that $p_z(\theta)$ is the PMF of the binomial model with probability parameter $\theta\in (0,1)$ and sample space $\mathcal{Z} = \{0,1,\dots, m\}$, $E$-estimators minimize the loss in (ref) by using the discrepancy function
where $c_1 < 0 \leq c_2$ are prespecified constants and the convention $0\log(0)=0$ is employed.\footnote{Due to a location shift, Pearson residuals in ruckstuhl2001 take values values in $[0,\infty)$ rather than $[-1,\infty)$. We have adapted their definition of the “Huberized” discrepancy function to our setup where Pearson residuals are defined on $[-1,\infty)$, which is reflected in (ref).} Figure (ref) provides an illustration. The idea behind this discrepancy function is to “Huberize” Pearson residuals in the sense that values below $c_1$ or above $c_2$ will only have a linear effect on the estimate, whereas values in $[c_1,c_2]$ have a superlinear effect. In particular, Pearson residuals falling in $[c_1,c_2]$ behave the same way as in maximum likelihood estimation. As such, this discrepancy function can be seen as a binomial analogue to the huber1964 loss function. However, just like Huber loss, this discrepancy function is not twice (or thrice) differentiable at its threshold points $c_1$ and $c_2$. Hence, $E$-estimators are not minimum disparity estimators since the latter require thrice differentiability. The lack of thrice differentiability is more than just a mere technicality, as it has important consequences for the properties with respect to robustness and asymptotics. ruckstuhl2001 show that as long as $c_2 \neq 0$, the $E$-estimator is fully efficient at the binomial model, but, consequently, is not first-order robust. However, when $c_2 = 0$, the estimator is, in fact, first-order robust, but has a non-Gaussian limit. Nevertheless, the first-order robustness makes the choice $c_2 = 0$ particularly interesting from a robustness perspective. $C$-estimators therefore allow for non-smooth discrepancy functions (Assumption (ref)). Moreover, unlike the class of $E$-estimators as in ruckstuhl2001, $C$-estimators are not restricted to the binomial model, but allows for general models for categorical data. To make this distinction explicit, we call any estimator that uses (ref) as discrepancy function a generalized $E$-estimator.
$C$-estimation does not cover $M$-estimators as well as the class of MGP estimators victoriafeser1997. $M$-estimators for categorical data were originally proposed in hampel1968 for the special cases of the binomial and Poisson model. In our setup of general categorical data, an $M$-estimator is the solution $\bm{\theta} = \widehat{\bm{\theta}}_N$ to the estimating equation
subject to the conditions
where $\boldsymbol{\phi}_b(\boldsymbol{x}) = \boldsymbol{x}\min\{1, b / \| \boldsymbol{x}\|\}, \boldsymbol{x}\in\mathbb{R}^d$,is the multivariate Huber function with tuning constant $b > 0$, $\boldsymbol{s}_{\boldsymbol{z}}\left(\bm{\theta}\right) = \nabla_{\bm{\theta}}\log(p_{\boldsymbol{z}}\left(\bm{\theta}\right))$ is the log-likelihood score function, $\boldsymbol{I}_{d\times d}$ is the $(d\times d)$ identity matrix, $\boldsymbol{H}\left( \bm{\theta} \right)$ is a $(d\times d)$ matrix, and $\boldsymbol{a}\left( \bm{\theta} \right)$ is a $d$-dimensional vector, with the latter two expressed as functions of $\bm{\theta}$ hampel1986. Optimality of this estimator for models for categorical data, among other properties, is proven in simpson1987Mestimator, where “optimal” means the best compromise between efficiency at the true model and first-order robustness.\footnote{More formally, in robust statistics an estimator is called “optimal” if it minimizes its asymptotic variance at the true model subject to a given bound on its influence function at the true model. In other words, an optimal estimator is the most efficient first-order-robust estimator in a certain class of estimators. This notion of optimality is due to hampel1968.}
The class of MGP estimators of victoriafeser1997 was originally devised for grouped data (particularly grouped income data) but can also be applied to general categorical data. For a constant $\gamma\in\mathbb{R}$ and a given function $\boldsymbol{\varphi}_{\boldsymbol{z}}(\bm{\theta}) = \boldsymbol{\varphi}\left( \boldsymbol{z}, \bm{\theta} \right)$ that maps from $\boldsymbol{\mathcal{Z}}\times\bm{\Theta}$ to $\mathbb{R}^d$ and is continuously differentiable with respect to $\bm{\theta}$, an MGP estimator is the solution $\bm{\theta} = \widehat{\bm{\theta}}_N$ to the estimating equation \[ \boldsymbol{0} = \sum_{\boldsymbol{z}\in\boldsymbol{\mathcal{Z}}}\boldsymbol{\varphi}_{\boldsymbol{z}}(\bm{\theta})\widehat{f}_{N}(\boldsymbol{z})^\gamma, \] subject to the Fisher consistency and normalization conditions in (ref). victoriafeser1997 show that the optimal MGP estimator is given by the choice \[ \boldsymbol{\varphi}_{\boldsymbol{z}}(\bm{\theta}) = \boldsymbol{\varphi}_{b,\boldsymbol{z}}(\bm{\theta}) = \boldsymbol{\phi}_b\Big(\boldsymbol{H}\left( \bm{\theta} \right) \big( \boldsymbol{s}_{\boldsymbol{z}}\left(\bm{\theta}\right) - \boldsymbol{a}\left( \bm{\theta} \right) \big)\Big)p_{\boldsymbol{z}}\left(\bm{\theta}\right)^{1-\gamma}, \qquad\boldsymbol{z}\in\boldsymbol{\mathcal{Z}}, \] where matrix $\boldsymbol{H}\left( \bm{\theta} \right)$ and vector $\boldsymbol{a}\left( \bm{\theta} \right)$ are implicitly defined according to (ref). Observe that for the choice $\boldsymbol{\varphi}_{\boldsymbol{z}}(\bm{\theta}) = p_{\boldsymbol{z}}\left(\bm{\theta}\right)^{-\gamma}\nabla_{\bm{\theta}}p_{\boldsymbol{z}}\left(\bm{\theta}\right)$, one can express as MGP estimator the MLE $(\gamma = 1)$ and the minimum Hellinger distance estimator $(\gamma = 0.5)$, which are also contained in the class of $C$-estimators. However, other $C$-estimators, in particular the generalized $E$-estimators, cannot be expressed as MGP estimator.
We proceed by studying the theoretical properties of $C$-estimators.
The estimand of a $C$-estimator $\widehat{\bm{\theta}}_N$ in (ref) is given by the parameter that minimizes the population risk associated with the empirical risk. Formally, for the population Pearson residual and population risk respectively defined by
the estimand equals the population risk minimizer
As such, $\bm{\theta}_0$ depends on the unobserved fraction and type of contamination, $\varepsilon$ and $h(\cdot)$, respectively, as well as discrepancy function $\rho(\cdot)$. Under the additional assumption that $\rho$ is strictly convex, $\bm{\theta}_0$ equals the true $\bm{\theta}_*$ in the absence of model misspecification $(\varepsilon = 0)$ csiszar1963. In other words, if the discrepancy function is strictly convex, then $\widehat{\bm{\theta}}_N$ is Fisher-consistent. Assumption (ref) establishes point-identification of $\bm{\theta}_0$ under compactness of the parameter space.
We proceed by studying the estimator's limit behavior.
The following theorem establishes consistency of $\widehat{\bm{\theta}}_N$ for estimand $\bm{\theta}_0$.
The proofs of this theorem and all subsequent mathematical statements in this paper are deferred to Appendix (ref).
We now turn to the asymptotic distribution of $C$-estimators, for which we introduce additional notation to accommodate two special cases.
First, we will frequently perform the operation $\boldsymbol{A}\boldsymbol{b}$ for a $(d\times m)$-dimensional matrix $\boldsymbol{A}$ and an $m$-dimensional vector $\boldsymbol{b}$. In the case where $m=\infty$ (like the Poisson model), we define the matrix operation $\boldsymbol{Ab}$ based on a certain inner product on an infinite-dimensional vector space. Specifically, in our setup, that inner product is given by an infinite sum over the inner product of the corresponding row in $\boldsymbol{A}$ and vector $\boldsymbol{b}$. With this definition, we have that \[ \boldsymbol{A}\boldsymbol{b} =
, \] so that the resulting vector is of finite dimension $d$. Similar reasoning applies to the operation $\boldsymbol{A}\boldsymbol{B}\boldsymbol{A}^\top$, where $\boldsymbol{B}$ is a $(m \times m)$ matrix. Throughout this paper, whenever we sum over the sample space $\boldsymbol{\mathcal{Z}}$ like above, a PMF of a discrete distribution is involved such that the summation is equal to an expectation over that distribution. We will later impose an assumption that ensures that this expectation exists and is finite when the sample space $\boldsymbol{\mathcal{Z}}$ is infinite and (possibly) $\varepsilon > 0$.
Second, in some objects relevant to the estimator's asymptotics, we often need to evaluate limits of the form
where $\delta : \bm{\Theta}\to [-1,\infty)$ denotes a generic Pearson residual, $\alpha > 0$ is a scalar, $\boldsymbol{w}\neq\boldsymbol{0}$ is a deterministic $d$-dimensional vector, and $\bm{\theta}\in\bm{\Theta}$ is arbitrary. If $\psi$ is differentiable at $\delta(\bm{\theta})$, then this limit equals $\psi'(\delta(\bm{\theta}))$. However, $\psi$ may not be differentiable at $\delta(\bm{\theta})$ (cf. Assumption (ref)). It is useful for the asymptotic analysis to introduce additional notation that reflects this possible non-differentiability.
In the univariate case ($d=1$), there are two directions along which one can evaluate the limit in (ref), governed by the sign of $w\neq 0$. If $w>0$, then (ref) equals the right limit $\lim_{x\downarrow \theta}\psi'(\delta(x)) =: \psi'(\delta(\theta+)) $, and if $w<0$, then it equals the left limit $\lim_{x\uparrow \theta}\psi'(\delta(x)) =: \psi'(\delta(\theta-))$. We now describe a notion of such directional limits that also covers the multidimensional case ($d>1$).
To start, note that for evaluating the limit in (ref), the exact value of the $d$-dimensional vector $\boldsymbol{w}$ is not important. Rather than $\boldsymbol{w}$ itself, the sign of each coordinate in $\boldsymbol{w} = (w_1,\dots, w_d)^\top$ determines the direction of the linear path $\bm{\theta} + \alpha \boldsymbol{w}$ as $\alpha\downarrow 0$. Defining the coordinatewise sign function $\bm{\kappa}: \mathbb{R}^d\to\{-1,0,1\}^d$ as \[ \bm{\kappa}(\boldsymbol{w}) = \Big( \textnormal{sign}\big(w_j\big) \Big)_{j=1}^d \] enables us to, without loss of generality, replace $\boldsymbol{w}$ by $\bm{\kappa}(\boldsymbol{w})$ in limit (ref). Since $\bm{\kappa}(\boldsymbol{w})\in\{-1,0,1\}^d\setminus\{\boldsymbol{0}\}$ for any $\boldsymbol{w}\neq \boldsymbol{0}$, this notation makes it explicit that there are $3^d-1$ possible directions along which one could evaluate the limit along the linear path, and $\boldsymbol{w}$ uniquely determines the choice of direction.\footnote{In addition to multiple directions of a linear path, there are infinitely many (possibly nonlinear) paths along which one could evaluate the limit of a multivariable function, but, in the asymptotic analysis it suffices to only consider the linear path in (ref) for the analysis of the function composition $\psi'\circ \delta$.}
With the coordinatewise sign function, we can define the shorthand
which equals the limit in (ref) for any $\boldsymbol{w}\neq\boldsymbol{0}$. The superscript $``(\bm{\kappa}(\boldsymbol{w}))\pm"$ reminds us that we do not evaluate the function composition $\psi'\circ \delta$ at $\bm{\theta}$, but along a linear path $\bm{\theta} + \alpha\bm{\kappa} \left( \boldsymbol{w}\right))$, $\alpha\downarrow 0$, whose direction is governed by the coordinatewise signs $\bm{\kappa} \left( \boldsymbol{w}\right)$. As such, the definition $\psi'\left( {\delta}\left({\bm{\theta}^{\left(\bm{\kappa} \left( \boldsymbol{w}\right) \right)\pm}}\right) \right)$ generalizes the well-known unidimensional left and right limits to higher dimensions. Indeed, if $d=1$, then $\psi'(\delta(\theta^{(\kappa(w))\pm })) = \lim_{x\uparrow \theta} \psi'(\delta(x)) = \psi'(\delta(\theta -))$ if $\kappa(w) = -1$, and $\psi'(\delta(\theta^{(\kappa(w))\pm })) = \lim_{x\downarrow \theta} \psi'(\delta(x)) = \psi'(\delta(\theta +))$ if $\kappa(w) = 1$, for any unidimensional $\theta$ and $w\neq 0$. Furthermore, if ${\psi'}\left({{\delta}\left({\bm{\theta}}\right)}\right)}\right)}$ exists, then $\psi'\left( {\delta}\left({\bm{\theta}^{\left(\bm{\kappa} \left( \boldsymbol{w}\right) \right)\pm}}\right) \right) = {\psi'}\left({{\delta}\left({\bm{\theta}}\right)}\right)}\right)}$ for all directions $\boldsymbol{w}$.
With the new notation, we are ready to define matrices that are important in the estimator's limit theory. Specifically, the asymptotic covariance matrix of $C$-estimators is a function of only two matrices, $\boldsymbol{M}$ and $\boldsymbol{U}$, which will be defined in the following.
For $\bm{\theta}\in\bm{\Theta}$, put the $(d\times m)$ matrix \[ \boldsymbol{W}(\bm{\theta})= \bigg( {A'}\left({\delta_{\boldsymbol{z}_1,\varepsilon}\left(\bm{\theta}\right)}\right)\boldsymbol{s}_{\boldsymbol{z}_{1}}\left(\bm{\theta}\right), \cdots, {A'}\left({\delta_{\boldsymbol{z}_m,\varepsilon}\left(\bm{\theta}\right)}\right)\boldsymbol{s}_{\boldsymbol{z}_{m}}\left(\bm{\theta}\right) \bigg), \] where $\boldsymbol{s}_{\boldsymbol{z}}\left(\bm{\theta}\right) = \nabla_{\bm{\theta}}\log(p_{\boldsymbol{z}}\left(\bm{\theta}\right)) = \nabla_{\bm{\theta}}p_{\boldsymbol{z}}\left(\bm{\theta}\right)/p_{\boldsymbol{z}}\left(\bm{\theta}\right)$ is the log-likelihood score function, and, for $\boldsymbol{f}_\varepsilon = (f_{\varepsilon} (\boldsymbol{z}))_{\boldsymbol{z}\in\boldsymbol{\mathcal{Z}}}$, define the $(m\times m)$ matrix \[ \boldsymbol{\Omega} = \textnormal{diag} \left( \boldsymbol{f}_\varepsilon\right) - \boldsymbol{f}_\varepsilon\boldsymbol{f}_\varepsilon^\top. \] The sole purpose of matrices $\boldsymbol{W}(\bm{\theta})$ and $\boldsymbol{\Omega}$ is to calculate the symmetric $(d\times d)$ matrix \[ \boldsymbol{U}\left( \bm{\theta} \right) = \boldsymbol{W}\left( \bm{\theta} \right) \boldsymbol{\Omega} \boldsymbol{W}\left( \bm{\theta} \right)^\top. \] Furthermore, for $\bm{\theta},\bm{\theta}'\in\bm{\Theta}$, define the symmetric $(d\times d)$ matrix \[ \boldsymbol{M}\left( \bm{\theta},\bm{\theta}' \right) = \sum_{\boldsymbol{z}\in\boldsymbol{\mathcal{Z}}} \bigg[ \boldsymbol{s}_{\boldsymbol{z}}\left(\bm{\theta}\right) \boldsymbol{s}_{\boldsymbol{z}}\left(\bm{\theta}\right)^\top f_{\varepsilon} (\boldsymbol{z}) \psi' \left( \delta_{\boldsymbol{z},\varepsilon}\left(\bm{\theta}'\right)\right)(\delta_{\boldsymbol{z},\varepsilon}\left(\bm{\theta}\right)+1) - \left(\nabla^2_{\bm{\theta}}p_{\boldsymbol{z}}\left(\bm{\theta}\right)\right) A \left( \delta_{\boldsymbol{z},\varepsilon}\left(\bm{\theta}\right)\right) \bigg] \] and put \[ \boldsymbol{M}\left( \bm{\theta} \right) = \boldsymbol{M}\left( \bm{\theta},\bm{\theta} \right) = \sum_{\boldsymbol{z}\in\boldsymbol{\mathcal{Z}}} \bigg[ \boldsymbol{s}_{\boldsymbol{z}}\left(\bm{\theta}\right) \boldsymbol{s}_{\boldsymbol{z}}\left(\bm{\theta}\right)^\top f_{\varepsilon} (\boldsymbol{z}) {A'}\left({\delta_{\boldsymbol{z},\varepsilon}\left(\bm{\theta}\right)}\right) - \left(\nabla^2_{\bm{\theta}}p_{\boldsymbol{z}}\left(\bm{\theta}\right)\right) A \left( \delta_{\boldsymbol{z},\varepsilon}\left(\bm{\theta}\right)\right) \bigg], \] where $\nabla^2_{\bm{\theta}}p_{\boldsymbol{z}}\left(\bm{\theta}\right) = \frac{\partial^2 }{\partial {\bm{\theta}}\partial {\bm{\theta}}^\top}p_{\boldsymbol{z}}\left(\bm{\theta}\right)$ denotes the Hessian matrix of $p_{\boldsymbol{z}}\left(\bm{\theta}\right)$.
There are two situations warranting special attention.
First, both $\boldsymbol{M}\left( \bm{\theta},\bm{\theta}' \right)$ and $\boldsymbol{U}\left( \bm{\theta} \right)$ depend on ${\psi'}\left({\delta_{\boldsymbol{z},\varepsilon}\left(\bm{\theta}\right)}\right)$ (the latter through $x\mapsto A'(x) = (x+1)\psi'(x)$), but $\psi'$ may not exist everywhere (cf. Assumption (ref)). To circumvent this issue, we sometimes use for a given direction $\boldsymbol{w} \neq \boldsymbol{0}$ the limit approximations $\boldsymbol{M}\left( \bm{\theta}, \bm{\theta}^{\left(\bm{\kappa} \left( \boldsymbol{w}\right) \right)\pm} \right)$ and $\boldsymbol{U}\left( \bm{\theta}^{\left(\bm{\kappa} \left( \boldsymbol{w}\right) \right)\pm} \right)$, which are to be understood in the sense of (ref).
Second, the matrix-valued functions $\boldsymbol{M}$ and $\boldsymbol{U}$ may not necessarily be finite across all possible contamination fractions ($\varepsilon > 0$) and types (PMF $h(\cdot)$). This is particularly a concern when the sample space is of infinite cardinality ($m = \infty$) since we evaluate infinite sums involving the population PMF $f_{\varepsilon}$. Hence, we impose a finiteness and invertibility assumption, similar to lindsay1994 and victoriafeser1997.
Observe that if ${\psi'}\left({\delta_{\boldsymbol{z},\varepsilon}\left(\bm{\theta}_0\right)}\right)$ exists for all $\boldsymbol{z}\in\boldsymbol{\mathcal{Z}}$, then this assumption requires that the matrices $\boldsymbol{M}\left( \bm{\theta}_0 \right)$ and $\boldsymbol{U}\left( \bm{\theta}_0 \right)$ exist and are finite as well as positive definite. A necessary condition for this assumption to be satisfied is that all ${\psi'}\left({\delta_{\boldsymbol{z},\varepsilon}\left(\bm{\theta}_0\right)}\right)$ are strictly positive if they exist.
The following assumption, which is our final one, imposes a certain notion of local convexity of population risk $L_\varepsilon (\bm{\theta})$ in a neighborhood of estimand $\bm{\theta}_0$.
Requirements like Assumptions (ref) on the local behavior of the population risk in a neighborhood of the estimand are common when studying the limit distribution of empirical risk minimizers. For instance, vandervaart1998 assumes a variant of Lipschitz continuity in a neighborhood of the population risk minimizer, lindsay1994 requires certain functions of the parameters to be locally dominated, and ruckstuhl2001 impose local convexity of the (in their case univariate) population risk.
We proceed with an intuitive explanation of a certain matrix that plays an important role in the estimator's limit distribution by uniquely determining its asymptotic covariance matrix (despite not being the latter's square root matrix).
Fix a $d$-dimensional $\boldsymbol{t}\neq \boldsymbol{0}$ and let $\boldsymbol{V}$ be a $(d\times d)$ matrix that satisfies the equation
We will later show that although $\boldsymbol{V}$ may not be symmetric, all of its eigenvalues are real and positive definite under our assumptions. To provide intuition for this expression and the matrix $\boldsymbol{V}$, consider first the unidimensional case $(d=1)$ so that $t\neq 0$ and $V>0$ are scalars. In this case, the direction along which the limit in (ref) is evaluated is fully determined by the sign of $t$ because $V$ is always positive, that is, $\kappa(Vt) = \kappa(t) = \textnormal{sign}(t)$. Indeed, if $t>0$, then equation (ref) is satisfied by the right-limit solution $V = \sqrt{{U}\left({\theta_0+}\right)} \big/ {M}\left({\theta_0,\theta_0+}\right)$, and if $t <0$ it is satisfied by the left-limit solution $V = \sqrt{{U}\left({\theta_0-}\right)} \big/ {M}\left({\theta_0,\theta_0-}\right)$. It follows that the solution $V$ depends on the sign of $t$. For the special case of the (unidimensional) binomial model, such a scalar $V>0$ that depends on the sign of $t$ and satisfies (ref) plays a crucial role in the asymptotic theory of ruckstuhl2001.
However, in the multidimensional case $(d > 1)$, finding a matrix $\boldsymbol{V}$ that satisfies equation (ref) is more involved because the direction along which the limit is evaluated not only depends on vector $\boldsymbol{t}$, but also on matrix $\boldsymbol{V}$ itself. Indeed, the $j$-th direction in the coordinatewise sign function $\bm{\kappa} \left( \boldsymbol{Vt}\right)$ is equal to the sign of the $j$-th coordinate of vector $\boldsymbol{Vt}$, $j=1,\dots,d$. Thus, in order to find the desired matrix $\boldsymbol{V}$, we need to find the direction $\boldsymbol{w}\in\{-1,0,1\}^d\setminus{\{\boldsymbol{0}\}}$ for which the candidate solution
satisfies (ref). Importantly, this solution depends on $\boldsymbol{t}$ and may change for different $\boldsymbol{t}$. Thus, to conclude, given $\boldsymbol{t}\neq\boldsymbol{0}$, we know the functional form of the solution matrix $\boldsymbol{V}$ (see previous display), but we do not a priori know the direction $\boldsymbol{w} = \bm{\kappa} \left( \boldsymbol{Vt}\right)$ associated with the solution matrix $\boldsymbol{V}$ satisfying (ref). Since there is only a finite number of directions here, namely $3^d-1$, we can for a given $\boldsymbol{t}$ calculate candidate solutions in (ref) for all feasible directions $\boldsymbol{w}$ and then evaluate which one satisfies the objective equation (ref). This solution exists under our assumptions because Assumption (ref) requires all candidate solutions in (ref) to exist. A later proposition will show that it is furthermore unique and all of its eigenvalues are real and strictly positive, and we will formalize the process of finding this solution by means of an algorithm. We stress that in dimension $d>1$, matrix $\boldsymbol{V}$ in (ref) need not be symmetric positive definite. As such, it is not equal to a square root matrix of a $C$-estimator's asymptotic covariance matrix, but is nevertheless useful for constructing said asymptotic covariance matrix.
We are now ready to state the asymptotic distribution of $C$-estimators.
An algorithm to find the implicitly defined matrix $\boldsymbol{V_t}$ in (ref) is described in Appendix (ref). The following proposition establishes the properties of this matrix.
The following corollary is an immediate consequence of Theorem (ref).
A consistent estimator of the unobserved asymptotic covariance matrix $\bm{\Sigma}\left({\bm{\theta}_0}\right)$ can be constructed as follows. Replace all population class probabilities $f_{\varepsilon} (\boldsymbol{z})$ by their corresponding empirical counterparts $\widehat{f}_{N}(\boldsymbol{z})$ in matrices $\boldsymbol{W}\left( \bm{\theta} \right), \boldsymbol{M}\left( \bm{\theta} \right)$, and $\boldsymbol{\Omega}$. Then exploit the plug-in principle and evaluate $\boldsymbol{U}\left( \bm{\theta} \right)$ and $\boldsymbol{M}\left( \bm{\theta} \right)$ at the point estimate $\widehat{\bm{\theta}}_N.$ Denote the ensuing plug-in estimator by $\bm{\Sigma}\left({\widehat{\bm{\theta}}_N}\right)$, which is consistent for $\bm{\Sigma}\left({\bm{\theta}_0}\right)$ by Theorem (ref) and the continuous mapping theorem.
It turns out that if the model is correctly specified $(\varepsilon = 0)$ and $\psi'(0)$ exists (which it does if the discrepancy function is twice differentiable), then the asymptotic covariance matrix of all $C$-estimators is equal to the inverse Fisher matrix at the model. We formalize this result in Lemma (ref) in the appendix. Hence, all $C$-estimators with a twice differentiable discrepancy function are fully efficient, which suggests that robustness can be achieved without having to sacrifice efficiency. To study this apparent lack of a robustness-efficiency tradeoff, the following section is concerned with influence functions.
The maximum likelihood estimator of the postulated model $\boldsymbol{p} \left( \bm{\theta} \right)$ is obtained by choosing $x\mapsto \rho(x) = (x+1)\log(x+1)$ in the empirical risk (ref). The influence function hampel1974 of a maximum likelihood estimator $\widehat{\bm{\theta}}_N^{\mathrm{\ MLE}}$ is given by \[ \mathrm{IF}\left(\boldsymbol{z}, \widehat{\bm{\theta}}_N^{\mathrm{\ MLE}}, \boldsymbol{p} \left( \bm{\theta}_* \right)\right) = \boldsymbol{J}\left( \bm{\theta}_* \right)^{-1}\boldsymbol{s}_{\boldsymbol{z}}\left(\bm{\theta}_*\right), \qquad\boldsymbol{z}\in\boldsymbol{\mathcal{Z}}, \] where $\boldsymbol{J}\left( \bm{\theta} \right) = \sum_{\boldsymbol{z}\in\boldsymbol{\mathcal{Z}}}p_{\boldsymbol{z}}\left(\bm{\theta}\right)\boldsymbol{s}_{\boldsymbol{z}}\left(\bm{\theta}\right)\boldsymbol{s}_{\boldsymbol{z}}\left(\bm{\theta}\right)^\top$ denotes the model's Fisher informatiom matrix. The influence function describes the bias caused by an infinitesimally small contamination on an estimate markatou1997.
The following theorem states the influence function $C$-estimators.
We close this section by discussing the derived theoretical properties.
If $\psi'(0)$ exists, then all $C$-estimators have the influence function of the MLE and are therefore fully efficient at the postulated model. On the other hand, an influence function equal to that of the MLE means a lack of first-order robustness. However, many choices of the discrepancy functions are designed to downweight the influence of poorly fitted classes (see Table (ref)), so one would expect a robustness gain over the MLE. In contrast to elementary robustness theory hampel1986, this suggests that the influence function does not carry all information on a procedure's robustness properties in the case of categorical data. Hence, when modeling categorical data, it is possible to construct fully efficient and robust estimators (cf. Theorem (ref)).
If $\psi'(0)$ does not exist, then it is possible to construct first-order robust estimators, that is, estimators with a smaller influence function than the MLE. For instance, this applies to estimators with the discrepancy function of ruckstuhl2001 in (ref), for the choice $c_2 = 0$. While first-order robustness is in principle appealing, it comes at the significant downside that $C$-estimators with non-existing $\psi'(0)$ do not converge in distribution at the postulated model (Theorem (ref)). Curiously, under contamination the same estimator does converge in distribution, and is asymptotically Gaussian. Therefore, inference is easier with contaminated samples than non-contaminated ones. This result ties in with ruckstuhl2001, who demonstrated it for the binomial model, and we have shown that it also holds true for general models of categorical data. In general, non-standard asymptotics for estimators with loss functions that are not everywhere twice differentiable have been described before in huber1964, pollard1985, simpson1987Mestimator, and huber2009 for location estimation.
In addition to robust estimation, it is often of interest to identify outliers in a given sample. It may not be intuitive what constitutes outlyingness in categorical data, which may only be supported on a finite set and/or non-numerical. We propose to adopt the philosophy of croux2007 and define a categorical outlier as “an observation which is unlikely to have been generated by the imposed model”. Equivalent definitions have been proposed in davies1993, lindsay1994, and markatou1997.
With this outlier philosophy, one may identify a categorical cell $\boldsymbol{z}\in\boldsymbol{\mathcal{Z}}$ as outlying if its population probability $f_{\varepsilon} (\boldsymbol{z})$ is statistically significantly different from the model probability at estimand $\bm{\theta}_0$, $p_{\boldsymbol{z}}\left(\bm{\theta}_0\right)$, which corresponds to the best fit that can be achieved with the postulated model for the chosen discrepancy function. In other words, for a given cell $\boldsymbol{z}$, we wish to perform the hypothesis test
where the alternative hypothesis may be replaced by a one-sided alternative depending on whether one expects outlyingness to manifest in $p_{\boldsymbol{z}}\left(\bm{\theta}_0\right)$ being larger or smaller than $f_{\varepsilon} (\boldsymbol{z})$. Such a priori beliefs are implicitly imposed by the choice of discrepancy function. While the discrepancy functions in Table (ref) are primarily designed to downweight classes with large Pearson residuals, there exist discrepancy functions that can downweight “inliers”, that is, classes whose Pearson residuals are structurally smaller than 0. Examples are the negative exponential function (see Table (ref), proposed by lindsay1994) and the discrepancy function in eq. (ref) with the choice $c_1 = 0$ (see ruckstuhl2001 for a discussion).
Ideally, a test for the hypotheses in (ref) will reject $\textnormal{H}_{0}$ if cell $\boldsymbol{z}$ is outlying, and sustain $\textnormal{H}_{0}$ if it is not. It turns out that a test statistic that satisfies these two desirable properties is given by
where $\sigma_{\boldsymbol{z}}^2\left(\bm{\theta}_0\right) = \nabla_{\bm{\theta}}p_{\boldsymbol{z}}\left(\bm{\theta}_0\right)^\top \bm{\Sigma}\left({\bm{\theta}_0}\right)\nabla_{\bm{\theta}}p_{\boldsymbol{z}}\left(\bm{\theta}_0\right)$. The following corollary of Theorem (ref) establishes the validity of the test statistic $T_N(\boldsymbol{z})$ for testing the null hypothesis $\textnormal{H}_{0}: \delta_{\boldsymbol{z},\varepsilon}\left(\bm{\theta}_0\right)=0$.
In practice, the test statistic is unobserved because neither $\sigma_{\boldsymbol{z}}^2\left(\bm{\theta}_0\right)$ nor $f_{\varepsilon} (\boldsymbol{z})$ are observed. However, these two quantities can be consistently estimated through $\sigma_{\boldsymbol{z}}^2\left(\widehat{\bm{\theta}}_N\right)$ (due to Theorem (ref) and the continuous mapping theorem) and $\widehat{f}_{N}(\boldsymbol{z})$, respectively. Hence, in practice, the feasible test statistic
which is consistent for $T_N(\boldsymbol{z})$ by the continuous mapping theorem, may used for testing the null hypothesis in (ref).
The simulation study in this section considers $C$-estimation for the problem of fitting the Poisson model with zero-inflated count data. The inflation of zeros is frequently encountered in count data in biomedical sciences, such as patient days in a hospital, the number of wisdom teeth extracted, or the number of episodes of hypoglycemia (low blood sugar) per year for diabetes patients lachin2014.
Denote by $p_{z}\left(\lambda\right) = \lambda^{z}\exp(-\lambda)/z!$ the Poisson PMF with sample space $\mathcal{Z} = \{0,1,2,\dots\}\ni z$ and rate parameter $\lambda > 0$. We draw a sample $Z_1,\dots, Z_N$ of size $N=1,000$ from the Poisson PMF with true rate parameter $\lambda_* = 4.2$. Then, we introduce contamination by randomly replacing a fraction $\varepsilon$ of the sample by zeros, where we consider $\varepsilon\in\{0,0.1,0.2, 0.3\}$. We generate 1,000 datasets with this process. A particular challenge for estimation is overlap between the contaminated observations and model-generated observations since zeros have a strictly positive probability under the model at the true parameter, $p_{0}\left(\lambda_*\right) \approx 0.015$.
Using the contaminated datasets, we fit the Poisson model with five types of $C$-estimators, namely maximum likelihood, Hellinger distance estimation, negative exponential estimation, and the generalized $E$-estimator (GE) (see Table (ref)). For the GE, we choose two configurations of tuning constants, namely $c_1 = -\infty$ and $c = c_2 = 0.6$ as well as $c_1 = -\infty$ and $c = c_2 = 0$ so that we have one case where the estimator is asymptotically Gaussian at the Poisson model (the former), and one where it is not (the latter). The choice $c_1= -\infty$ follows a suggestion by ruckstuhl2001.
Figure (ref) visualizes the simulation results by means of boxplots of the bias $\widehat{\lambda}_N-\lambda_*$ of each estimator for every considered contamination fraction.\footnote{It is worth mentioning that the whole simulation study executed in less than one second on a standard laptop with 32GB RAM and 22 Ultra 7 165H CPUs, running on Ubuntu 22.}
In the no-contamination case ($\varepsilon = 0$), all estimators (except the GE for $c=0$) accurately estimate the true $\lambda_*$ with approximately equal variance. This result is expected because these estimators share the same influence function (Theorem (ref)). The GE for $c=0$ seems to exhibit a small bias at the true model, though. We explain in Appendix (ref) that this is a simple finite sample issue that vanishes in larger sample sizes.
When contamination is present ($\varepsilon > 0$), the estimators deviate noticeably. The MLE is substantially biased already at $\varepsilon = 0.1$, and further deteriorates at increasing contamination levels. The Hellinger distance estimator is an improvement to the MLE in terms of robustness, but is still considerably biased, especially at high contamination levels (e.g., bias of about $-0.7$ at $\varepsilon = 0.3$). The negative exponential is more robust, but still more biased than the two generalized $E$-estimators (GEs), which display excellent robustness properties. In particular, the GE with $c=0$ exhibits almost no bias even at the high contamination level of $\varepsilon = 0.3$. Also the second GE estimator with the less robust choice $c=0.6$ is very robust and only slightly more biased than with $c=0$.
We extend this simulation to the diagnostic test of the previous section in Appendix (ref). In summary, the test has excellent power for the detection of outliers, but its feasible version $\widehat{T}_N(\boldsymbol{z})$ occasionally overrejects at the true model, which is likely due to the additional estimation uncertainty stemming from the use of estimate $\widehat{f}_{N}(\boldsymbol{z})$. We discuss these results in more detail in Appendix (ref).
While the minimum disparity estimators (Hellinger and negative exponential) offer enhanced robustness compared to MLE, we conclude that the generalized $E$-estimators are most robust in this simulation. In fact, even at 30% contamination, the two considered generalized $E$-estimators have almost no bias. Thus, robust $C$-estimators, particularly generalized $E$-estimators, could be an alternative to existing approaches to zero-inflated count data. Existing approaches, such as the estimator of lambert1992, generally model the inflation of zeros explicitly through mixtures of two Poisson distributions. In contrast, $C$-estimators do not explicitly model contamination. As such, they are designed to be robust against arbitrary specification errors of the Poisson model, which may go beyond inflated zero counts.
We have proposed a new class of robust estimators for general models of categorical data, called $C$-estimators. $C$-estimators extend the class of minimum disparity estimators of lindsay1994 by allowing for risk functions that are not everywhere twice differentiable, which can result in more robust estimators ruckstuhl2001. It turns out that such risk functions lead to enhanced first-order robustness if the risk is not twice differentiable at the origin. However, $C$-estimators with such a risk function are not asymptotically Gaussian at the postulated model, but, quite surprisingly, are asymptotically Gaussian when contamination is present. In contrast, $C$-estimators whose risk is twice differentiable at the origin are asymptotically Gaussian and fully efficient at the model, suggesting that robustness can be achieved without sacrificing efficiency. In addition, we propose a diagnostic test to identify categorical outliers.
A particularly relevant area of further research are robustness measures for estimators of models for categorical data, like $C$-estimators. The fact that robust estimators can have the same influence function as non-robust maximum likelihood suggests that the influence function may not be an appropriate robustness measure. This aspect is discussed in lindsay1994, who argues that the second derivative of the residual adjustment function---a crucial function in estimation with categorical data---carries information on the robustness of such estimators. lindsay1994 therefore suggests using the second derivative of the residual adjustment function evaluated at 0 as robustness measure, with large negative values indicating greater robustness. However, lindsay1994 cautions in Remark E that this criterion, while sufficient, may not be a necessary condition for robustness. In addition, this criterion does not exist for first-order robust estimators, whose risk is not twice differentiable at 0. Hence, robustness measures that allow for the comparison of different $C$-estimators are needed. Potentially fruitful approaches could be the local shift sensitivity of hampel1974 or angular breakdown values zhao2018.
From a practical perspective, the theory developed in this paper might be particularly useful for robustifying the analyses of questionnaire response data against inattentive responding or bot responses, which is an increasing concern in internet-collected data ilagan2023,arias2020,meade2012. We apply our theory to models of response data in the companion paper welz2024polycor and demonstrate that it yields substantial improvements in terms of robustness against erroneous responses while also helping identify them. Hence, $C$-estimators could aid in making the analysis of categorical data, particularly rating data, less dependent on model assumptions.