Negative Binomial distribution

Question

I am trying to compare average calving intervals (number of days between two consecutive times that a cow gives birth) between different breeds of cattle. The variable called "calving interval" is a discrete variable, and I believe it has a negative binomial distribution.

What would be a simple way to test if a variable really has a neg. binomial distribution, using R?

Answer 1

Use a Chi-squared test as explained at https://stats.stackexchange.com/a/17148/919. The R code below implements such a test, with defaults appropriate for calving data.

The Chi-squared test is appropriate for such discrete datasets, as explained at that link. To see that it might work to decide whether a particular dataset is consistent with a Negative Binomial distribution, here are the results of simulating a thousand datasets of 180 independent values.

The first two datasets are shown in the scatterplot at left (the pairings are arbitrary). A histogram of the Chi-squared p-values is shown next. Its small deviations from a uniform distribution (shown by the horizontal gray line) are attributable to chance, strongly indicating this test provides the correct p-values when the null hypothesis (of a Negative Binomial distribution) is true.

The power of this test is its ability to discriminate Negative Binomial from other distributions. For typical calving data, the Negative Binomial is close to Normal (allowing for rounding to the nearest day). So are other distributions, such as Poisson distributions with appropriate parameters. Thus, we shouldn't expect much of this test (or any such test). The distributions of p-values from simulated data with Poisson and Normal distributions appear in the right two histograms. Because there is a tendency for p-values to be smaller with these alternatives, the test has some power to detect the difference. But because the p-values aren't very small very often, the power is low: with a dataset of 180, it will be difficult to distinguish Negative Binomial from Poisson from Normal data. This suggests that the question whether the data are consistent with a Negative Binomial distribution might have little inherent meaning or usefulness.

---

The parameters for this example come from Werth, Azzam and Kinder, Calving intervals in beef cows at 2, 3, and 4 years of age when breeding is not restricted after calving. J. Animal Sci. 1996, 74:593-596. Because this paper does not provide adequate descriptive statistics, I estimated the mean and variance (and set the breaks for the chi-squared test) from this figure:

---

This is the R code to implement the calculations and plots shown here. It's not bulletproof: before applying any of these functions to other datasets, it would be prudent to test them and perhaps include code to verify the maximum likelihood estimates are correct.

library(MASS) # rnegbin
#
# Specify parameters to generate data.
#
mu <- 360 # Mean days in interval
v <- 30^2 # Variance of days: must exceed mu^2
n <- 18000  # Sample size
n.sim <- 3e2 # Simulation size
#
# Functions to fit a negative binomial to data.
#
pnegbin <- function(k, mu, theta) {
  v <- mu + mu^2/theta         # Variance
  p <- 1 - mu / v              # "Success probability"
  r <- mu * (1-p) / p          # "Number of failures until the experiment is stopped"
  pbeta(p, k+1, r, lower.tail=FALSE)
}
# #
# # Test `pnegbin` by comparing it to randomly generated data.
# #
# z <- rnegbin(1e3, mu, theta)
# plot(ecdf(z))
# curve(pnegbin(x, mu, theta), add=TRUE, col="Red", lwd=2)
#
# Maximum likelihood fitting of data based on counts in predefined bins.
# Returns the fit and chi-squared statistics.
#
negbin.fit <- function(x, breaks) {
  if (missing(breaks))
      breaks <- c(-1, seq(-40, 30, by=10) + 365, Inf)
  observed <- table(cut(x, breaks))
  n <- length(x)

  counts.expected <- function(n, mu, theta) 
    n * diff(pnegbin(breaks, mu, theta))

  log.lik.m <- function(parameters) {
    mu <- parameters[1]
    theta <- parameters[2]
    -sum(observed * log(diff(pnegbin(breaks, mu, theta))))
  }

  v <- var(x)
  m <- mean(x)
  if (v > m) theta <- m^2 / (v - m) else theta <- 1e6 * m^2
  parameters <- c(m, theta)
  fit <- optim(parameters, log.lik.m)

  expected <- counts.expected(n, fit$par[1], fit$par[2])
  chi.square <- sum(res <- (observed - expected)^2 / expected)
  df <- length(observed) - length(parameters) - 1
  p.value <- pchisq(chi.square, df, lower.tail=FALSE)
  return(list(fit=fit, chi.square=chi.square, df=df, p.value=p.value, 
              data=x, breaks=breaks, observed=observed, expected=expected,
              residuals=res))
}
#
# Test on randomly generated data.
#
# set.seed(17)
sim <- replicate(n.sim, negbin.fit(rnegbin(n, mu, theta))$p.value)
#
# Generate data for illustration.
#
theta <- mu^2 / (v - mu)
x <- rnegbin(n, mu, theta)
y <- rnegbin(n, mu-4.3, theta)
#
# Display data and simulation.
#
par(mfrow=c(1,4))
plot(x-365, y-365, pch=15, col="#00000040",
     xlab="First calving interval", ylab="Second calving interval",
     main="Simulated Data")
