8
$\begingroup$

I have a set of numbers which are assumed to be coming from a Poisson distribution. The set has some outliers also and because of that, maximum likelihood estimates are badly affected. I heard that robust estimation procedures can help in such a situation. Can any one explain how to do this? I am not a statistics student.

I found that the glmrob function in R can be used for this. (I am quite new to R). But I could not figure out, how to use that despite reading the manual pages. In particular I am unable to understand how to get a forumula which is the first argument to the glmrob function.

Thanks.

$\endgroup$
  • 1
    $\begingroup$ Have you read the glmrob documentation? There is sample code for a robust Poisson fit at p. 23. $\endgroup$ – whuber Apr 9 '12 at 22:01
  • $\begingroup$ I did look at them. In that example in page 23, there is a formula like this sumY ~ Age10 + Base4*Trt. With the data I have, I am unable to come up with any such formulae. I just have a bunch of numbers assumed to be coming from a Poisson distribution. As I understand, I do not have a dependent variable and independent variable and may be that is why I am unable to come up with a formula. $\endgroup$ – suresh Apr 9 '12 at 22:45
  • 2
    $\begingroup$ You implicitly have an intercept term. $\endgroup$ – cardinal Apr 10 '12 at 0:05
  • 2
    $\begingroup$ You have poisson counts that are your responses and your "predictors" consist of just an intercept term, i.e., you want to model it as if all of your responses come from the same distribution. Does that make more sense? :) $\endgroup$ – cardinal Apr 10 '12 at 0:16
  • 1
    $\begingroup$ Example formula: y~1. $\endgroup$ – cardinal Apr 10 '12 at 0:26
17
$\begingroup$

@cardinal has telegraphed an answer in comments. Let's flesh it out. His point is that although general linear models (such as implemented by lm and, in this case, glmRob) appear intended to evaluate relationships among variables, they can be powerful tools for studying a single variable, too. The trick relies on the fact that regressing data against a constant is just another way of estimating its average value ("location").

As an example, generate some Poisson-distributed data:

set.seed(17)
x <- rpois(10, lambda=2)

In this case, R will produce the vector $(1,5,2,3,2,2,1,1,3,1)$ of values for x from a Poisson distribution of mean $2$. Estimate its location with glmRob:

library(robust)
glmrob(x ~ 1, family=poisson())

The response tells us the intercept is estimated at $0.7268$. Of course, anyone using a statistical method needs to know how it works: when you use generalized linear models with the Poisson family, the standard "link" function is the logarithm. This means the intercept is the logarithm of the estimated location. So we compute

exp(0.7268)

The result, $2.0685$, is comfortably close to $2$: the procedure seems to work. To see what it is doing, plot the data:

plot(x, ylim=c(0, max(x)))
abline(exp(0.7268), 0, col="red")

Plot with fitted line

The fitted line is purely horizontal and therefore estimates the middle of the vertical values: our data. That's all that's going on.

To check robustness, let's create a bad outlier by tacking a few zeros onto the first value of x:

x[1] <- 100

This time, for greater flexibility in post-processing, we will save the output of glmRob:

m <- glmrob(x ~ 1, family=poisson())

To obtain the estimated average we can request

exp(m$coefficients)

The value this time equals $2.496$: a little off, but not too far off, given that the average value of x (obtained as mean(x)) is $12$. That is the sense in which this procedure is "robust." More information can be obtained via

summary(m)

Its output shows us, among other things, that the weight associated with the outlying value of $100$ in x[1] is just $0.02179$, almost $0$, pinpointing the suspected outlier.

$\endgroup$
  • 3
    $\begingroup$ (+1) Nice presentation, as usual. :) $\endgroup$ – cardinal Apr 10 '12 at 14:13
  • $\begingroup$ Thanks a lot. Such an answer is very important for me now because I am completely new to all this (terms like predictor, intercept, and even the package R). $\endgroup$ – suresh Apr 10 '12 at 20:20
  • $\begingroup$ Can you point to a similar tool for a bivariate Poisson distribution? $\endgroup$ – Diogo Santos Aug 13 at 7:55
  • $\begingroup$ @Diogo It depends on what that tool is intended to do and precisely what form of dependence among the marginals is being used. $\endgroup$ – whuber Aug 13 at 13:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.