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