# R/gldrmPIT.R In gldrm: Generalized Linear Density Ratio Models

#### Documented in gldrmPIT

#' 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,
#'
#' # 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)
}


## Try the gldrm package in your browser

Any scripts or data that you put into this service are public.

gldrm documentation built on April 13, 2018, 9:04 a.m.