R/summary_ztpreg.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

#' @title Object summaries for Zero-truncated Poisson and Negative Binomial
#' Regression models and its dispaly
#'
#' @description \code{summary} is a generic function for summarizing results
#' produced from Zero-truncated Poisson and Negative Binomial models
#' with maximum likelihood estimation.
#' Using the estimates of fitted object \code{obj}, the size of
#' unknown population \code{N} and its \code{100*(1-alpha)} confidence
#' interval are calculated.
#'  
#' @param object a fitted model object,
#' @param alpha a confidence level required (\code{alpha=0.05} by
#' default), 
#' @param k an integer, the penalty per parameter to be used; the
#' \code{k=2} is used for classical AIC (Akaike's information
#' criteria), 
#' @param j an integer, for Lagrange multiplier test
#' @param ... other arguments
#' @param LM.test Lagrange multiplier test
#' @param HT.est Horvitz-Thompson estimator
#'
#' @return An object of class \code{summary.ztpr} with components
#' including
#' \item{tab}{Formatted object to be printed including coefficients,
#' standard errors, etc.  Significance star is also printed with stars,}
#' \item{aic}{Akaike Information Criteria,}
#' \item{N}{Estimate of unknown population size,}
#' \item{VarN}{Variance of \code{N},}
#' \item{ciu}{Upper bound of confidence interval,}
#' \item{cil}{Lower bound of confidence interval.}
#'
#' @seealso \code{\link{AIC}}, \code{\link{extractAIC}}, \code{\link{printCoefmat}}
#' @name summary.ztpr
#' @aliases print.summary.ztpr
#'
#' @note \code{print.summary.ztpr} is a generic function that displays
#' all summaries of object generated from \code{summary.ztpr}.
#' 
#' @author Chel Hee Lee <\email{gnustats@@gmail.com}>
#' @method summary ztpreg
#' @S3method summary ztpreg
summary.ztpreg <- function(object, alpha=0.05, k=2, j=1, LM.test=FALSE, HT.est=FALSE, ...){
  obj <- object

#  if(se.check) se <- sqrt(diag(vcv))
#  else stop("SE of coefficients are not all positive")

  X <- as.matrix(obj$X)
  obj$n <- n <- nrow(X)
  obj$p <- p <- ncol(X)
  y <- as.vector(obj$y)
  dist <- obj$dist 
  cfs <- obj$cfs
  lshape <- cfs[p+1]
  vcv <- as.matrix(obj$vcv)
  se.check <- obj$se.check
  se <- suppressWarnings(sqrt(diag(vcv)))
  # se <- sqrt(diag(vcv))
  llk <- obj$llk
  
  # stopifnot( (dist=="nbinom" & !LM.test) ) 
  if(dist=="nbinom" & LM.test) stop("'LM.test=TRUE' cannot be used with 'dist=nbinom'")

  z <- cfs/se
  pval <- 2*pnorm(-abs(z))
  tab <- cbind(cfs, se, z, pval)
  colnames(tab) <- c("Estimate", "SE", "z-score", "Pr(>|z|)")
  obj$tab <- tab
  
  obj$aic <- aic <- -2*llk + k*(p+(dist=="nbinom"))
  obj$df <- df <- p+(dist=="nbinom")
  
  cfs <- cfs[seq_len(p)]
  vcv <- vcv[seq_len(p), seq_len(p)]
  
  obj$mu <- mu <- as.vector(exp(X%*%cfs))
  scaler <- switch(dist, 
    "poisson"=ppois(q=0, lambda=mu, lower.tail=FALSE, log.p=FALSE), 
    "nbinom" =pnbinom(q=0, size=exp(lshape), mu=mu, lower.tail=FALSE, log.p=FALSE)
  )
  
  ## Test of overdispersion (chi-square with df=1)
  obj$LM.test <- LM.test
  if(LM.test){
    obj$yhat <- yhat <- mu/scaler
    obj$res <- res <- y - yhat

    af <- mu/(exp(mu)-1)
    XX <- lapply(as.list(as.data.frame(t(X))), function(x) outer(x,x))

    bb <- Reduce("+", mapply("*", yhat*(1-af), XX, SIMPLIFY=FALSE))
    ba <- as.vector(0.5*colSums(mu*yhat*af*X))
    ab <- 0.5*Reduce("+", mapply("*", mu^j*yhat*af, as.list(as.data.frame(t(X))), SIMPLIFY=FALSE))
    aa <- sum(0.5*mu^(2*j-1)*yhat*(1-0.5*mu*af))

    num <- sum(mu^(j-1)*(res^2-y+(res+y)*af))

    t <- (num/(2*(aa - ab%*%ginv(bb)%*%ba)^0.5))^2
    obj$LM.chisq <- c(t)
  }

  ## Estimation of a population size
  obj$HT.est <- HT.est
  if(HT.est){
    obj$N <- N <- sum(1/scaler)

    w <- exp(log(mu)-mu)/scaler^2
    cXw <- colSums(X*w)
    Var1 <- t(cXw) %*% vcv %*% cXw
    Var2 <- switch(dist, 
        "poisson"=sum(ppois(q=0, lambda=mu, lower.tail=TRUE, log.p=FALSE)/scaler^2), 
        "nbinom"=sum(pnbinom(q=0, size=exp(lshape), mu=mu, lower.tail=TRUE, log.p=FALSE)/scaler^2)
        )  
    obj$VarN <- VarN <- Var1 + Var2
    seN <- suppressWarnings(sqrt(VarN)) # due to NaN in Hessian ('nbinom')
    obj$ciu <- N + qnorm(1-alpha/2)*seN
    obj$cil <- N + qnorm(alpha/2)*seN
    obj$alpha <- alpha
  }
  
  class(obj) <- "summary.ztpreg"
  return(obj)
}
NULL