R/qb_pi.R

In predint: Prediction Intervals

#' Simple uncalibrated prediction intervals for quasi-binomial data
#'
#' \code{qb_pi()} is a helper function that is internally called by  \code{quasi_bin_pi()}. It
#' calculates simple uncalibrated prediction intervals for binary
#' data with constant overdispersion (quasi-binomial assumption).
#'
#' @param newsize number of experimental units in the historical clusters.
#' @param histsize number of experimental units in the future clusters.
#' @param pi binomial proportion
#' @param phi dispersion parameter
#' @param q quantile used for interval calculation
#' @param alternative either "both", "upper" or "lower"
#' \code{alternative} specifies, if a prediction interval or
#' an upper or a lower prediction limit should be computed
#' @param histdat additional argument to specify the historical data set
#' @param newdat additional argument to specify the current data set
#' @param algorithm used to define the algorithm for calibration if called via
#' \code{quasi_bin_pi}. This argument is not of interest for the calculation
#' of simple uncalibrated intervals
#'
#' @details This function returns a simple uncalibrated prediction interval
#' \deqn{[l,u]_m = n^*_m \hat{\pi} \pm q \sqrt{\hat{\phi} n^*_m \hat{\pi} (1- \hat{\pi}) +
#' \frac{\hat{\phi} n^{*2}_m \hat{\pi} (1- \hat{\pi})}{\sum_h n_h}}}
#'
#' with \eqn{n^*_m} as the number of experimental units in the \eqn{m=1, 2, ... , M} future clusters,
#' \eqn{\hat{\pi}} as the estimate for the binomial proportion obtained from the
#' historical data, \eqn{\hat{\phi}} as the estimate for the dispersion parameter
#' and \eqn{n_h} as the number of experimental units per historical cluster. \cr
#'
#' The direct application of this uncalibrated prediction interval to real life data
#' is not recommended. Please use the \code{beta_bin_pi()} functions for real life applications. \cr
#'
#' @return \code{qb_pi} returns an object of class \code{c("predint", "quasiBinomailPI")}.
#' @export
#'
#' @importFrom stats qnorm
#'
#' @examples
#' qb_pred <- qb_pi(newsize=50, pi=0.3, phi=3, histsize=c(50, 50, 30), q=qnorm(1-0.05/2))
#' summary(qb_pred)
#'
qb_pi <- function(newsize,
                  histsize,
                  pi,
                  phi,
                  q=qnorm(1-0.05/2),
                  alternative="both",
                  newdat=NULL,
                  histdat=NULL,
                  algorithm=NULL){

        # histsize must be numeric or integer
        if(!(is.numeric(histsize) | is.integer(histsize))){
                stop("histsize must be numeric or integer")
        }

        # newsize must be numeric or integer
        if(!(is.numeric(newsize) | is.integer(newsize))){
                stop("newsize must be numeric or integer")
        }

        # pi must be be smaller than 1 and bigger than
        if(pi<0){
                stop("pi must be bigger than 0 and smaller than 1")
        }

        if(pi>1){
                stop("pi must be bigger than 0 and smaller than 1")
        }

        # Phi must be numeric or integer
        if(!(is.numeric(phi) | is.integer(phi))){
                stop("phi must be numeric or integer")
        }

        # q must be numeric or integer
        if(!(is.numeric(q) | is.integer(q))){
                stop("q must be numeric or integer")
        }

        if(length(q) > 2){
                stop("length(q) > 2")
        }

        # alternative must be defined
        if(isTRUE(alternative!="both" && alternative!="lower" && alternative!="upper")){
                stop("alternative must be either both, lower or upper")
        }

        # check algorithm
        if(!is.null(algorithm)){
                if(algorithm != "MS22"){
                        if(algorithm != "MS22mod"){
                                stop("algoritm must be either NULL, MS22 of MS22mod")
                        }
                }
        }

        #-----------------------------------------------------------------------

        # historical number of experimental units
        histn <- sum(histsize)

        # Expected future observations
        y_star_hat <- newsize * pi

        # var for future random variable
        var_y_star <- phi * newsize * pi *(1-pi)

        # var for the expected future observations
        var_y_star_hat <- phi * newsize^2 * pi *(1-pi) / histn

        # SE for prediction
        pred_se <- sqrt(var_y_star + var_y_star_hat)

        #-----------------------------------------------------------------------
        ### Calculate the interval

        if(alternative=="both"){

                if(length(q) ==1){
                        lower <- y_star_hat - q * pred_se
                        upper <- y_star_hat + q * pred_se

                        out <- data.frame(lower,
                                          upper)
                }

                if(length(q) ==2){
                        lower <- y_star_hat - q[1] * pred_se
                        upper <- y_star_hat + q[2] * pred_se

                        out <- data.frame(lower,
                                          upper)
                }

        }

        if(alternative=="lower"){

                lower <- y_star_hat - q * pred_se

                out <- data.frame(lower)
        }

        if(alternative=="upper"){

                upper <- y_star_hat + q * pred_se

                out <- data.frame(upper)
        }

        # Output has to be an S3 object

        out_list <- list("prediction"=out,
                         "newsize"=newsize,
                         "newdat"=newdat,
                         "histsize"=histsize,
                         "histdat"=histdat,
                         "y_star_hat"=y_star_hat,
                         "pred_se"=pred_se,
                         "alternative"=alternative,
                         "q"=q,
                         "pi"=pi,
                         "phi"=phi,
                         "algorithm"=algorithm)

        out_s3 <- structure(out_list,
                            class=c("predint", "quasiBinomialPI"))

        return(out_s3)
}