abline(h=0)
abline(v=0)

hist(sim, freq=FALSE, xlab="p-values", ylab="Frequency", 
     main="Histogram of Simulated P-values",
     sub="Negative Binomial Data")
abline(h=1, col="Gray", lty=3)
#
# Simulate non-Negative Binomial data for comparison.
#
sim.2 <- replicate(n.sim, negbin.fit(rpois(n, mu))$p.value)
hist(sim.2, freq=FALSE, xlab="p-values", ylab="Frequency", 
     main="Histogram of Simulated P-values",
     sub="Poisson Data")
abline(h=1, col="Gray", lty=3)

sim.3 <- replicate(n.sim, negbin.fit(floor(rnorm(n, mu, sqrt(mu))))$p.value)
hist(sim.3, freq=FALSE, xlab="p-values", ylab="Frequency", 
     main="Histogram of Simulated P-values",
     sub="Normal Data")
abline(h=1, col="Gray", lty=3)
par(mfrow=c(1,1))

Answer 2

I don't think that negative binomial is a reasonable first choice of the distribution for this variable. Yes, the number of days is a discrete number, but the true interval between the events is continuous: cows do not give birth exactly at given hour of the day. It just happens so that you measure the interval in days. Therefore, there is not reason to start with discrete distributions. The underlying distribution is certainly not discrete. If you were measuring number of births a cow given in 5 years, that would be inherently discrete quantity, and would ask for a discrete probability distribution.

In your case, my first guess would be to try something like an exponential distribution.

Answer 3

If you are set on using a discrete distribution then think of the days as a count variable and try a Poisson instead of a negative binomial. Then you can just use a glm in R. (I know people will get angry and say use a continuous, but maybe running both you can see why).

so do something like this

modp<- glm(Y ~ X1 + X2, family = poisson, data)

then if you are really set on the negative binomial you can load the MASS package and use:

modnb <- glm.nb(Y ~ X1 + X2, data)

Some comments:

Some ways to see if the form you chose after the poisson model is correct:

run summary(modp) and look at the residual deviance. If it is greater than the residual deviance degrees of freedom then you have a bad fit. You will need to do a few things:

First, check for outliers using halfnorm(residuals(modp) If there aren't any surprises there, then try looking at the variance. You can do something like

plot(log(fitted(modp)), log((data$Y-fitted(modp))^2), xlab=
expression(hat(mu)),ylab=expression((y-hat(mu))^2))
abline(0,1)

Make sure the variance is proportional (moves with the mean, since Poisson only has one parameter). You will also want the points to look randomly thrown in around the abline. So if they are all over or all under, then you have a dispersion problem.

To fix the dispersion problem, you can use family = quasipoisson. so make a new modp1 <- glm(Y ~ X1 + X2, family = quasipoisson, data) and then look at your summary again: summary(modp1)

if your residual deviance is still too high, then you need to transform your variables or, more likely, you have a specification error, ie, wrong distribution.

To check for bad fit here, you can use a chi-square test

pchisq(deviance(modp1),df.residual(modp1),lower=FALSE)

You will want this to be greater than a level of significance, so something like 0.05 or 0.1. If it is smaller than this, then you still have a bad fit and you can try the negative binomial code above.