Distribución Binomial Negativa

Question

Estoy tratando de comparar los intervalos promedio de parto (número de días entre dos veces consecutivas que una vaca da a luz) entre diferentes razas de ganado. La variable llamada "intervalo de parto" es una variable discreta, y creo que tiene una distribución binomial negativa.

¿Cuál sería una forma sencilla de comprobar si una variable tiene realmente una distribución binomial negativa, usando R?

Answer 1

Usar una prueba de Chi-cuadrado como se explica en https://stats.stackexchange.com/a/17148/919 . El R El código que figura a continuación implementa dicha prueba, con los valores predeterminados apropiados para los datos de parto.

La prueba de Chi cuadrado es apropiada para tales conjuntos de datos discretos, como se explica en ese enlace. Para ver que podría funcionar para decidir si un conjunto de datos en particular es consistente con una distribución Binomial Negativa, aquí están los resultados de la simulación de mil conjuntos de datos de 180 valores independientes.

Los dos primeros conjuntos de datos se muestran en el gráfico de dispersión de la izquierda (los emparejamientos son arbitrarios). A continuación se muestra un histograma de los valores p de Chi cuadrado. Sus pequeñas desviaciones de una distribución uniforme (mostradas por la línea gris horizontal) son atribuibles al azar, lo que indica claramente que esta prueba proporciona los valores p correctos cuando la hipótesis nula (de una distribución binomial negativa) es cierta.

El poder de esta prueba es su capacidad para discriminar el Binomio Negativo de otras distribuciones. Para los datos típicos de parto, el Binomio Negativo es cercano al Normal (permitiendo el redondeo al día más cercano). También lo son otras distribuciones, como las distribuciones de Poisson con parámetros apropiados. Así pues, no deberíamos esperar mucho de esta prueba (o cualquier otra prueba de este tipo). Las distribuciones de los valores p de los datos simulados con las distribuciones de Poisson y Normal aparecen en los dos histogramas de la derecha. Debido a que hay una tendencia a que los valores p sean menores con estas alternativas, la prueba tiene cierto poder para detectar la diferencia. Pero debido a que los valores p no son muy pequeños muy a menudo, la potencia es baja: con un conjunto de datos de 180, será difícil distinguir el Binomio Negativo de los datos de Poisson de los datos de Normal. Esto sugiere que la cuestión de si los datos son consistentes con una distribución de Binomio Negativo podría tener poco significado o utilidad inherente.

---

Los parámetros de este ejemplo provienen de Werth, Azzam y Kinder, Intervalos de parto en vacas de carne a los 2, 3 y 4 años de edad cuando la cría no está restringida después del parto. J. Animal Sci. 1996, 74:593-596 . Debido a que este documento no proporciona estadísticas descriptivas adecuadas, he estimado la media y la varianza (y he establecido las rupturas para la prueba de ji cuadrado) a partir de esta cifra:

---

Este es el R para implementar los cálculos y gráficos que se muestran aquí. No es a prueba de balas: antes de aplicar cualquiera de estas funciones a otros conjuntos de datos, sería prudente probarlos y quizás incluir código para verificar que las estimaciones de máxima probabilidad son correctas.

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

Tienes razón, el wolframio está mal. Podría suceder...

Sólo usted debe corregir su exposición de la definición. Usted dice:

By definition, if blah blah, then bleh bleh

deberías decir:

By definition, blah blah, if bleh bleh

De hecho, demuestras que bleh bleh tiene bla bla.

Answer 3

Si está decidido a utilizar una distribución discreta, piense en los días como una variable de conteo y pruebe un Poisson en lugar de un binomio negativo. Entonces puedes usar un glm en R. (Sé que la gente se enfadará y dirá que uses un continuo, pero quizás corriendo ambos puedas ver por qué).

así que haz algo como esto

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

entonces, si estás realmente fijado en el binomio negativo, puedes cargar el paquete MASS y usarlo:

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

Algunos comentarios:

Algunas formas de ver si la forma que eligió después del modelo de Poisson es correcta:

dirigir summary(modp) y mira la desviación residual. Si es mayor que los grados de libertad de la desviación residual, entonces tienes un mal ajuste. Necesitarás hacer algunas cosas:

Primero, comprueba si hay valores atípicos usando halfnorm(residuals(modp) Si no hay ninguna sorpresa allí, entonces intenta mirar la variación. Puedes hacer algo como

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

Asegúrate de que la varianza es proporcional (se mueve con la media, ya que Poisson sólo tiene un parámetro). También querrás que los puntos se vean lanzados al azar alrededor del abline. Así que si están por encima o por debajo, entonces tienes un problema de dispersión.

Para solucionar el problema de la dispersión, puedes usar family = quasipoisson . así que haz un nuevo modp1 <- glm(Y ~ X1 + X2, family = quasipoisson, data) y luego mira tu resumen de nuevo: summary(modp1)

si tu desviación residual sigue siendo demasiado alta, entonces necesitas transformar tus variables o, más probablemente, tienes un error de especificación, es decir, una distribución errónea.

Para comprobar si hay un mal ajuste aquí, puedes usar una prueba de chi-cuadrado

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

Querrás que esto sea más grande que un nivel de importancia, así que algo así como 0.05 o 0.1. Si es más pequeño que esto, entonces todavía tienes un mal ajuste y puedes probar el código binomial negativo de arriba.