R/summary_imprecise.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

#' @aliases summary summary.imprecise
#'
#' @title 
#' Summarizing Imprecise Class Objcects
#'
#' @description
#' The function \code{summary.imprecise} is the \code{imprecise} method of 
#' generic function \code{summary} which is used for producing summaries of
#' the results from various model fitting functions performed on the imprecise
#' inferential framework.
#' See \sQuote{Details}.
#'
#' @param object 
#' the object of class \code{imprecise}.
#'
# @param prob 
# a numeric vector of probabilities with values in [0,1].
#' 
#' @param HT.est 
#' a logical value; Defaults to \code{FALSE}.
#' Would like to see Horvitz-Thompson estimator?
#' 
#' @param silent
#' a logical value; Would like to see the report of progress messages 
#' (including warnings and error messages) generated from numerical methods
#' that are used for a numerical approximation?
#' Defaults to \code{FALSE}.
#'
#' @details
#' Object summaries depend on the option chosen in the \code{method}.
#' The commonly returned object is \code{est} which is a collection of 
#' all Bayes estimates from the imprecise posterior.
#' 
#' The Horvitz-Thompson estimator is the implementation of that in the study
#' of \cite{Heijden 2003}.
#' \code{HT.est=TRUE} is not allowed when \code{ztrunc=TRUE}.
#'
#' @return 
#' A list with the following compoentns:
#' \item{est}{a list of vectors or matrices of the estimates of canonical 
#' parameters.}
#' \item{N}{Horvitz-Thompson estimates if \code{HT.est=TRUE}}
#' \item{mui}{Estimate of the mean parameter if \code{xi} is supplied.}
#' \item{inf}{a numeric vector of length \code{p}, having an infimum of each 
#' canonical parameter.}
#' \item{sup}{a numeric vector of length \code{p}, having a supremum of each
#' canonical parameter.}
#' \item{inf.idx}{an index vector corresponding to \code{inf} in \code{est}.}
#' \item{sup.idx}{an index vector corresponding to \code{sup} in \code{est}.}
#' \item{delta}{the vector of imprecision}
#' 
#' @references
#' Peter GM van der Heijden, Rami Bustami, Maarten JLF
#' Cruyff, Godfried Engbersen and Hans C van Houwelingen (2003), Point
#' and interval estimation of the population size using the truncated
#' Poisson regression model, \emph{Statistical Modelling}, \bold{3},
#' pp. 305-322. 
#'
#' Lee (2013) ``Imprecise inferential framework'', PhD thesis.
#' 
#' Walley P (1991). _Statistical reasoning with imprecise probabilities_,
#' series Monographs on statistics and applied probability. Chapman and
#' Hall. ISBN 9780412286605, <URL:
#' http://books.google.ca/books?id=Nk9Qons1kHsC>.
#'
#' @keywords 
#' imprecision, Horvitz-Thompson estimator
#' 
#' @concept
#' imprecision, prior ignorance, prior modelling, imprecise prior,
#' quantification of prior ignorance, indeterminate decision
#' 
#' @author Chel Hee Lee <\email{gnustats@@gmail.com}>
#' @method summary imprecise
#' @S3method summary imprecise
summary.imprecise <- function(object, HT.est=FALSE, silent=FALSE, ...){
  
  mc <- match.call()
                              
  # check and update the stage on the inferential framework
  if (object$stage != "update") {
    stop("Not correct sequence of imprecise inferential framework.\n",
         sQuote(mc[[1]]), " is followed by ", 
         sQoute("update()"), ".\n")
  }
  object$stage <- "summary"

  # sanity check
  stopifnot(!missing(object), is.logical(HT.est))
#  stopifnot(is.numeric(prob), is.vector(prob))
#  stopifnot(is.numeric(alpha), is.vector(alpha))
  stopifnot(is.logical(HT.est), length(HT.est)==1)
  stopifnot(is.logical(silent), length(silent)==1)
  
  # naming conventions
  xtms <- object$xtms
  xreg <- object$xreg
  ztrunc <- object$ztrunc
  method <- object$method
  apriori <- object$apriori
  y <- object$y
  X <- object$X
  m1 <- object$m1
  xi <- object$xi
  m0shape <- object$m0shape
  formula <- object$formula
  object$HT.est <- HT.est
    
  if (!ztrunc & HT.est) {
			stop("HT.est=", sQuote(HT.est), "cannot be used with ztrunc=", 
         sQuote(ztrunc))
  }

  
  if (!silent) {
    message("Input sanity check .................... PASS! \n",
            "in ", mc[[1]], ".\n")
  }
    
  # definitions for estimating 'mui' and 'N'
  fnxregmui <- function(x, ...){
    # Computing the value of a linear predictor 'mu_i' 
    #
    # Args:
    #   x: The estimate of 'beta'
    #
    # Return:
    #   The evaluated scalar 'mu_i'
    
    beta <- as.vector(x)
    val <- exp(xi%*%beta)
    return(val)
  }

  fnxregN <- function(x, ...){
    # Computing the Horvitz-Thompson estimate of 'N' in the context of 
    # a generalized linear model
    # 
    # Args: 
    #   x: The estimate of 'beta'
    #
    # Return:
    #   The evaluated scalar 'N'
    
    beta <- as.vector(x)
    mu <- exp(X%*%beta)
    val <- sum(1/ppois(q=0, lambda=mu, lower.tail=FALSE, log.p=FALSE))
    return(val)
  }

  fnN <- function(t, ...) {
    # Computing the Horvtiz-Thompson estimate of 'N' in the context of 
    # models that do not involve with covariates
    #
    # Args:
    #   t: The estimate of a canonical parameter 'theta'
    # 
    # Return:
    #   The evaluated scalar 'N'
    
    val <- length(y)/ppois(q=0, lambda=exp(t), lower.tail=FALSE, log.p=FALSE)
    return(val)
  }
    
  # tracking the correct extreme point 
  xid <- do.call(c, lapply(m1, "[[", "xid"))
  # nxid <- length(xid)
  
  # General summarization for all numerical methods
  if (xreg) { # 'nrow:xtms' x 'p:cfs'
    est <- as.data.frame(do.call(rbind, lapply(m1, "[[", "cfs")))
    inf <- do.call(c, lapply(est, min))
    sup <- do.call(c, lapply(est, max))
    inf.idx <- do.call(c, lapply(est, which.min))
    sup.idx <- do.call(c, lapply(est, which.max))
    tslpars <- est
  } else {   # 'length:xtms' x '1:theta'
    est <- do.call(c, lapply(m1, "[[", "theta"))
    names(est) <- xid
    inf <- min(est)
    sup <- max(est)
    inf.idx <- which.min(est)
    sup.idx <- which.max(est)
    tslpars <- do.call(rbind, lapply(m1, "[[", "eta"))
  }
  
  delta <- abs(sup-inf)
  object <- c(object, list(est=est, inf=inf, sup=sup, inf.idx=inf.idx, 
                           sup.idx=sup.idx, delta=delta, tslpars=tslpars))
  
  if (!is.null(object$qcLA)) {
    message("Minimal summaries are kept for imprecise inferential framework.")
    invisible(object)
  }
  
  # additional summaries for 'MH'
  if (method == "MH" & is.null(object$qcLA)) {
  
    if (xreg) {
      MHchain <- lapply(lapply(m1, "[[", "MHchain"), as.data.frame)
    } else {
      MHchain <- as.data.frame(do.call(cbind, lapply(m1, "[[", "MHchain")))
    }
    accept.rate <- do.call(c, lapply(m1, "[[", "accept.rate"))
    object$MHchain <- MHchain
    object$accept.rate <- accept.rate
    
    if (ztrunc) {
      Nchain <- muichain <- list()
      
      for (i in xid) {
        if (!silent) {
          message(xid[i], " is processing now ...")
        }
        mci <- as.matrix(MHchain[[i]])
        
        if (xreg) {
          Nchain[[i]] <- apply(mci, 1, fnxregN)
          muichain[[i]] <- apply(mci, 1, fnxregmui)
        } else {
          Nchain[[i]] <- fnN(mci)
          muichain[[i]] <- exp(mci)
        }
      }
      
      N <- do.call(c, lapply(Nchain, mean))
      mui <- do.call(c, lapply(muichain, mean))
      
      object$N <- N
      object$mui <- mui
      object$muichain <- muichain
      object$Nchain <- Nchain
    }
  }

  # additional summaries for 'IS' 
  if(method == "IS" & is.null(object$qcLA)){
  
    if (xreg) {
      isample <- lapply(m1, "[[", "isample")
      lweight <- lapply(m1, "[[", "lweight")
    } else {
      isample <- as.data.frame(do.call(cbind, lapply(m1, "[[", "isample")))
      lweight <- as.data.frame(do.call(cbind, lapply(m1, "[[", "lweight")))
    }
    ess <- do.call(c, lapply(m1, "[[", "ess"))
    
    object$isample <- isample
    object$lweight <- lweight
    object$ess <- ess
    
    if (ztrunc) {
      Nsmp <- muismp <- list()
      
      for(i in xid) {
        if (!silent) {
          message(xid[i], "is processing now ...")
        }
        isi <- as.matrix(isample[[i]])
        lwi <- as.vector(lweight[[i]])
        
        if (xreg) {
          Nsmp[[i]] <- apply(isi, 1, fnxregN)
          muismp[[i]] <- apply(isi, 1, fnxregmui)
        } else {
          Nsmp[[i]] <- fnN(isi)
          muismp[[i]] <- isi*exp(lwi)/sum(exp(lwi))
        }
      }
      
      N <- do.call(c, lapply(Nsmp, function(x){
        val <- sum(x*exp(lwi)/sum(exp(lwi)))
        return(val)
        }))
      mui <- do.call(c, lapply(muismp, function(x){
        val <- sum(x*exp(lwi)/sum(exp(lwi)))
        return(val)
        }))
      
      object$N <- N
      object$mui <- mui
      object$Nsmp <- Nsmp
      object$muismp <- muismp
    }
  }

  # additional summaries for 'LA'
  if (method=="LA" & is.null(object$qcLA)) {
    
    sratio <- do.call(rbind, lapply(m1, "[[", "sratio"))
    fdiff <- do.call(rbind, lapply(m1, "[[", "fdiff"))
    niter <- do.call(rbind, lapply(m1, "[[", "niter"))
    colnames(niter) <- colnames(fdiff) <- colnames(sratio) <- colnames(X)
    object <- c(object, list(sratio=sratio, fdiff=fdiff, niter=niter))
    
    if (!is.null(xi)) {
      if (xreg) {
        object$mui <- do.call(c, lapply(m1, "[[", "mui"))
      } else {
        object$mui <- NA # TODO
      }
    }
    
    if (ztrunc & HT.est) {
      if (xreg) {
        object$N <- do.call(c, lapply(m1, "[[", "N"))
      } else {
        object$N <- fn(est)
      }
    }  
  }

  class(object) <- c("summary.imprecise")
  return(object)
}
NULL