Robust estimators for count data

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

• Which parameters or distribution characteristics do you want to estimate?
– Dave
Commented Jan 4, 2023 at 14:01
• It depends on the distribution. For now I assume Poisson and calculate the mean of the counts. But in the literature they use negative binomial for similar problems (as over-dispersed Poisson.) Commented Jan 4, 2023 at 14:10
• But what do you want to estimate? Mean? Variance? Cubed kurtosis divided by IQR?
– Dave
Commented Jan 4, 2023 at 14:12
• Then you're not estimating the parameter of a Poisson distribution. Please clarify preferably in an edit of the original question, what you want to estimate. I think you want to estimate the mean, but only you can say for sure what you want to do.
– Dave
Commented Jan 4, 2023 at 14:50
• You say you want to estimate some parameters. What parameters are those? You will (or at least can) use different approaches to estimate the mean than you would to estimate the variance.
– Dave
Commented Jan 4, 2023 at 14:58

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)$$.

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)$$.

• Thanks for the tip. I will look into this direction. Commented Jan 4, 2023 at 15:52
• Do you have any intuition why this estimator would be robust against mixture and outlier relative to the UMVUE (sample mean)? Commented Jan 4, 2023 at 18:03
• @AdamO - yes, and I'll edit in a longer response to that excellent question in a bit. Commented Jan 4, 2023 at 18:34
• @AdamO - Version 1 of the update is in place! Commented Jan 4, 2023 at 19:35

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)

• Thanks. People working on the same problem do use trimmed and winsorized means, as if the counts were normally distributed. However, it seems to me problematic 1) when the data points are few, and 2) because anomalously small counts and anomalously big may have different effect and come for different reasons. So you hit the right nail regarding the bias. Commented Jan 4, 2023 at 18:16
• @RogerVadim you'd have to explore the scenario with simulation to know for any given data analysis, it would depend on the extent of mixing (1%? 10%?), and the extremity of the outlying observations. I showed another example of the 5/50 contamination having 10% infiltation. Commented Jan 4, 2023 at 18:31
• Although easy to compute, trimmed estimators are typically non differentiable and the calculation of the variance is difficult. Commented Jan 4, 2023 at 20:08
• @utobi the Winsorized estimate of variance is a jackknife, which is fairly easy to implement. My method is actually even more complicated. I'm not sure of the asymptotic variance, although I'd recommend a bootstrap which would become computationally expensive. Commented Jan 4, 2023 at 20:32

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.