Robust estimators for count data

Question

I am looking for robust estimators for parameters of a processes producing count data: $$ (n_1,...,n_K), n_i\in\mathbb{N} $$ that is the underlying distribution is something like Poisson or Negative binomial, but possibly contaminated - with some samples exhibiting an excessively low or excessively high counts.

I have in mind something like an estimator along the lines of truncated mean or Winsorized mean (that is discarding some data), or an M-estimator along the lines of Huber loss giving lower weight to outliers.

Remarks

• Thread Robust estimation of Poisson distribution discusses solving a similar problem in R (robust package description). I am however more interested in the underlying math (which I intend to implement myself in Python.)
• Other R packages: poissonMT, robeth, robustbase

Answer 1

Yes, there are. To name just one, which I've had good experiences with, you can minimize the Cramer-von Mises distance between the empirical distribution and the theoretical distribution with estimated parameters, possibly throwing out a specified percentage of the lowest and highest data points (or not; see below.) The C-vM distance is:

$$\omega^2 = \int_{-\infty}^{\infty} [F_n(x) - F(x;\theta)]^2\,\mathrm{d}F(x;\theta)$$

Obviously, in the case of a discrete distribution, you would use the sum, and you would integrate / sum over a prespecified subset of the data, making the appropriate adjustments in the calculation of $F(x;\theta)$.

Some general links are below:

The Consistency and Robustness of Modified Cramer-Von Mises and Kolmogorov-Cramer Estimators

Robust Weighted Cramer-von Mises estimators of location...

Robust estimation vis minimum distance methods

Robust Estimation with Exponentially Tilted Hellinger Distance

To see why this estimator is robust against small amounts of contamination (we will address the larger amounts of contamination case below), consider the function itself. The range of values of the square term is $[0,1]$; no matter how big an outlier in a sample of size $N$ we see, the influence of it on the squared term is limited by the $1/N$ change it can inflict on the empirical CDF. Likewise, as an outlier becomes more extreme, $\mathrm{d}F(x;\theta)$ goes to zero, further limiting its influence. There still is some influence, as the outlier reduces the empirical CDF of all the non-outlier data points by a factor of $(N-1)/N$.

To illustrate this, we construct an example. We generate 100 data points from a $\mathrm{Poisson(1)}$ distribution, and replace the $100^{th}$ with an increasing outlier. We calculate the C-vM estimate of the parameter $\theta$, and plot how the estimate changes as the outlier increases:

cvm.calc <- function(theta, x) {
  ub <- max(qpois(0.9999, theta)+1, length(x))
  if (ub > length(x)) {
    x <- c(x, rep(1, ub-length(x)))
  }
  
  sum_over_range <- 0:(ub-1)
  sum((x-ppois(sum_over_range, theta))^2 * dpois(sum_over_range, theta))
}

x <- rpois(100, 1)
estimate <- rep(0,51)
for (outlier in 0:50) {
  x[100] <- outlier
  xtbl <- table(factor(x, 0:max(x)))
  x_edf <- cumsum(xtbl) / sum(xtbl)
  estimate[outlier+1] <- optimize(cvm.calc, interval=c(0.1,5), x=x_edf)$minimum
}

plot(estimate~I(0:50), type="b", pch=16,
     main="Sensitivity Curve - 1% Contaminated Poisson(1)",
     xlab="Outlier (1 % contamination)", 
     ylab="C-vM Estimate")

mean(x[1:99])  # The sample mean of the first 99 observations
[1] 1.030303
mean(x)        # The sample mean with the 100th observation = 50
[1] 1.52

with the result:

A similar plot, this time with 10% of the data replaced by outliers, results in:

with a sample mean in the outlier <- 50 case of 5.9.

In the more general case, where we expect larger amounts of contamination, we can always use a truncated version of both the empirical distribution function and the theoretical distribution function, e.g., dropping the top $\alpha \%$ of the data and adjusting the theoretical distribution accordingly.

An obvious question is: how efficient is it in the zero-outlier case? A simulation will help to answer that in the $\mathrm{Poisson(1)}$ case:

cvm_estimate <- rep(0,10000)
for (i in seq_along(cvm_estimate)) {
  x <- rpois(100, 1)
  xtbl <- table(factor(x, 0:max(x)))
  x_edf <- cumsum(xtbl) / sum(xtbl)
  cvm_estimate[i] <- optimize(cvm.calc, interval=c(0.1,5), x=x_edf)$minimum
} 

mean((cvm_estimate - 1)^2)
[1] 0.01102528

which is only slightly worse than the MSE of the sample mean $(0.01)$.

Answer 2

If the issue merely boils down to very high or very low observations, one would be tempted to just use a trimmed mean. The problem with that of course, is that your estimate may be biased. You could say bias shmias! We often trade an unbiased estimator for a biased one that achieves much better variance when the overall MSE is superior. In my toy example:

set.seed(123)
out <- replicate(1e3, {
  y <- rpois(1e3, 20)
  y[1:50] <- rpois(50, 5)
  y[51:100] <- rpois(50, 50)
  q <- quantile(y, c(0.025, 0.975))
  ysub <- y[y > q[1] & y < q[2]]
  c(mean(y), mean(ysub))
})

The untransformed mean outperforms the trimmed mean though the trimmed mean does have a lower bias

> rowMeans(out)
[1] 20.74856 20.32148
> rowMeans((out-20)^2)
[1] 0.5816081 0.1286018

If an unbiased estimate is truly what you're after, a neat trick could be to trim the data, and then estimate the resulting Poisson density as a truncated one using maximum likelihood. Consider the example below:

set.seed(123)
y <- rpois(1e4, 20)
 
## only 1% of the data formally outlying
y[1:50] <- rpois(50, 5)
y[51:100] <- rpois(50, 50)

## trim 5% and treat it like a truncated Poisson
q <- quantile(y, c(0.025, 0.975))
ysub <- y[y > q[1] & y < q[2]]


negloglik <- function(lambda, x, q) {
  -sum(dpois(x=x, lambda=lambda, log=T) - 
    log(ppois(q=q[2], lambda=lambda) - ppois(q=q[1], lambda=lambda))
  )
}

nlm(negloglik, 10, x=ysub, q=q)

Answer 3

For a Poisson distribution $f_\theta(y)=\frac{e^{_\theta}\theta^y}{y!},y\in\mathbb{N}_0$, an M-estimator is given by $$ \sum_{i=1}^n\psi(y_i, \hat{\theta})=0, $$ where $$ \psi_{k,a}(y_i,\theta)=\frac{y_i-\theta-a\sqrt{\theta}}{\sqrt{\theta}} \min\left[1,\frac{k\sqrt{\theta}}{y_i-\theta-a\sqrt{\theta}}\right]. $$ This estimator is equivalent to the one based on Huber loss in case of estimation of location for a Normal distribution.

Parameter $a$ is necessary to correct for bias (the problem nicely discussed by @AdamO), and satisfies $$ E\psi_{k,a}(Y,\theta)=0 $$ In practice, the $a$ is calculated using the sample average and the current estimate of $\theta$, in the course of iteratively reweighted least squares procedure (IRWLS), commonly used for M-estimators.

This estimator is implemented in the R packages cited in the OP. The particular form of the estimator is due to Cadigan and Chen, Properties of robust m-estimators for Poisson and negative binomial data. A clear presentation, accessible to non-statisticians, can be found in thesis Robust modelling of count data by Elsaied, which also proposes a similar generalization of Tukey's M-estimator.