我正在尝试比较不同品种的牛之间的平均产犊间隔(母牛连续两次分娩之间的天数)。称为“产犊间隔”的变量是一个离散变量,我相信它具有负二项分布。
什么是测试变量是否真的有否定的简单方法。二项分布,使用 R?
我正在尝试比较不同品种的牛之间的平均产犊间隔(母牛连续两次分娩之间的天数)。称为“产犊间隔”的变量是一个离散变量,我相信它具有负二项分布。
什么是测试变量是否真的有否定的简单方法。二项分布,使用 R?
使用卡方检验,如https://stats.stackexchange.com/a/17148/919中所述。下面的R
代码实现了这样的测试,默认值适用于产犊数据。
如该链接所述,卡方检验适用于此类离散数据集。要查看确定特定数据集是否与负二项分布一致可能有效,以下是模拟 180 个独立值的一千个数据集的结果。
前两个数据集显示在左侧的散点图中(配对是任意的)。接下来显示卡方 p 值的直方图。它与均匀分布(由水平灰线显示)的小偏差可归因于偶然性,强烈表明当零假设(负二项分布)为真时,该检验提供了正确的 p 值。
该测试的强大之处在于它能够将负二项式与其他分布区分开来。对于典型的产犊数据,负二项式接近正常(允许四舍五入到最近的日期)。其他分布也是如此,例如具有适当参数的泊松分布。因此,我们不应该对这个测试抱有太多期望(或任何此类测试)。来自具有泊松分布和正态分布的模拟数据的 p 值分布出现在右侧的两个直方图中。由于这些备选方案的 p 值趋于较小,因此该检验具有一定的能力来检测差异。但是因为 p 值通常不是很小,所以功效很低:对于 180 的数据集,很难区分负二项式和泊松与正态数据。这表明数据是否与负二项分布一致的问题可能没有多少内在意义或有用性。
此示例的参数来自 Werth、Azzam 和 Kinder,在产犊后不限制繁殖的情况下,肉牛在 2、3 和 4 岁时的产犊间隔。 J. 动物科学。1996, 74:593-596。因为这篇论文没有提供足够的描述性统计数据,我从这个图中估计了均值和方差(并设置了卡方检验的间隔):
这是R
实现此处显示的计算和绘图的代码。这不是万无一失的:在将这些函数中的任何一个应用于其他数据集之前,谨慎的做法是对其进行测试,并可能包含验证最大似然估计值是否正确的代码。
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))
我不认为负二项式是该变量分布的合理首选。是的,天数是一个离散的数字,但事件之间的真实间隔是连续的:奶牛不会在一天中的特定时间准确分娩。它只是发生,以便您以天为单位测量间隔。因此,没有理由从离散分布开始。底层分布当然不是离散的。如果您要测量一头母牛在 5 年内的出生次数,那将是固有的离散量,并且会要求离散概率分布。
在你的情况下,我的第一个猜测是尝试指数分布。
如果您打算使用离散分布,则将天数视为计数变量并尝试使用泊松而不是负二项式。然后你可以在 R 中使用 glm 。(我知道人们会生气并说使用连续,但也许同时运行你会明白为什么)。
所以做这样的事情
modp<- glm(Y ~ X1 + X2, family = poisson, data)
那么如果你真的设置了负二项式,你可以加载 MASS 包并使用:
modnb <- glm.nb(Y ~ X1 + X2, data)
一些评论:
看看你在泊松模型后选择的形式是否正确的一些方法:
运行summary(modp)
并查看剩余偏差。如果它大于剩余偏差自由度,那么您的拟合度会很差。你需要做几件事:
halfnorm(residuals(modp)
首先,使用如果没有任何意外,检查异常值,然后尝试查看方差。你可以做类似的事情
plot(log(fitted(modp)), log((data$Y-fitted(modp))^2), xlab=
expression(hat(mu)),ylab=expression((y-hat(mu))^2))
abline(0,1)
确保方差成比例(随均值移动,因为泊松只有一个参数)。您还希望这些点看起来随机分布在 abline 周围。因此,如果它们全部超过或全部低于,那么您就会遇到分散问题。
要解决色散问题,您可以使用family = quasipoisson
. 所以做一个新的modp1 <- glm(Y ~ X1 + X2, family = quasipoisson, data)
,然后再看看你的总结:summary(modp1)
如果您的剩余偏差仍然太高,那么您需要转换您的变量,或者更有可能的是,您有一个规范错误,即错误的分布。
要检查此处是否不合适,您可以使用卡方检验
pchisq(deviance(modp1),df.residual(modp1),lower=FALSE)
您会希望它大于一个显着性水平,例如 0.05 或 0.1。如果它小于这个值,那么你仍然不合适,你可以尝试上面的负二项式代码。