5
$\begingroup$

I have some count data I need to fit to a distribution. I have tried several distributions but none of them start out steep enough.

Here is the fit of a geometric distribution:

enter image description here

My question is: Which distribution should I fit here? (And (how) is that possible in R?)

In case you want to play:

Due to its size I put the data on pastebin. You can read it into R using the pastebin package:

## library("devtools")
## devtools::install_github("hrbrmstr/pastebin")
library("pastebin")
ttt <- get_paste("hkuUg0m3")$text
ttt <- as.numeric(strsplit(ttt, "\r\n")[[1]])

I fit the geometric distribution using the fitdistrplus package:

library("fitdistrplus")
fit <- fitdist(ttt, distr = "geom")

And this is what I use for visualization:

viz_fitdist(fit, distr = "geom")

using this function:

viz_fitdist <- function(fit, distr="geom", annotation_colour="#326a97")
{
  library("ggplot2")
  library("grid")
  library("gridExtra")

  data <- data.frame(dat=fit$data)

  plotdata <- as.data.frame(table(data))
  plotdata <- as.data.frame(table(data)/nrow(data))
  plotdata$data <- as.numeric(as.character(plotdata$data))
  p <- plotdata %>% ggplot(aes(x=data, y=Freq)) + geom_bar(stat="identity")

  if (distr=="geom") {
    xsupport <- 1:max(data$dat)
    } else {
      xsupport <- unique(data$dat)
    }
  if (distr!="none") {
    dfun <- paste0("d", distr)
    y <- do.call(dfun, c(list(x = xsupport), as.list(fit$estimate)))
  }

  if (distr!="none") {
    fitdata <- data.frame(x=xsupport, y=y)

    p <- p +
      geom_point(aes(x=x, y=y), data=fitdata, colour=annotation_colour)

    p <- p + geom_line(aes(x=x, y=y), data=fitdata, colour=annotation_colour, size=1)

    p <- p +
      annotation_custom(grob=textGrob(sprintf("Parameter: %.2f", fit$estimate),
                                      x=1-0.07,
                                      y=1-0.07,
                                      just=c("right", "top")))
  }

  p + scale_x_continuous("Data", breaks=xsupport) + ylab("Frequency")
}

EDIT: more distributions

Here is the fit of a negative binomial:

ttt %>%
  fitdist(distr = "nbinom") %>%
  viz_fitdist(distr = "nbinom")

enter image description here

And a zero inflated poisson:

library("extraDistr")
ttt %>%
  fitdist(distr = "zip", start = list(lambda = 1.6, pi = 0.5), lower = c(-Inf, 0), upper = c(Inf, 1)) %>%
  viz_fitdist(distr = "zip")

enter image description here

And a Weibull:

ttt %>%
  fitdist(distr = "weibull") %>%
  viz_fitdist(distr = "weibull")

enter image description here

$\endgroup$
  • 2
    $\begingroup$ Welcome to Cross Validated! Please take a moment to view our tour. There are several possible answers to this. How closely does the model need to match your data? What are you using as your measure of "fit?" $\endgroup$ – Tavrock Mar 4 '17 at 21:37
  • $\begingroup$ One purpose of the fit is to have a means to quantify 'unexpectedness' for high counts. For that I guess the model does not need to be very close. My current measure of 'fit' are my eyes... $\endgroup$ – Andreas Mar 4 '17 at 21:46
  • 2
    $\begingroup$ The point of fitting a distribution to your data is obscure. How exactly will you use it? What do these data represent? $\endgroup$ – whuber Mar 4 '17 at 22:23
  • $\begingroup$ I want to identify unexpectedly high counts in the spirit of alphaOutlier for which I need a probability model. $\endgroup$ – Andreas Mar 5 '17 at 11:48
5
$\begingroup$

Let me open by saying that searching for and identifying a model that appears to describe the data is not where I would normally begin my modelling; it begins with the information whuber and I were seeking on how the model will be used, with what is being measured - and further, to consider any associated theoretical or expert practical understanding of its properties, the performance of models on similar data sets, and so on. Considerations like these are essential to the choice of an appropriate model.

Nevertheless, as requested, I will attempt to describe some steps one could take in examining data like this and seeing if it could suggest a model.

If you're going to use the resulting distributional model for inference (hypothesis tests, or confidence intervals, or prediction intervals for example), there's a problem that you chose the model based on your data so the nominal properties of the inference will no longer hold. (If you really can't use the characteristics of the variable and subject area knowledge, having no other basis than the sample itself, then it's better to pull off a random subset and use that for model identification.)

