6
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ by eye. Fit a distribution to your data and look if it fits good. You can use K-S test, but in practice 1) goodness-of-fit of your distribution will depend on your fitted parameters, 2) usually no real data follows exactly the theoretical distribution. $\endgroup$ – German Demidov Oct 13 '17 at 15:49
  • $\begingroup$ stats.stackexchange.com/questions/20648/… $\endgroup$ – German Demidov Oct 13 '17 at 15:50
  • 4
    $\begingroup$ @German You link to a thread about continuous variables, with the K-S test as the prominent recommendation. This test does not apply to discrete variables because the sampling distribution of the test statistic differs from the one for continuous variables. People usually perform goodness of fit tests using GLMs or (more simply, if there are sufficient data) chi-squared tests. $\endgroup$ – whuber Oct 13 '17 at 16:17
  • $\begingroup$ @whuber thank you! Will google for it, never thought about that. $\endgroup$ – German Demidov Oct 13 '17 at 16:41
  • $\begingroup$ Having grown up on a farm I would be very surprised if the distribution is negative binomial! My initial guess at a simple distribution would be const_1 + const_2*Geometric $\endgroup$ – P.Windridge Oct 13 '17 at 17:58
4
$\begingroup$

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.

Figure

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:

Figure 1 from the paper


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))
$\endgroup$
1
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ If the births are recorded only to the day, then the data you ask for simply don't exist. What would you propose doing then? $\endgroup$ – whuber Oct 13 '17 at 20:36
  • $\begingroup$ @whuber, I meant that days count is an approximation. So, using the the continuous distribution is actually going to be closer to the underlying process. This is not a good case where you say that continuous distribution is an approximation of discrete. It's actually the reverse case: the discrete distribution would be an approximation $\endgroup$ – Aksakal Oct 13 '17 at 20:38
  • 2
    $\begingroup$ Regardless, insisting that we use a continuous distribution to describe discrete data seems ill-supported. (It can be justified as an approximation to the discrete distribution when no individual value has an appreciable probability.) And if you were to do that, you ought to choose something more appropriate: the negative binomial is an analog of a Gamma distribution. Exponential distributions are not even plausible for modeling elapsed times between births: a pregnancy has to intervene, while Exponential distributions assign the highest probability density to an immediate second birth! $\endgroup$ – whuber Oct 13 '17 at 20:41
  • $\begingroup$ @whuber, FP numbers in a computer are discrete. So, when modeling a continuous Gaussian variable with FP numbers in a computer nobody would argue that the Binomial distribution must be used. It's exactly the same argument here you're making. It doesn't look like Neg Binomial is used in the field for calving intervals. I see people using the standard survival analysis tools: hazards, kaplan meyer etc. $\endgroup$ – Aksakal Oct 13 '17 at 21:00
  • 1
    $\begingroup$ Let' not confuse how we do our calculations with how we conceive of and model our problems. We shouldn't care (much) either whether a particular distribution is conventionally used in a field--if we did, we would never make any progress. Why not address the question that was asked? Let's explain how to test whether these data are consistent with some negative binomial distribution. $\endgroup$ – whuber Oct 13 '17 at 22:42
0
$\begingroup$

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.

$\endgroup$

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.