R/gldrmPIT.R

In gldrm: Generalized Linear Density Ratio Models

#' Confidence intervals for gldrm coefficients
#'
#' Plots and returns the randomized probability inverse transform of a 
#' fitted gldrm.
#'
#' @param gldrmFit A gldrm model fit. Must be an S3 object of class "gldrm",
#' returned from the \code{gldrm} function. The matrix of semiparametric
#' tilted probabilities must be returned, which is done by fitting gldrm with
#' \code{gldrmControl = gldrm.control(returnfTiltMatrix = TRUE)}.
#' @param nbreaks Number of breaks in the histogram.
#' @param cex.main Text size for main titles.
#' @param cex.lab Text size for axis labels.
#' @param cex.axis Text size for axis numbers.
#' 
#' @return Randomized robability inverse transform as a vector. Also plots the 
#' histogram and uniform QQ plot.
#' 
#' @details
#' The probability inverse transform is defined generally as \eqn{\hat{F}(y|x)},
#' which is the fitted conditional cdf of each observation evaluated at the
#' observed response value. In the case of gldrm, the fitted cdf is descrete, so
#' we draw a random value from a uniform distribution on the interval 
#' (\eqn{\hat{F}(y|x)}, \eqn{\hat{F}(y-|x)}), where \eqn{y-} is the next largest 
#' observed support less than \eqn{y} (or -Infinity if \eqn{y} is the minimum 
#' support value). The output and plots generated by this function will vary
#' slightly each time it is called (unless the random number generator seed is 
#' set beforehand).
#' 
#' @examples
#' data(iris, package="datasets")
#'
#' ### Fit gldrm and return fTiltMatrix
#' fit <- gldrm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Species,
#'              data=iris, link="log")
#'              
#' # Probability inverse transform plot
#' gldrmPIT(fit)
#' 
#' @export
gldrmPIT <- function(gldrmFit, nbreaks=7, cex.main=NULL, cex.lab=NULL, cex.axis=NULL) {
    
  if (class(gldrmFit) != "gldrm")
      stop(paste0("gldrmFit must be an S3 object of class 'gldrm'",
                  " (returned from the gldrm function)"))
  if (!("fTiltMatrix" %in% names(gldrmFit)))
      stop("gldrmFit must be run with gldrmControl = gldrm.control(returnfTiltMatrix = TRUE)")
    
  # Extract response variable - it is always the first column of $data
  y <- gldrmFit$data[, 1]
  n <- length(y)
  
  # observed support and estimated fTiltMatrix
  spt <- gldrmFit$spt
  fTiltMatrix <- gldrmFit$fTiltMatrix
  
  # compute upper and lower cumulative probabilities
  Fy <- vector(length=n)
  Fy.minus.one <- vector(length = n)
  for(i in 1:n) {
    Fy[i] <- sum(fTiltMatrix[i,spt <= y[i]])
    Fy.minus.one[i] <- sum(fTiltMatrix[i,spt < y[i]])
  }
  
  # randomised prob inverse transform (Smith, 1985)
  PIT <- runif(n, min=Fy.minus.one, max=Fy)
  par(mfrow=c(1,2))
  hist(PIT, breaks=nbreaks, xlab="Prob. Inv. Transform", 
       cex.main=cex.main, cex.lab=cex.lab, cex.axis=cex.axis)
  plot(sort(PIT),1:n/n, ylab = "Uniform quantiles", xlab = "Prob. Inv. Transform", 
       main = "PIT-Uniform QQ plot", 
       cex.main=cex.main, cex.lab=cex.lab, cex.axis=cex.axis)
  abline(0,1)
  return(PIT)
}