Given the heavy tail, one option is to look at the way the log of the empirical survivor function relates to $x$. My first plots were to plot $\log(\hat{S}(x))$ against $\log(x)$ (not shown; linear at the top but then somewhat curved near the middle) and then against $\sqrt x$. The appearance of these then suggested to me that a fair first model might be a zeta distribution. However, they can also suggest a number of other possibilities (such as perhaps directly modelling the seeming linear relationship in the first plot below; it might provide a somewhat better fit, at the expense of being less familiar)

So I then looked at the empirical pmf, plotting $\log(\hat{p}(x))$ against $\log{x}$ as a check on the potential for fitting a zeta function.

enter image description here

That indeed suggests a fairly straight line relationship, and it looks like a parameter in the ballpark of 2 or perhaps a bit more would be a reasonable first model -- note that the most precisely estimated proportions are the first few, and fitting those accurately will matter more for a lot of applications -- but not all. However, it's not a perfect fit - with a large sample size like you have, I believe you would easily reject a pure zeta model (or indeed any other simple model).

enter image description here

(NB the parameter value here has not been optimized; it's just a reasonable guess at a suitable value)

So if I had to suggest a simple model on the basis of only looking at the data, the zeta distribution would be an obvious candidate; if an initial guess at the parameter value is needed, I'd start at 2.

On fitting power laws, such as the zeta distribution, also see Cosma Shalizi's essay So you think you have a power law? which debunks some of the enthusiasm for power laws (a lot of things will look more or less power-law like) and the paper it discusses, which has a lot of useful information.

$\endgroup$
  • $\begingroup$ Thanks for this answer! I appreciate your comment about fitting meaningful distributions. But even more I appreciate you answering my question of which simple distribution might fit the given data. And most of all, I appreciate your very educational explanations on how you got there. So, +3 would be appropriate. Thanks! $\endgroup$ – Andreas Mar 12 '17 at 20:18
0
$\begingroup$

Have you tried a Negative Binomial distribution or a zero-inflated Poisson? These distributions are often used to take care of overdispersion.

You have been suggested below a Weibull distribution, which may fit the data but it is not appropriate for count data. A Weibull is a distribution for continuous measurements, and not for discrete data. Negative Binomial distributions are routinely used for RNAseq and Microbiome data, among other types of data. So are the zero-inflated Poisson distribution. They are standard tools and are easily implemented in R packages See Cameron-Trivedi (2013, Regression Analysis of Count Data). Now, there are also Discrete Weibull distributions, but they are not so much used. See the R package "DiscreteWeibull" and related manuscripts.

The negative binomial fit above is similar to the Weibull and it is more appropriate, since it takes into account the nature of the data.

If you really wanna be flexible, mixtures of Poisson-Gamma may help you to capture overdispersion in a Bayesian fit. They are of course not trivial. But if the question is about fitting the data, then all the options should be proposed. Now, your data show that histogram because there are probably bumps at the right end of the x-axis. No parametric model can capture those bumps correctly. That's why I was suggesting mixtures of parametric models. And in such case, the Bayesian fit is probably the way to go.

The right model depends of course on the objective of your inference, which is not clear from your question.

$\endgroup$
  • $\begingroup$ A well chosen Weibul could be made to fit this as well. $\endgroup$ – Tavrock Mar 4 '17 at 21:48
  • 3
    $\begingroup$ A Weibull distribution (which is continuous) makes little sense for counted data! $\endgroup$ – whuber Mar 4 '17 at 22:21
  • 2
    $\begingroup$ Many of those categories have probably zero frequency. Except a few. Again, these type of situations/needs for model fit are standard in the analysis of sequencing data. But we haven't been told what these data represent. $\endgroup$ – M G Mar 5 '17 at 1:24
  • 3
    $\begingroup$ @Tavrock Weibull is continuous. Data is (mostly small) positive integers (>90% of values are 1,2 or 3). Makes no sense to use Weibull for that. You could perhaps make an argument for modeling the top 2 or 3 percent of the distribution as continuous but that seems unnecessary when so much of it is very small integers. $\endgroup$ – Glen_b Mar 5 '17 at 2:09
  • 2
    $\begingroup$ @Andreas It's quite easy to find a discrete distribution that fits this data pretty well (a couple of quick plots was sufficient), but that's going about the problem backwards. You need to talk about what you're modelling (what did you count?*) and what the model will be used for. $\quad$ * while I could probably make a fair guess at the general kind of thing you're counting, we shouldn't guess $\endgroup$ – Glen_b Mar 5 '17 at 2:31

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.