contrib/benchmarking_code_block/Rdhist.R

   1 #-------------------------------------------------------------------------------
   2 # Copyright (c) 2012 University of Illinois, NCSA.
   3 # All rights reserved. This program and the accompanying materials
   4 # are made available under the terms of the
   5 # University of Illinois/NCSA Open Source License
   6 # which accompanies this distribution, and is available at
   7 # http://opensource.ncsa.illinois.edu/license.html
   8 #-------------------------------------------------------------------------------
   9 ##--------------------------------------------------------------------------------------------------#
  10 ##' Variable-width (dagonally cut) histogram
  11 ##'
  12 ##'
  13 ##' When constructing a histogram, it is common to make all bars the same width.
  14 ##' One could also choose to make them all have the same area.
  15 ##' These two options have complementary strengths and weaknesses; the equal-width histogram oversmooths in regions of high density, and is poor at identifying sharp peaks; the equal-area histogram oversmooths in regions of low density, and so does not identify outliers.
  16 ##' We describe a compromise approach which avoids both of these defects. We regard the histogram as an exploratory device, rather than as an estimate of a density.
  17 ##' @name dhist
  18 ##' @title Diagonally Cut Histogram
  19 ##' @param x is a numeric vector (the data)
  20 ##' @param a is the scaling factor, default is 5 * IQR
  21 ##' @param nbins is the number of bins, default is assigned by the Stuges method
  22 ##' @param rx is the range used for the left of the left-most bin to the right of the right-most bin
  23 ##' @param eps used to set artificial bound on min width / max height of bins as described in Denby and Mallows (2009) on page 24.
  24 ##' @param xlab is label for the x axis
  25 ##' @param plot = TRUE produces the plot, FALSE returns the heights, breaks and counts
  26 ##' @param lab.spikes = TRUE labels the \% of data in the spikes
  27 ##' @return list with two elements, heights of length n and breaks of length n+1 indicating the heights and break points of the histogram bars.
  28 ##' @author Lorraine Denby, Colin Mallows
  29 ##' @references Lorraine Denby, Colin Mallows. Journal of Computational and Graphical Statistics. March 1, 2009, 18(1): 21-31. doi:10.1198/jcgs.2009.0002.
  30 dhist <- function(x, a=5*iqr(x),
  31                 nbins=nclass.Sturges(x), rx = range(x,na.rm = TRUE),
  32                 eps=.15, xlab = "x", plot = TRUE,lab.spikes = TRUE)
  33 {
  34
  35   if(is.character(nbins))
  36     nbins <- switch(casefold(nbins),
  37                     sturges = nclass.Sturges(x),
  38                     fd = nclass.FD(x),
  39                     scott = nclass.scott(x),
  40                     stop("Nclass method not recognized"))
  41   else if(is.function(nbins))
  42     nbins <- nbins(x)
  43
  44   x <- sort(x[!is.na(x)])
  45   if(a == 0)
  46     a <- diff(range(x))/100000000
  47   if(a != 0 & a != Inf) {
  48     n <- length(x)
  49     h <- (rx[2] + a - rx[1])/nbins
  50     ybr <- rx[1] + h * (0:nbins)
  51     yupper <- x + (a * (1:n))/n
  52                                         # upper and lower corners in the ecdf
  53     ylower <- yupper - a/n
  54                                         #
  55     cmtx <- cbind(cut(yupper, breaks = ybr), cut(yupper, breaks =
  56                                 ybr, left.include = TRUE), cut(ylower, breaks = ybr),
  57                   cut(ylower, breaks = ybr, left.include = TRUE))
  58     cmtx[1, 3] <- cmtx[1, 4] <- 1
  59                                         # to replace NAs when default r is used
  60     cmtx[n, 1] <- cmtx[n, 2] <- nbins
  61                                         #
  62                                         #checksum <- apply(cmtx, 1, sum) %% 4
  63     checksum <- (cmtx[, 1] + cmtx[, 2] + cmtx[, 3] + cmtx[, 4]) %%
  64     4
  65                                         # will be 2 for obs. that straddle two bins
  66     straddlers <- (1:n)[checksum == 2]
  67                                         # to allow for zero counts
  68     if(length(straddlers) > 0) {
  69       counts <- table(c(1:nbins, cmtx[ - straddlers, 1]))
  70     } else {
  71       counts <- table(c(1:nbins, cmtx[, 1]))
  72     }
  73     counts <- counts - 1
  74                                         #
  75     if(length(straddlers) > 0) {
  76       for(i in straddlers) {
  77         binno <- cmtx[i, 1]
  78         theta <- ((yupper[i] - ybr[binno]) * n)/a
  79         counts[binno - 1] <- counts[binno - 1] + (
  80                                                   1 - theta)
  81         counts[binno] <- counts[binno] + theta
  82       }
  83     }
  84     xbr <- ybr
  85     xbr[-1] <- ybr[-1] - (a * cumsum(counts))/n
  86     spike<-eps*diff(rx)/nbins
  87     flag.vec<-c(diff(xbr)<spike,F)
  88     if ( sum(abs(diff(xbr))<=spike) >1) {
  89       xbr.new<-xbr
  90       counts.new<-counts
  91       diff.xbr<-abs(diff(xbr))
  92       amt.spike<-diff.xbr[length(diff.xbr)]
  93       for (i in rev(2:length(diff.xbr))) {
  94         if (diff.xbr[i-1] <= spike&diff.xbr[i] <= spike &
  95             !is.na(diff.xbr[i])) {
  96           amt.spike <- amt.spike+diff.xbr[i-1]
  97           counts.new[i-1] <- counts.new[i-1]+counts.new[i]
  98           xbr.new[i] <- NA
  99           counts.new[i] <- NA
 100           flag.vec[i-1] <- T
 101         }
 102         else amt.spike<-diff.xbr[i-1]
 103       }
 104       flag.vec<-flag.vec[!is.na(xbr.new)]
 105       flag.vec<-flag.vec[-length(flag.vec)]
 106       counts<-counts.new[!is.na(counts.new)]
 107       xbr<-xbr.new[!is.na(xbr.new)]
 108
 109     }
 110     else flag.vec<-flag.vec[-length(flag.vec)]
 111     widths <- abs(diff(xbr))
 112     ## N.B. argument "widths" in barplot must be xbr
 113     heights <- counts/widths
 114   }
 115   bin.size <- length(x)/nbins
 116   cut.pt <- unique(c(min(x) - abs(min(x))/1000,
 117                      approx(seq(length(x)), x, (1:(nbins - 1)) * bin.size, rule = 2)$y, max(x)))
 118   aa <- hist(x, breaks = cut.pt, plot = FALSE, probability = TRUE)
 119   if(a == Inf) {
 120     heights <- aa$counts
 121     xbr <- aa$breaks
 122   }
 123   amt.height<-3
 124   q75<-quantile(heights,.75)
 125   if (sum(flag.vec)!=0) {
 126     amt<-max(heights[!flag.vec])
 127     ylim.height<-amt*amt.height
 128     ind.h<-flag.vec&heights> ylim.height
 129     flag.vec[heights<ylim.height*(amt.height-1)/amt.height]<-F
 130     heights[ind.h] <- ylim.height
 131   }
 132   amt.txt<-0
 133   end.y<-(-10000)
 134   if(plot) {
 135     barplot(heights, abs(diff(xbr)), space = 0, density = -1, xlab =
 136             xlab, plot = TRUE, xaxt = "n",yaxt='n')
 137     at <- pretty(xbr)
 138     axis(1, at = at - xbr[1], labels = as.character(at))
 139     if (lab.spikes) {
 140       if (sum(flag.vec)>=1) {
 141         usr<-par('usr')
 142         for ( i in seq(length(xbr)-1)) {
 143           if (!flag.vec[i]) {
 144             amt.txt<-0
 145             if (xbr[i]-xbr[1]<end.y) amt.txt<-1
 146           }
 147           else {
 148             amt.txt<-amt.txt+1
 149             end.y<-xbr[i]-xbr[1]+3*par('cxy')[1]
 150           }
 151           if (flag.vec[i]) {
 152             txt<-paste(' ',format(round(counts[i]/
 153                                         sum(counts)*100)),'%',sep='')
 154             par(xpd = TRUE)
 155             text(xbr[i+1]-xbr[1],ylim.height-par('cxy')[2]*(amt.txt-1),txt, adj=0)
 156           }}
 157       }
 158       else print('no spikes or more than one spike')
 159     }
 160     invisible(list(heights = heights, xbr = xbr))
 161   }
 162   else {
 163     return(list(heights = heights, xbr = xbr,counts=counts))
 164   }
 165 }
 166 #==================================================================================================#
 167
 168
 169 #--------------------------------------------------------------------------------------------------#
 170 ##' Calculate interquartile range
 171 ##'
 172 ##' Calculates the 25th and 75th quantiles given a vector x; used in function \link{dhist}.
 173 ##' @name iqr
 174 ##' @title Interquartile range
 175 ##' @param x vector
 176 ##' @return numeric vector of length 2, with the 25th and 75th quantiles of input vector x.
 177 iqr <- function(x){
 178   return(diff(quantile(x, c(0.25, 0.75), na.rm = TRUE)))
 179 }
 180 ##==================================================================================================#
 181
 182 ##--------------------------------------------------------------------------------------------------#
 183 ##' Creates empty ggplot object
 184 ##'
 185 ##' An empty base plot to which layers created by other functions
 186 ##' (\code{\link{plot.data}}, \code{\link{plot.prior.density}},
 187 ##' \code{\link{plot.posterior.density}}) can be added.
 188 ##' @name create.base.plot
 189 ##' @title Create Base Plot
 190 ##' @return empty ggplot object
 191 ##' @export
 192 ##' @author David LeBauer
 193 create.base.plot <- function() {
 194   base.plot <- ggplot()
 195   return(base.plot)
 196 }
 197 #==================================================================================================#
 198
 199
 200
 201
 202 ##--------------------------------------------------------------------------------------------------#
 203 ##' Add data to an existing plot or create a new one from \code{\link{create.base.plot}}
 204 ##'
 205 ##' Used to add raw data or summary statistics to the plot of a distribution.
 206 ##' The height of Y is arbitrary, and can be set to optimize visualization.
 207 ##' If SE estimates are available, tehse wil be plotted
 208 ##' @name plot.data
 209 ##' @title Add data to plot
 210 ##' @param trait.data data to be plotted
 211 ##' @param base.plot a ggplot object (grob),
 212 ##' created by \code{\link{create.base.plot}} if none provided
 213 ##' @param ymax maximum height of y
 214 ##' @seealso \code{\link{create.base.plot}}
 215 ##' @return updated plot object
 216 ##' @author David LeBauer
 217 ##' @export
 218 ##' @examples
 219 ##' \dontrun{plot.data(data.frame(Y = c(1, 2), se = c(1,2)), base.plot = NULL, ymax = 10)}
 220 plot.data <- function(trait.data, base.plot = NULL, ymax, color = 'black') {
 221   if(is.null(base.plot)) base.plot <- create.base.plot()
 222   n.pts <- nrow(trait.data)
 223   if(n.pts == 1){
 224     ymax <- ymax/8
 225   } else if (n.pts < 5) {
 226     ymax <- ymax / 4
 227   } else {
 228     ymax <- ymax / 2
 229   }
 230   y.pts <- seq(0, ymax, length.out = 1 + n.pts)[-1]
 231   plot.data <- data.frame(x = trait.data$Y,
 232                           y = y.pts,
 233                           se = trait.data$se,
 234                           control = !trait.data$trt == 1 & trait.data$ghs == 1)
 235   new.plot <- base.plot +
 236     geom_point(data = plot.data,
 237                aes(x = x, y = y,
 238                color = control)) +
 239                  geom_segment(data = plot.data,
 240                               aes(x = x - se, y = y, xend = x + se, yend = y,
 241                                   color = control)) +
 242                                     scale_color_manual(values = c('black', 'grey')) +
 243                                       opts(legend_position = "none")
 244   return(new.plot)
 245 }
 246 ##==================================================================================================#
 247
 248
 249 #--------------------------------------------------------------------------------------------------#
 250 ##' Add borders to .. content for \description{} (no empty lines) ..
 251 ##'
 252 ##' Has ggplot2 display only specified borders, e.g. ("L"-shaped) borders, rather than a rectangle or no border. Note that the order can be significant; for example, if you specify the L border option and then a theme, the theme settings will override the border option, so you need to specify the theme (if any) before the border option, as above.
 253 ##' @name theme_border
 254 ##' @title Theme border for plot
 255 ##' @param type
 256 ##' @param colour
 257 ##' @param size
 258 ##' @param linetype
 259 ##' @return adds borders to ggplot as a side effect
 260 ##' @author Rudolf Cardinal
 261 ##' @author \url{ggplot2 google group}{https://groups.google.com/forum/?fromgroups#!topic/ggplot2/-ZjRE2OL8lE}
 262 ##' @examples
 263 ##' \dontrun{
 264 ##' df = data.frame( x=c(1,2,3), y=c(4,5,6) )
 265 ##' ggplot(data=df, aes(x=x, y=y)) + geom_point() + theme_bw() +
 266 ##' opts(panel.border = theme_border(c("bottom","left")) )
 267 ##' ggplot(data=df, aes(x=x, y=y)) + geom_point() + theme_bw() +
 268 ##' opts(panel.border = theme_border(c("b","l")) )
 269 ##' }
 270 theme_border <- function(type = c("left", "right", "bottom", "top",
 271                            "none"), colour = "black", size = 1, linetype = 1) {
 272   type <- match.arg(type, several.ok=TRUE)
 273   structure(
 274             function(x = 0, y = 0, width = 1, height = 1, ...) {
 275               xlist <- c()
 276               ylist <- c()
 277               idlist <- c()
 278               if ("bottom" %in% type) { # bottom
 279                 xlist <- append(xlist, c(x, x+width))
 280                 ylist <- append(ylist, c(y, y))
 281                 idlist <- append(idlist, c(1,1))
 282               }
 283               if ("top" %in% type) { # top
 284                 xlist <- append(xlist, c(x, x+width))
 285                 ylist <- append(ylist, c(y+height, y+height))
 286                 idlist <- append(idlist, c(2,2))
 287               }
 288               if ("left" %in% type) { # left
 289                 xlist <- append(xlist, c(x, x))
 290                 ylist <- append(ylist, c(y, y+height))
 291                 idlist <- append(idlist, c(3,3))
 292               }
 293               if ("right" %in% type) { # right
 294                 xlist <- append(xlist, c(x+width, x+width))
 295                 ylist <- append(ylist, c(y, y+height))
 296                 idlist <- append(idlist, c(4,4))
 297               }
 298               polylineGrob(
 299                            x=xlist, y=ylist, id=idlist, ..., default.units = "npc",
 300                            gp=gpar(lwd=size, col=colour, lty=linetype),
 301                            )
 302             },
 303             class = "theme",
 304             type = "box",
 305             call = match.call()
 306             )
 307 }
 308 #==================================================================================================#
 309
 310
 311 ####################################################################################################
 312 ### EOF. End of R script file.
 313 ####################################################################################################