1# Copyright 2004 Virgilio Gómez Rubio and Roger Bivand
2
3test.nb.pois <- function(x.nb, x.glm){
4	if (!(inherits(x.nb, "negbin"))) stop("not a negbin object")
5	if (!(inherits(x.glm, "glm"))) stop("not a glm object")
6	zscore <- -log(x.nb$theta)*x.nb$theta/x.nb$SE.theta
7	u <- pnorm(zscore)
8	pz <- 2*min(u, 1-u)
9	probs <- x.nb$theta/(x.nb$theta+x.nb$fitted.values)
10	lrt <- 2*(sum(dnbinom(x.nb$y, x.nb$theta, probs, log=TRUE))-
11		sum(dpois(x.nb$y, x.glm$fitted.values, log=TRUE)))
12	names(lrt) <- "LR"
13	pchi <- pchisq(lrt, df=1, lower.tail=FALSE)
14	vec <- c(zscore, pz)
15	names(vec) <- c("zscore", "p.mayor.modZ")
16	res <- list(estimate=vec, statistic=lrt, p.value=pchi, parameter=1,
17		method="Likelihood ratio test for overdispersion",
18		data.name=paste(deparse(substitute(x.nb)), ":",
19		deparse(substitute(x.glm))))
20	class(res) <- "htest"
21	res
22}
23
24DeanB <- function(x.glm, alternative="greater") {
25	alternative <- match.arg(alternative, c("less", "greater", "two.sided"))
26	if (!(inherits(x.glm, "glm"))) stop("not a glm object")
27	y <- model.response(model.frame(x.glm))
28	mu <- fitted(x.glm)
29	Pb <- sum((y-mu)^2-y)/sqrt(2*sum(mu^2))
30	names(Pb) <- "P_B"
31	Pv <- NA
32	if (is.finite(Pb)) {
33			if (alternative == "two.sided")
34				Pv <- 2 * pnorm(abs(Pb), lower.tail=FALSE)
35			else if (alternative == "greater")
36				Pv <- pnorm(Pb, lower.tail=FALSE)
37			else Pv <- pnorm(Pb)
38	}
39	res <- list(statistic=Pb, p.value=Pv, alternative=alternative,
40		method="Dean's P_B test for overdispersion",
41		data.name=deparse(substitute(x.glm)))
42	class(res) <- "htest"
43	res
44}
45
46DeanB2 <- function(x.glm, alternative="greater"){
47
48	y <- model.response(model.frame(x.glm))
49	mu <- fitted(x.glm)
50	h <- hatvalues(x.glm)
51	Pb2 <- sum((y-mu)^2-y+h*mu)/sqrt(2*sum(mu^2))
52	names(Pb2) <- "P'_B"
53	Pv <- NA
54	if (is.finite(Pb2)) {
55			if (alternative == "two.sided")
56				Pv <- 2 * pnorm(abs(Pb2), lower.tail=FALSE)
57			else if (alternative == "greater")
58				Pv <- pnorm(Pb2, lower.tail=FALSE)
59			else Pv <- pnorm(Pb2)
60	}
61	res <- list(statistic=Pb2, p.value=Pv, alternative=alternative,
62		method="Dean's P'_B test for overdispersion",
63		data.name=deparse(substitute(x.glm)))
64	class(res) <- "htest"
65	res
66}
67