Correções e adições para exibição dos resultados

 - Corrige teste F para familia quasipoisson, as diferencas de
   deviances são escolanadas pelo respectivo parametro de
   dispersão do modelo mais completo (dois a dois)
 - Adciona função apra teste de hipoteses sobre o parametro de
   dispersão da familia quasipoisson.
 - Adiciona função para extração e comporação dos coeficientes
   estimados sob os modelos Poisson, CMP, BinNeg e Quasi-Poisson

docs/_setup.R

@@ -127,8 +127,9 @@ myanova <- function(...) {
             "P(>Chisq)" = c(NA, pvs))
     }
     if (test == "F") {
-        sigma <- summary(model.list[[which.max(nps)]])$dispersion
-        df.sigma <- model.list[[which.max(nps)]]$df.residual
+        ## ver pág. 187 expresssão 7.10 em Modern Applied, Ripley
+        sigma <- sapply(model.list, function(x) summary(x)$dispersion)
+        df.sigma <- sapply(model.list, function(x) x$df.residual)
         dev <- sapply(model.list, function(x) deviance(x))
         dfs <- c(NA, -diff(dev))
         dnp <- c(NA, diff(nps))
@@ -145,6 +146,24 @@ myanova <- function(...) {
     return(tab)
 }
 
+##======================================================================
+## Testa o parametro de dispersão de um modelo quasi poisson
+quasitest <- function(...) {
+    quasi.list <- list(...)
+    cls <- sapply(quasi.list, function(model) model$family$family)
+    if (sum(!cls %in% "quasipoisson") > 0) {
+        stop("O modelo não é quasipoisson")
+    }
+    ##-------------------------------------------
+    dev <- sapply(quasi.list, function(model) deviance(model))
+    dfs <- sapply(quasi.list, function(model) model$df.residual)
+    sig <- sapply(quasi.list, function(model) summary(model)$dispersion)
+    pvs <- pchisq(q = dev, df = dfs)*2
+    out <- cbind("sigma2" = sig, `P(>Chisq)` = pvs)
+    rownames(out) <- rep("sigma == 1", length(quasi.list))
+    return(out)
+}
+
 ##======================================================================
 ## Definições de plot.profile.mle2 para perfil de log-verossimilhança
 
@@ -211,3 +230,55 @@ mycorrplot <- function(Corr, ...) {
         diag = "l", tl.pos = "lt", tl.col = "black",
         tl.cex = 0.8, col = brewer.pal(9, "Greys")[-(1:3)], ...)
 }
+
+##======================================================================
+## Tabela de comparação dos coeficientes dos modelos para dados de
+## contagem
+
+##-------------------------------------------
+## Extrai os coeficientes dos modelos
+getCoefs <- function(model, rownames = NULL) {
+    cls <- class(model)[1]
+    if (!cls %in% c("glm", "negbin", "mle2")) {
+        stop("Classe de modelos não contemplada")
+    }
+    est <- switch(
+        cls,
+        "glm" = {
+            fam <- model$family$family
+            if (!fam %in% c("poisson", "quasipoisson")) {
+                stop("Família de distribuições não contemplada")
+            }
+            switch(
+                fam,
+                "poisson" = {
+                    rbind(c(NA, NA), summary(model)$coef[, c(1, 3)])
+                },
+                "quasipoisson" = {
+                    rbind(c(summary(model)$dispersion, NA),
+                          summary(model)$coef[, c(1, 3)])
+                }
+            )
+        },
+        "negbin" = {
+            rbind(c(model$theta, model$SE.theta),
+                  summary(model)$coef[, c(1, 3)])
+        },
+        "mle2" = {
+            summary(model)@coef[, c(1, 3)]
+        }
+    )
+    colnames <- c("estimate", "est/std.error")
+    if (!is.null(rownames)) {
+        rownames(est) <- rownames
+    }
+    return(est)
+}
+##-------------------------------------------
+## Faz a tabela de comparação
+coeftab <- function(..., rownames = NULL) {
+    model.list <- list(...)
+    coefs <- lapply(model.list, getCoefs, rownames = rownames)
+    tab <- do.call(cbind, coefs)
+    return(tab)
+}