R/histogram.R

In HWxtest: Exact Tests for Hardy-Weinberg Proportions

# Plot histogram of HW test statistics
# (c) William R. Engels, 2014


#' Functions to plot a histogram of test statistic
#' 
#' Running the function \code{\link{hwx.test}} can create data for a frequency distribution plot of one of the four test statistics provided the parameter \code{histobins} is positive. These plotting data are contained within the \dQuote{hwtest} object generated by \code{\link{hwx.test}}. When that object is printed, a plot is drawn by \code{ggplot2}. If the user wishes to capture the \code{ggplot2} results, such as to use in making a composite figure, the \dQuote{gg} object can be obtained by calling \code{makeHistogram}.
#' 

#' @param x output from \code{\link{hwx.test}}
#' @param color1 The color for outcomes fitting the null distribution better than the observed
#' @param color2 The color for outcomes deviating from the null at leasst as much as observed. Area of \code{color2} is the P value.
#' @param curveColor color for the asymptotic distribution curve
#' @param ostats Observed statistics for the 4 test measures, \code{LLR}, \code{Prob}, \code{U} and \code{Chisq}.
#' @param statID Value 1-4 indicating which statistic to use for the plot.
#' @param m vector of allele counts
#' 
#' @return \code{defaultHistobounds} returns a vector containing the left and right boundaries for the x axis. This function is not normally called by the user
#' 
#' @return \code{makeHistobram} A graphic object of class \dQuote{gg} or \dQuote{ggplot} from \code{ggplot2}.
#' 
#' @seealso \code{\link{hwx.test}}


#' @rdname histogramFunctions
#' @export
defaultHistobounds <- 
function(ostats, statID, m) {
	k <- length(m);
	b <- double(2);
	n <- sum(m)/2;
	df <- k *(k-1)/2;
	if(statID==1 || statID==4) b[2] <- qchisq(.999, df);
	if(statID==2) {
		# find maximum probability
		ae <- matrix(0,k,k);
		for(i in 1:k) for(j in 1:k) ae[i,j] <- m[i]*m[j]/(2*n);
		for(i in 1:k)ae[i,i] <- ae[i,i]/2;
		b[1] <- (-2)*log(observedProb(ae));
		b[2] <- b[[1]] + qchisq(.999, df);
	}
	if(statID==3) {
		seu <- sqrt(n * (k-1));
		b[[1]] <- -4 * seu;
		b[[2]] <- 4 * seu;
	}
	return(b);	
}


#' @rdname histogramFunctions
#' @export
makeHistogram <- function(x, curveColor="blue", color1 = "gray40", color2 = "lightcoral") {
	requireNamespace("ggplot2")
	k <- length(x$alleles)
	labels <- c("-2 ln(LR)", "-2 ln(Probability)", "U Score (test for homozygote excess)", "Pearson's Chi Squared")
	statID <- x$statID
	showCurve <- x$showCurve
	os <- x$observed[statID]
	histobins <- x$histobins
	histobounds <- x$histobounds
	
	if (statID == 1) 
		os <- -2 * os
	if (statID == 2) 
		os <- -2 * log(os)
	nleft <- as.integer(histobins * (os - histobounds[1])/(histobounds[2] - histobounds[1]))
	if (nleft < 0) 
		nleft <- 0
	nright <- histobins - nleft + 1
	if (nright < 0) 
		nright <- 0
	if (statID == 3 && os < 0) { #swap colors
		temp <- color1
		color1 <- color2
		color2 <- temp
		labels[3] <- "U Score (test for heterozygote excess)"
	}
	statvals <- seq(histobounds[1], histobounds[2], length.out = histobins + 1)

	if (showCurve && (statID == 1 || statID == 4)) {
		# if (ntrials == 0) 
			# ntrials <- 1
		chi <- dchisq(statvals, k * (k - 1)/2)
	} else chi = 0
	colrs <- c(rep(I(color1), nleft), rep(color2, nright))
	if(length(colrs) > histobins+1) length(colrs) <- histobins+1
	lwd <- 200/histobins

	d <- data.frame(statvals = statvals, histoData = x$histoData, colrs = colrs, chi = chi)
	g <- ggplot2::ggplot(d, ggplot2::aes_string(x = 'statvals', y = 'histoData', color = 'colrs')) + ggplot2::geom_segment(ggplot2::aes(xend = statvals, yend = 0), size = lwd, 
		show.legend = F) + ggplot2::scale_color_identity() + ggplot2::ylab("Density") + ggplot2::xlab(labels[statID])

	if (showCurve && (statID == 1 || statID == 4)) {
		g <- g + ggplot2::geom_line(ggplot2::aes(y = chi), color = curveColor, size = 1)
	}
	if(!is.null(x$sampleName)) g <- g + ggplot2::ggtitle(x$sampleName)
	return(g)
}