PageRenderTime 79ms CodeModel.GetById 24ms app.highlight 47ms RepoModel.GetById 1ms app.codeStats 1ms

/scripts/helper_functions.R

http://github.com/hpiwowar/alt-metrics_stats
R | 566 lines | 526 code | 15 blank | 25 comment | 213 complexity | fa3674c5578411e06bf11a56aee11afe MD5 | raw file
  1
  2#### HELPER FUNCTIONS
  3
  4# transform for count data due to
  5# http://pareonline.net/getvn.asp?v=8&n=6
  6# if no numbers between 0 and 1, not going to add a +1 offset
  7tr.0 = function(x) return(sqrt(x))
  8
  9# using sqrt with minimum value of 1, as per advice at
 10# http://www.webcitation.org/query?url=http%3A%2F%2Fpareonline.net%2Fgetvn.asp%3Fv%3D8%26n%3D6&date=2010-02-11
 11tr = function(x) return(sqrt(1 + x))
 12
 13log.tr = function(x) return(log(1 + x))
 14
 15#### Some of my variables are extracted by looking to see if MEDLINE records have a MeSH term.
 16# This is fine, but sometimes MEDLINE records are incomplete, they are in process or aren't on the list
 17# of journals indexed by MeSH indexers.  This means a lack of relevant MeSH term does not mean the MeSH
 18# term does not apply.  Thus, must replace these 0s with NAs to indicate the data is indeed missing
 19medline.values = function(raw_values, medline_status) {
 20    values = raw_values
 21    values[medline_status!="indexed for MEDLINE"] = NA
 22    return(values)
 23}
 24
 25print.thresh = function(x, thresh, decreasing=TRUE) {
 26    response = ""
 27    if (decreasing) {
 28        x.thresh = x[which(x>=thresh)]
 29    } else {
 30        x.thresh = x[which(x<=thresh)]
 31    }
 32    x.sorted = sort(x.thresh, decreasing)   
 33    for (ii in 1:length(x.sorted)) {
 34        response = cat(response, sprintf("%f %s\n", x.sorted[ii], names(x.sorted)[ii]))
 35    }
 36    if (length(response>0)) {
 37        print(response)
 38    }
 39}
 40
 41####Get scores:
 42# based on middle of factanal()
 43factor.scores.bartlett = function(x, fa.fit, na.action=NULL) {
 44    Lambda <- fa.fit$loadings
 45    z <- as.matrix(x)
 46    if (!is.numeric(z)) 
 47#        stop("factor analysis applies only to numerical variables")
 48        z = as.matrix(colwise(as.numeric)(x))
 49    zz <- scale(z, TRUE, TRUE)
 50    d <- 1/fa.fit$uniquenesses
 51    tmp <- t(Lambda * d)
 52    scores <- t(solve(tmp %*% Lambda, tmp %*% t(zz)))
 53    rownames(scores) <- rownames(z)
 54    colnames(scores) <- colnames(Lambda)
 55    if (!is.null(na.action)) 
 56        scores <- napredict(na.act, scores)
 57    scores
 58}
 59
 60# Extracted then modified from hetcor.data.frame from polycor
 61adjust.to.positive.definite = function(inputcor) {
 62    min.eigen = min(eigen(inputcor, only.values=TRUE)$values)
 63    if (min.eigen < 0){
 64        print("will try to make correlation matrix positive-definite")
 65        cor.corrected <- nearcor(inputcor)  # also could try nearPD
 66        if (!cor.corrected$converged) {
 67            stop("attempt to make correlation matrix positive-definite failed")
 68        }
 69        print("the correlation matrix has been adjusted to make it positive-definite")
 70        outputcor = cor.corrected$cor
 71        rownames(outputcor) <- rownames(inputcor)
 72        colnames(outputcor) <- colnames(inputcor)
 73    } else {
 74        print("Correlation matrix already positive-definite, so no adjustment needed")
 75        outputcor = inputcor
 76    }
 77    outputcor
 78}
 79
 80
 81# Based on hetcore.data.frame from polycor
 82# Modified slightly to call different correlation function for continuous correlations
 83# and print out updates
 84"hetcor.modified" <-
 85function(data, ML=FALSE, std.err=TRUE, use=c("complete.obs", "pairwise.complete.obs"),
 86  bins=4, pd=TRUE, type=c("pearson", "spearman"), ...){
 87  se.r <- function(r, n){
 88    rho <- r*(1 + (1 - r^2)/(2*(n - 3))) # approx. unbiased estimator
 89    v <- (((1 - rho^2)^2)/(n + 6))*(1 + (14 + 11*rho^2)/(2*(n + 6)))
 90    sqrt(v)
 91    }
 92  use <- match.arg(use)
 93  if (class(data) != "data.frame") stop("argument must be a data frame.")
 94  if (use == "complete.obs") data <- na.omit(data)
 95  p <- length(data)
 96  if (p < 2) stop("fewer than 2 variables.")
 97  R <- matrix(1, p, p)
 98  Type <- matrix("", p, p)
 99  SE <- matrix(0, p, p)
100  N <- matrix(0, p, p)
101  Test <- matrix(0, p, p)
102  diag(N) <- if (use == "complete.obs") nrow(data)
103             else sapply(data, function(x) sum(!is.na(x)))
104  for (i in 2:p) {
105    print(i)
106    for (j in 1:(i-1)){
107      x <- data[[i]]
108      y <- data[[j]]
109      if (inherits(x, c("numeric", "integer")) && inherits(y, c("numeric", "integer"))) {
110#         r <- cor(x, y, use="complete.obs")
111#         r <- rcorr(x, y, type="pearson")$r[1,2]
112         r <- rcorr(x, y, type=type)$r[1,2]
113#         Type[i, j] <- Type[j, i] <- "Pearson"
114#         Type[i, j] <- Type[j, i] <- "rcorr Pearson"
115		 if (type=="spearman")
116         	Type[i, j] <- Type[j, i] <- "rcorr Spearman"
117		 else
118      		Type[i, j] <- Type[j, i] <- "rcorr Pearson"
119         R[i, j] <- R[j, i] <- r
120         if (std.err) {
121           n <- sum(complete.cases(x, y))
122           SE[i, j] <- SE[j, i] <- se.r(r, n)
123           N[i, j] <- N[j, i] <- n
124           Test[i, j] <- pchisq(chisq(x, y, r, bins=bins), bins^2 - 2, lower.tail=FALSE)
125           }
126         }
127      else if (inherits(x, "factor") && inherits(y, "factor")) {
128         Type[i, j] <- Type[j, i] <- "Polychoric"
129         result <- polychor(x, y, ML=ML, std.err=std.err)
130         if (std.err){
131           n <- sum(complete.cases(x, y))
132           R[i, j] <- R[j, i] <- result$rho
133           SE[i, j] <- SE[j, i] <- sqrt(result$var[1,1])
134           N[i, j] <- N[j, i] <- n
135           Test[i, j] <- if (result$df > 0)
136                pchisq(result$chisq, result$df, lower.tail=FALSE)
137                else NA
138           }
139         else R[i, j] <- R[j, i] <- result
140         }
141       else {
142         if (inherits(x, "factor") && inherits(y, c("numeric", "integer")))
143           result <- polyserial(y, x, ML=ML, std.err=std.err, bins=bins)
144         else if (inherits(x, c("numeric", "integer")) && inherits(y, "factor"))
145           result <- polyserial(x, y, ML=ML, std.err=std.err, bins=bins)
146         else {
147             stop("columns must be numeric or factors.")
148             }
149         Type[i, j] <- Type[j, i] <- "Polyserial"
150         if (std.err){
151           n <- sum(complete.cases(x, y))
152           R[i, j] <- R[j, i] <- result$rho
153           SE[i, j] <- SE[j, i] <- sqrt(result$var[1,1])
154           N[i, j] <- N[j, i] <- n
155           Test[i, j] <- pchisq(result$chisq, result$df, lower.tail=FALSE)
156           }
157         else R[i, j] <- R[j, i] <- result
158         }
159       }
160     }
161   if (pd) {
162       if (min(eigen(R, only.values=TRUE)$values) < 0){
163            cor <- nearcor(R)
164            if (!cor$converged) stop("attempt to make correlation matrix positive-definite failed")
165            warning("the correlation matrix has been adjusted to make it positive-definite")
166            R <- cor$cor
167        }
168    }
169   rownames(R) <- colnames(R) <- names(data)
170   result <- list(correlations=R, type=Type, NA.method=use, ML=ML)
171   if (std.err) {
172     rownames(SE) <- colnames(SE) <- names(data)
173     rownames(N) <- colnames(N) <- names(N)
174     rownames(Test) <- colnames(Test) <- names(data)
175     result$std.errors <- SE
176     result$n <- if (use == "complete.obs") n else N
177     result$tests <- Test
178     }
179   class(result) <- "hetcor"
180   result
181   }
182
183
184# based on heatmap.2
185# Changed so that labels are on a different axis, close to clustering
186heatmap.3 = function (x, Rowv = TRUE, Colv = if (symm) "Rowv" else TRUE, 
187    distfun = dist, hclustfun = hclust, dendrogram = c("both", 
188        "row", "column", "none"), symm = FALSE, scale = c("none", 
189        "row", "column"), na.rm = TRUE, revC = identical(Colv, 
190        "Rowv"), add.expr, breaks, symbreaks = min(x < 0, na.rm = TRUE) || 
191        scale != "none", col = "heat.colors", colsep, rowsep, 
192    sepcolor = "white", sepwidth = c(0.05, 0.05), cellnote, notecex = 1, 
193    notecol = "cyan", na.color = par("bg"), trace = c("column", 
194        "row", "both", "none"), tracecol = "cyan", hline = median(breaks), 
195    vline = median(breaks), linecol = tracecol, margins = c(5, 
196        5), ColSideColors, RowSideColors, cexRow = 0.2 + 1/log10(nr), 
197    cexCol = 0.2 + 1/log10(nc), labRow = NULL, labCol = NULL, 
198    key = TRUE, keysize = 1.5, density.info = c("histogram", 
199        "density", "none"), denscol = tracecol, symkey = min(x < 
200        0, na.rm = TRUE) || symbreaks, densadj = 0.25, main = NULL, 
201    xlab = NULL, ylab = NULL, lmat = NULL, lhei = NULL, lwid = NULL, 
202    ...) 
203{
204    scale01 <- function(x, low = min(x), high = max(x)) {
205        x <- (x - low)/(high - low)
206        x
207    }
208    retval <- list()
209    scale <- if (symm && missing(scale)) 
210        "none"
211    else match.arg(scale)
212    dendrogram <- match.arg(dendrogram)
213    trace <- match.arg(trace)
214    density.info <- match.arg(density.info)
215    if (length(col) == 1 && is.character(col)) 
216        col <- get(col, mode = "function")
217    if (!missing(breaks) && (scale != "none")) 
218        warning("Using scale=\"row\" or scale=\"column\" when breaks are", 
219            "specified can produce unpredictable results.", "Please consider using only one or the other.")
220    if (is.null(Rowv) || is.na(Rowv)) 
221        Rowv <- FALSE
222    if (is.null(Colv) || is.na(Colv)) 
223        Colv <- FALSE
224    else if (Colv == "Rowv" && !isTRUE(Rowv)) 
225        Colv <- FALSE
226    if (length(di <- dim(x)) != 2 || !is.numeric(x)) 
227        stop("`x' must be a numeric matrix")
228    nr <- di[1]
229    nc <- di[2]
230    if (nr <= 1 || nc <= 1) 
231        stop("`x' must have at least 2 rows and 2 columns")
232    if (!is.numeric(margins) || length(margins) != 2) 
233        stop("`margins' must be a numeric vector of length 2")
234    if (missing(cellnote)) 
235        cellnote <- matrix("", ncol = ncol(x), nrow = nrow(x))
236    if (!inherits(Rowv, "dendrogram")) {
237        if (((!isTRUE(Rowv)) || (is.null(Rowv))) && (dendrogram %in% 
238            c("both", "row"))) {
239            if (is.logical(Colv) && (Colv)) 
240                dendrogram <- "column"
241            else dedrogram <- "none"
242            warning("Discrepancy: Rowv is FALSE, while dendrogram is `", 
243                dendrogram, "'. Omitting row dendogram.")
244        }
245    }
246    if (!inherits(Colv, "dendrogram")) {
247        if (((!isTRUE(Colv)) || (is.null(Colv))) && (dendrogram %in% 
248            c("both", "column"))) {
249            if (is.logical(Rowv) && (Rowv)) 
250                dendrogram <- "row"
251            else dendrogram <- "none"
252            warning("Discrepancy: Colv is FALSE, while dendrogram is `", 
253                dendrogram, "'. Omitting column dendogram.")
254        }
255    }
256    if (inherits(Rowv, "dendrogram")) {
257        ddr <- Rowv
258        rowInd <- order.dendrogram(ddr)
259    }
260    else if (is.integer(Rowv)) {
261        hcr <- hclustfun(distfun(x))
262        ddr <- as.dendrogram(hcr)
263        ddr <- reorder(ddr, Rowv)
264        rowInd <- order.dendrogram(ddr)
265        if (nr != length(rowInd)) 
266            stop("row dendrogram ordering gave index of wrong length")
267    }
268    else if (isTRUE(Rowv)) {
269        Rowv <- rowMeans(x, na.rm = na.rm)
270        hcr <- hclustfun(distfun(x))
271        ddr <- as.dendrogram(hcr)
272        ddr <- reorder(ddr, Rowv)
273        rowInd <- order.dendrogram(ddr)
274        if (nr != length(rowInd)) 
275            stop("row dendrogram ordering gave index of wrong length")
276    }
277    else {
278        rowInd <- nr:1
279    }
280    if (inherits(Colv, "dendrogram")) {
281        ddc <- Colv
282        colInd <- order.dendrogram(ddc)
283    }
284    else if (identical(Colv, "Rowv")) {
285        if (nr != nc) 
286            stop("Colv = \"Rowv\" but nrow(x) != ncol(x)")
287        if (exists("ddr")) {
288            ddc <- ddr
289            colInd <- order.dendrogram(ddc)
290        }
291        else colInd <- rowInd
292    }
293    else if (is.integer(Colv)) {
294        hcc <- hclustfun(distfun(if (symm) 
295            x
296        else t(x)))
297        ddc <- as.dendrogram(hcc)
298        ddc <- reorder(ddc, Colv)
299        colInd <- order.dendrogram(ddc)
300        if (nc != length(colInd)) 
301            stop("column dendrogram ordering gave index of wrong length")
302    }
303    else if (isTRUE(Colv)) {
304        Colv <- colMeans(x, na.rm = na.rm)
305        hcc <- hclustfun(distfun(if (symm) 
306            x
307        else t(x)))
308        ddc <- as.dendrogram(hcc)
309        ddc <- reorder(ddc, Colv)
310        colInd <- order.dendrogram(ddc)
311        if (nc != length(colInd)) 
312            stop("column dendrogram ordering gave index of wrong length")
313    }
314    else {
315        colInd <- 1:nc
316    }
317    retval$rowInd <- rowInd
318    retval$colInd <- colInd
319    retval$call <- match.call()
320    x <- x[rowInd, colInd]
321    x.unscaled <- x
322    cellnote <- cellnote[rowInd, colInd]
323    if (is.null(labRow)) 
324        labRow <- if (is.null(rownames(x))) 
325            (1:nr)[rowInd]
326        else rownames(x)
327    else labRow <- labRow[rowInd]
328    if (is.null(labCol)) 
329        labCol <- if (is.null(colnames(x))) 
330            (1:nc)[colInd]
331        else colnames(x)
332    else labCol <- labCol[colInd]
333    if (scale == "row") {
334        retval$rowMeans <- rm <- rowMeans(x, na.rm = na.rm)
335        x <- sweep(x, 1, rm)
336        retval$rowSDs <- sx <- apply(x, 1, sd, na.rm = na.rm)
337        x <- sweep(x, 1, sx, "/")
338    }
339    else if (scale == "column") {
340        retval$colMeans <- rm <- colMeans(x, na.rm = na.rm)
341        x <- sweep(x, 2, rm)
342        retval$colSDs <- sx <- apply(x, 2, sd, na.rm = na.rm)
343        x <- sweep(x, 2, sx, "/")
344    }
345    if (missing(breaks) || is.null(breaks) || length(breaks) < 
346        1) {
347        if (missing(col) || is.function(col)) 
348            breaks <- 16
349        else breaks <- length(col) + 1
350    }
351    if (length(breaks) == 1) {
352        if (!symbreaks) 
353            breaks <- seq(min(x, na.rm = na.rm), max(x, na.rm = na.rm), 
354                length = breaks)
355        else {
356            extreme <- max(abs(x), na.rm = TRUE)
357            breaks <- seq(-extreme, extreme, length = breaks)
358        }
359    }
360    nbr <- length(breaks)
361    ncol <- length(breaks) - 1
362    if (class(col) == "function") 
363        col <- col(ncol)
364    min.breaks <- min(breaks)
365    max.breaks <- max(breaks)
366    x[x < min.breaks] <- min.breaks
367    x[x > max.breaks] <- max.breaks
368    if (missing(lhei) || is.null(lhei)) 
369        lhei <- c(keysize, 4)
370    if (missing(lwid) || is.null(lwid)) 
371        lwid <- c(keysize, 4)
372    if (missing(lmat) || is.null(lmat)) {
373        lmat <- rbind(4:3, 2:1)
374        if (!missing(ColSideColors)) {
375            if (!is.character(ColSideColors) || length(ColSideColors) != 
376                nc) 
377                stop("'ColSideColors' must be a character vector of length ncol(x)")
378            lmat <- rbind(lmat[1, ] + 1, c(NA, 1), lmat[2, ] + 
379                1)
380            lhei <- c(lhei[1], 0.2, lhei[2])
381        }
382        if (!missing(RowSideColors)) {
383            if (!is.character(RowSideColors) || length(RowSideColors) != 
384                nr) 
385                stop("'RowSideColors' must be a character vector of length nrow(x)")
386            lmat <- cbind(lmat[, 1] + 1, c(rep(NA, nrow(lmat) - 
387                1), 1), lmat[, 2] + 1)
388            lwid <- c(lwid[1], 0.2, lwid[2])
389        }
390        lmat[is.na(lmat)] <- 0
391    }
392    if (length(lhei) != nrow(lmat)) 
393        stop("lhei must have length = nrow(lmat) = ", nrow(lmat))
394    if (length(lwid) != ncol(lmat)) 
395        stop("lwid must have length = ncol(lmat) =", ncol(lmat))
396    op <- par(no.readonly = TRUE)
397    on.exit(par(op))
398    layout(lmat, widths = lwid, heights = lhei, respect = FALSE)
399    if (!missing(RowSideColors)) {
400        par(mar = c(margins[1], 0, 0, 0.5))
401        image(rbind(1:nr), col = RowSideColors[rowInd], axes = FALSE)
402    }
403    if (!missing(ColSideColors)) {
404        par(mar = c(0.5, 0, 0, margins[2]))
405        image(cbind(1:nc), col = ColSideColors[colInd], axes = FALSE)
406    }
407    par(mar = c(margins[1], 0, 0, margins[2]))
408    x <- t(x)
409    cellnote <- t(cellnote)
410    if (revC) {
411        iy <- nr:1
412        if (exists("ddr")) 
413            ddr <- rev(ddr)
414        x <- x[, iy]
415        cellnote <- cellnote[, iy]
416    }
417    else iy <- 1:nr
418    image(1:nc, 1:nr, x, xlim = 0.5 + c(0, nc), ylim = 0.5 + 
419        c(0, nr), axes = FALSE, xlab = "", ylab = "", col = col, 
420        breaks = breaks, ...)
421    retval$carpet <- x
422    if (exists("ddr")) 
423        retval$rowDendrogram <- ddr
424    if (exists("ddc")) 
425        retval$colDendrogram <- ddc
426    retval$breaks <- breaks
427    retval$col <- col
428    if (!invalid(na.color) & any(is.na(x))) {
429        mmat <- ifelse(is.na(x), 1, NA)
430        image(1:nc, 1:nr, mmat, axes = FALSE, xlab = "", ylab = "", 
431            col = na.color, add = TRUE)
432    }
433    axis(1, 1:nc, labels = labCol, las = 2, line = -0.5, tick = 0, 
434        cex.axis = cexCol)
435    # Also put them on the other side    
436    axis(3, 1:nc, labels = labCol, las = 2, line = -0.5, tick = 0, 
437        cex.axis = cexCol)
438    if (!is.null(xlab)) 
439        mtext(xlab, side = 1, line = margins[1] - 1.25)
440    if (!is.null(xlab)) 
441        mtext(xlab, side = 3, line = margins[1] - 1.25)
442    # Also put them on the top axis
443    axis(2, iy, labels = labRow, las = 2, line = -0.5, tick = 0, 
444        cex.axis = cexRow)
445    axis(4, iy, labels = labRow, las = 2, line = -0.5, tick = 0, 
446        cex.axis = cexRow)
447    if (!is.null(ylab)) 
448        mtext(ylab, side = 2, line = margins[2] - 1.25)
449    if (!is.null(ylab)) 
450        mtext(ylab, side = 4, line = margins[2] - 1.25)
451    if (!missing(add.expr)) 
452        eval(substitute(add.expr))
453    if (!missing(colsep)) 
454        for (csep in colsep) rect(xleft = csep + 0.5, ybottom = rep(0, 
455            length(csep)), xright = csep + 0.5 + sepwidth[1], 
456            ytop = rep(ncol(x) + 1, csep), lty = 1, lwd = 1, 
457            col = sepcolor, border = sepcolor)
458    if (!missing(rowsep)) 
459        for (rsep in rowsep) rect(xleft = 0, ybottom = (ncol(x) + 
460            1 - rsep) - 0.5, xright = nrow(x) + 1, ytop = (ncol(x) + 
461            1 - rsep) - 0.5 - sepwidth[2], lty = 1, lwd = 1, 
462            col = sepcolor, border = sepcolor)
463    min.scale <- min(breaks)
464    max.scale <- max(breaks)
465    x.scaled <- scale01(t(x), min.scale, max.scale)
466    if (trace %in% c("both", "column")) {
467        retval$vline <- vline
468        vline.vals <- scale01(vline, min.scale, max.scale)
469        for (i in colInd) {
470            if (!is.null(vline)) {
471                abline(v = i - 0.5 + vline.vals, col = linecol, 
472                  lty = 2)
473            }
474            xv <- rep(i, nrow(x.scaled)) + x.scaled[, i] - 0.5
475            xv <- c(xv[1], xv)
476            yv <- 1:length(xv) - 0.5
477            lines(x = xv, y = yv, lwd = 1, col = tracecol, type = "s")
478        }
479    }
480    if (trace %in% c("both", "row")) {
481        retval$hline <- hline
482        hline.vals <- scale01(hline, min.scale, max.scale)
483        for (i in rowInd) {
484            if (!is.null(hline)) {
485                abline(h = i + hline, col = linecol, lty = 2)
486            }
487            yv <- rep(i, ncol(x.scaled)) + x.scaled[i, ] - 0.5
488            yv <- rev(c(yv[1], yv))
489            xv <- length(yv):1 - 0.5
490            lines(x = xv, y = yv, lwd = 1, col = tracecol, type = "s")
491        }
492    }
493    if (!missing(cellnote)) 
494        text(x = c(row(cellnote)), y = c(col(cellnote)), labels = c(cellnote), 
495            col = notecol, cex = notecex)
496    par(mar = c(margins[1], 0, 0, 0))
497    if (dendrogram %in% c("both", "row")) {
498        plot(ddr, horiz = TRUE, axes = FALSE, yaxs = "i", leaflab = "none")
499    }
500    else plot.new()
501    par(mar = c(0, 0, if (!is.null(main)) 5 else 0, margins[2]))
502    if (dendrogram %in% c("both", "column")) {
503        plot(ddc, axes = FALSE, xaxs = "i", leaflab = "none")
504    }
505    else plot.new()
506    if (!is.null(main)) 
507        title(main, cex.main = 1.5 * op[["cex.main"]])
508    if (key) {
509        par(mar = c(5, 4, 2, 1), cex = 0.75)
510        tmpbreaks <- breaks
511        if (symkey) {
512            max.raw <- max(abs(c(x, breaks)), na.rm = TRUE)
513            min.raw <- -max.raw
514            tmpbreaks[1] <- -max(abs(x))
515            tmpbreaks[length(tmpbreaks)] <- max(abs(x))
516        }
517        else {
518            min.raw <- min(x, na.rm = TRUE)
519            max.raw <- max(x, na.rm = TRUE)
520        }
521        z <- seq(min.raw, max.raw, length = length(col))
522        image(z = matrix(z, ncol = 1), col = col, breaks = tmpbreaks, 
523            xaxt = "n", yaxt = "n")
524        par(usr = c(0, 1, 0, 1))
525        lv <- pretty(breaks)
526        xv <- scale01(as.numeric(lv), min.raw, max.raw)
527        axis(1, at = xv, labels = lv)
528        if (scale == "row") 
529            mtext(side = 1, "Row Z-Score", line = 2)
530        else if (scale == "column") 
531            mtext(side = 1, "Column Z-Score", line = 2)
532        else mtext(side = 1, "Value", line = 2)
533        if (density.info == "density") {
534            dens <- density(x, adjust = densadj, na.rm = TRUE)
535            omit <- dens$x < min(breaks) | dens$x > max(breaks)
536            dens$x <- dens$x[-omit]
537            dens$y <- dens$y[-omit]
538            dens$x <- scale01(dens$x, min.raw, max.raw)
539            lines(dens$x, dens$y/max(dens$y) * 0.95, col = denscol, 
540                lwd = 1)
541            axis(2, at = pretty(dens$y)/max(dens$y) * 0.95, pretty(dens$y))
542            title("Color Key\nand Density Plot")
543            par(cex = 0.5)
544            mtext(side = 2, "Density", line = 2)
545        }
546        else if (density.info == "histogram") {
547            h <- hist(x, plot = FALSE, breaks = breaks)
548            hx <- scale01(breaks, min.raw, max.raw)
549            hy <- c(h$counts, h$counts[length(h$counts)])
550            lines(hx, hy/max(hy) * 0.95, lwd = 1, type = "s", 
551                col = denscol)
552            axis(2, at = pretty(hy)/max(hy) * 0.95, pretty(hy))
553            title("Color Key\nand Histogram")
554            par(cex = 0.5)
555            mtext(side = 2, "Count", line = 2)
556        }
557        else title("Color Key")
558    }
559    else plot.new()
560    retval$colorTable <- data.frame(low = retval$breaks[-length(retval$breaks)], 
561        high = retval$breaks[-1], color = retval$col)
562    invisible(retval)
563}
564
565
566