Heatmap.Rmd

---
title: "Heatmap for homegardens"
author: "Cory Whitney"
output: 
  html_document:
    toc: true
    toc_float: true
bibliography: packages.bib
---

<!-- Spelling -->
<!-- The ABC √ option (upper right on the Rmarkdown console)-->

<!-- Grammar -->
<!-- devtools::install_github("ropenscilabs/gramr") -->
<!-- run_grammar_checker("HighDimensionalData.rmd”) -->

<!-- Print pdf and word versions -->
<!-- rmarkdown::render("Heatmap.Rmd", output_format = "pdf_document") -->
<!-- rmarkdown::render("HighDimensionalData.Rmd", output_format = "word_document") -->

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
#packages in alphabetic order
library(bayesplot)
library(dplyr)
library(ethnobotanyR)
library(ggplot2)
library(ggpubr)
library(ggridges)
library(ggstance)
library(plyr)
library(RColorBrewer)
library(reshape)
library(tidyverse) #tidyverse includes a number of useful packages
```

```{r packages.bib, warning=FALSE, include = FALSE}
#Automatically write R package citation entries to a .bib file
knitr::write_bib(c(.packages(), 
                   'bayesplot', 
                   'decisionSupport',
                   'dplyr',
                   'ethnobotanyR',
                   'ggplot2', 
                   'ggridges', 
                   'plyr', 
                   'RColorBrewer',
                   'reshape', 
                   'tidyverse'), 'packages.bib')
```

Homegarden studies often attempt to compare a collection of related variables. Heatmaps and cluster diagrams can be useful for this. Here is a way to get an overview of the relationships by generating intergrated plots in the [R programming language](https://www.r-project.org/) [@R-base].

Needless to say any aims objectives and hypotheses should be determined before any data is collected. Data visualization is a good time to get a clear sense for how your data looks, but is not the time to start making up hypotheses about it.

## Heatmap function

We use a garden heatmap function `heatmap.function` after [obigriffith's heatmap.function](https://www.biostars.org/p/18211/). We use `grDevices` `terrain.colors` for filling the heatmap. Other colors can be chose from the `RColorBrewer` library [@R-RColorBrewer].

```{r heatmap.function, include=FALSE}
heatmap.function <- function(x,
                             row_values = TRUE,
                             cols_values = if (symm)
                               "row_values"
                             else
                               TRUE,
                             distfun = dist,
                             hclustfun = hclust,
                             dendrogram = c("both", "row", "column", "none"),
                             symm = FALSE,
                             scale = c("none", "row", "column"),
                             na.rm = TRUE,
                             revC = identical(cols_values, "row_values"),
                             add.expr,
                             breaks,
                             symbreaks = max(x < 0, na.rm = TRUE) ||
                               scale != "none",
                             col = "terrain.colors",
                             colsep,
                             rowsep,
                             sepcolor = "white",
                             sepwidth = c(0.05, 0.05),
                             cellnote,
                             notecex = 1,
                             notecol = "cyan",
                             na.color = par("bg"),
                             trace = c("none", "column", "row", "both"),
                             tracecol = "cyan",
                             hline = median(breaks),
                             vline = median(breaks),
                             linecol = tracecol,
                             margins = c(5, 5),
                             ColSideColors,
                             RowSideColors,
                             side.height.fraction = 0.3,
                             cexRow = 0.2 + 1 / log10(nr),
                             cexCol = 0.2 + 1 / log10(nc),
                             labRow = NULL,
                             labCol = NULL,
                             key = TRUE,
                             keysize = 1.5,
                             density.info = c("none", "histogram", "density"),
                             denscol = tracecol,
                             symkey = max(x < 0, na.rm = TRUE) ||
                               symbreaks,
                             densadj = 0.25,
                             main = NULL,
                             xlab = NULL,
                             ylab = NULL,
                             lmat = NULL,
                             lhei = NULL,
                             lwid = NULL,
                             ColSideColorsSize = 1,
                             RowSideColorsSize = 1,
                             KeyValueName = "Value",
                             ...) {
  invalid <- function (x) {
    if (missing(x) || is.null(x) || length(x) == 0)
      return(TRUE)
    if (is.list(x))
      return(all(sapply(x, invalid)))
    else if (is.vector(x))
      return(all(is.na(x)))
    else
      return(FALSE)
  }
  
  x <- as.matrix(x)
  scale01 <- function(x,
                      low = min(x),
                      high = max(x)) {
    x <- (x - low) / (high - low)
    x
  }
  retval <- list()
  scale <- if (symm && missing(scale))
    "none"
  else
    match.arg(scale)
  dendrogram <- match.arg(dendrogram)
  trace <- match.arg(trace)
  density.info <- match.arg(density.info)
  if (length(col) == 1 && is.character(col))
    col <- get(col, mode = "function")
  if (!missing(breaks) && (scale != "none"))
    warning(
      "Using scale=\"row\" or scale=\"column\" when breaks are",
      "specified can produce unpredictable results.",
      "Please consider using only one or the other."
    )
  if (is.null(row_values) || is.na(row_values))
    row_values <- FALSE
  if (is.null(cols_values) || is.na(cols_values))
    cols_values <- FALSE
  else if (cols_values == "row_values" && !isTRUE(row_values))
    cols_values <- FALSE
  if (length(di <- dim(x)) != 2 || !is.numeric(x))
    stop("`x' must be a numeric matrix")
  nr <- di[1]
  nc <- di[2]
  if (nr <= 1 || nc <= 1)
    stop("`x' must have at least 2 rows and 2 columns")
  if (!is.numeric(margins) || length(margins) != 2)
    stop("`margins' must be a numeric vector of length 2")
  if (missing(cellnote))
    cellnote <- matrix("", ncol = ncol(x), nrow = nrow(x))
  if (!inherits(row_values, "dendrogram")) {
    if (((!isTRUE(row_values)) || (is.null(row_values))) && (dendrogram %in%
                                                 c("both", "row"))) {
      if (is.logical(cols_values) && (cols_values))
        dendrogram <- "column"
      else
        dedrogram <- "none"
      warning(
        "Discrepancy: row_values is FALSE, while dendrogram is `",
        dendrogram,
        "'. Omitting row dendogram."
      )
    }
  }
  if (!inherits(cols_values, "dendrogram")) {
    if (((!isTRUE(cols_values)) || (is.null(cols_values))) && (dendrogram %in%
                                                 c("both", "column"))) {
      if (is.logical(row_values) && (row_values))
        dendrogram <- "row"
      else
        dendrogram <- "none"
      warning(
        "Discrepancy: cols_values is FALSE, while dendrogram is `",
        dendrogram,
        "'. Omitting column dendogram."
      )
    }
  }
  if (inherits(row_values, "dendrogram")) {
    ddr <- row_values
    rowInd <- order.dendrogram(ddr)
  }
  else if (is.integer(row_values)) {
    hcr <- hclustfun(distfun(x))
    ddr <- as.dendrogram(hcr)
    ddr <- reorder(ddr, row_values)
    rowInd <- order.dendrogram(ddr)
    if (nr != length(rowInd))
      stop("row dendrogram ordering gave index of wrong length")
  }
  else if (isTRUE(row_values)) {
    row_values <- rowMeans(x, na.rm = na.rm)
    hcr <- hclustfun(distfun(x))
    ddr <- as.dendrogram(hcr)
    ddr <- reorder(ddr, row_values)
    rowInd <- order.dendrogram(ddr)
    if (nr != length(rowInd))
      stop("row dendrogram ordering gave index of wrong length")
  }
  else {
    rowInd <- nr:1
  }
  if (inherits(cols_values, "dendrogram")) {
    ddc <- cols_values
    colInd <- order.dendrogram(ddc)
  }
  else if (identical(cols_values, "row_values")) {
    if (nr != nc)
      stop("cols_values = \"row_values\" but nrow(x) != ncol(x)")
    if (exists("ddr")) {
      ddc <- ddr
      colInd <- order.dendrogram(ddc)
    }
    else
      colInd <- rowInd
  }
  else if (is.integer(cols_values)) {
    hcc <- hclustfun(distfun(if (symm)
      x
      else
        t(x)))
    ddc <- as.dendrogram(hcc)
    ddc <- reorder(ddc, cols_values)
    colInd <- order.dendrogram(ddc)
    if (nc != length(colInd))
      stop("column dendrogram ordering gave index of wrong length")
  }
  else if (isTRUE(cols_values)) {
    cols_values <- colMeans(x, na.rm = na.rm)
    hcc <- hclustfun(distfun(if (symm)
      x
      else
        t(x)))
    ddc <- as.dendrogram(hcc)
    ddc <- reorder(ddc, cols_values)
    colInd <- order.dendrogram(ddc)
    if (nc != length(colInd))
      stop("column dendrogram ordering gave index of wrong length")
  }
  else {
    colInd <- 1:nc
  }
  retval$rowInd <- rowInd
  retval$colInd <- colInd
  retval$call <- match.call()
  x <- x[rowInd, colInd]
  x.unscaled <- x
  cellnote <- cellnote[rowInd, colInd]
  if (is.null(labRow))
    labRow <- if (is.null(rownames(x)))
      (1:nr)[rowInd]
  else
    rownames(x)
  else
    labRow <- labRow[rowInd]
  if (is.null(labCol))
    labCol <- if (is.null(colnames(x)))
      (1:nc)[colInd]
  else
    colnames(x)
  else
    labCol <- labCol[colInd]
  if (scale == "row") {
    retval$rowMeans <- rm <- rowMeans(x, na.rm = na.rm)
    x <- sweep(x, 1, rm)
    retval$rowSDs <- sx <- apply(x, 1, sd, na.rm = na.rm)
    x <- sweep(x, 1, sx, "/")
  }
  else if (scale == "column") {
    retval$colMeans <- rm <- colMeans(x, na.rm = na.rm)
    x <- sweep(x, 2, rm)
    retval$colSDs <- sx <- apply(x, 2, sd, na.rm = na.rm)
    x <- sweep(x, 2, sx, "/")
  }
  if (missing(breaks) || is.null(breaks) || length(breaks) < 1) {
    if (missing(col) || is.function(col))
      breaks <- 16
    else
      breaks <- length(col) + 1
  }
  if (length(breaks) == 1) {
    if (!symbreaks)
      breaks <-
        seq(min(x, na.rm = na.rm), max(x, na.rm = na.rm),
            length = breaks)
    else {
      extreme <- max(abs(x), na.rm = TRUE)
      breaks <- seq(-extreme, extreme, length = breaks)
    }
  }
  nbr <- length(breaks)
  ncol <- length(breaks) - 1
  if (class(col) == "function")
    col <- col(ncol)
  min.breaks <- min(breaks)
  max.breaks <- max(breaks)
  x[x < min.breaks] <- min.breaks
  x[x > max.breaks] <- max.breaks
  if (missing(lhei) || is.null(lhei))
    lhei <- c(keysize, 4)
  if (missing(lwid) || is.null(lwid))
    lwid <- c(keysize, 4)
  if (missing(lmat) || is.null(lmat)) {
    lmat <- rbind(4:3, 2:1)
    
    if (!missing(ColSideColors)) {
      #if (!is.matrix(ColSideColors))
      #stop("'ColSideColors' must be a matrix")
      if (!is.character(ColSideColors) ||
          nrow(ColSideColors) != nc)
        stop("'ColSideColors' must be a matrix of nrow(x) rows")
      lmat <- rbind(lmat[1,] + 1, c(NA, 1), lmat[2,] + 1)
      #lhei <- c(lhei[1], 0.2, lhei[2])
      lhei = c(lhei[1], side.height.fraction * ColSideColorsSize /
                 2, lhei[2])
    }
    
    if (!missing(RowSideColors)) {
      #if (!is.matrix(RowSideColors))
      #stop("'RowSideColors' must be a matrix")
      if (!is.character(RowSideColors) ||
          ncol(RowSideColors) != nr)
        stop("'RowSideColors' must be a matrix of ncol(x) columns")
      lmat <-
        cbind(lmat[, 1] + 1, c(rep(NA, nrow(lmat) - 1), 1), lmat[, 2] + 1)
      #lwid <- c(lwid[1], 0.2, lwid[2])
      lwid <-
        c(lwid[1], side.height.fraction * RowSideColorsSize / 2, lwid[2])
    }
    lmat[is.na(lmat)] <- 0
  }
  
  if (length(lhei) != nrow(lmat))
    stop("lhei must have length = nrow(lmat) = ", nrow(lmat))
  if (length(lwid) != ncol(lmat))
    stop("lwid must have length = ncol(lmat) =", ncol(lmat))
  op <- par(no.readonly = TRUE)
  on.exit(par(op))
  
  layout(lmat,
         widths = lwid,
         heights = lhei,
         respect = FALSE)
  
  if (!missing(RowSideColors)) {
    if (!is.matrix(RowSideColors)) {
      par(mar = c(margins[1], 0, 0, 0.5))
      image(rbind(1:nr), col = RowSideColors[rowInd], axes = FALSE)
    } else {
      par(mar = c(margins[1], 0, 0, 0.5))
      rsc = t(RowSideColors[, rowInd, drop = F])
      rsc.colors = matrix()
      rsc.names = names(table(rsc))
      rsc.i = 1
      for (rsc.name in rsc.names) {
        rsc.colors[rsc.i] = rsc.name
        rsc[rsc == rsc.name] = rsc.i
        rsc.i = rsc.i + 1
      }
      rsc = matrix(as.numeric(rsc), nrow = dim(rsc)[1])
      image(t(rsc), col = as.vector(rsc.colors), axes = FALSE)
      if (length(rownames(RowSideColors)) > 0) {
        axis(
          1,
          0:(dim(rsc)[2] - 1) / max(1, (dim(rsc)[2] - 1)),
          rownames(RowSideColors),
          las = 2,
          tick = FALSE
        )
      }
    }
  }
  
  if (!missing(ColSideColors)) {
    if (!is.matrix(ColSideColors)) {
      par(mar = c(0.5, 0, 0, margins[2]))
      image(cbind(1:nc), col = ColSideColors[colInd], axes = FALSE)
    } else {
      par(mar = c(0.5, 0, 0, margins[2]))
      csc = ColSideColors[colInd, , drop = F]
      csc.colors = matrix()
      csc.names = names(table(csc))
      csc.i = 1
      for (csc.name in csc.names) {
        csc.colors[csc.i] = csc.name
        csc[csc == csc.name] = csc.i
        csc.i = csc.i + 1
      }
      csc = matrix(as.numeric(csc), nrow = dim(csc)[1])
      image(csc, col = as.vector(csc.colors), axes = FALSE)
      if (length(colnames(ColSideColors)) > 0) {
        axis(
          2,
          0:(dim(csc)[2] - 1) / max(1, (dim(csc)[2] - 1)),
          colnames(ColSideColors),
          las = 2,
          tick = FALSE
        )
      }
    }
  }
  
  par(mar = c(margins[1], 0, 0, margins[2]))
  x <- t(x)
  cellnote <- t(cellnote)
  if (revC) {
    iy <- nr:1
    if (exists("ddr"))
      ddr <- ddr
    x <- x[, iy]
    cellnote <- cellnote[, iy]
  }
  else
    iy <- 1:nr
  image(
    1:nc,
    1:nr,
    x,
    xlim = 0.5 + c(0, nc),
    ylim = 0.5 + c(0, nr),
    axes = FALSE,
    xlab = "",
    ylab = "",
    col = col,
    breaks = breaks,
    ...
  )
  retval$carpet <- x
  if (exists("ddr"))
    retval$rowDendrogram <- ddr
  if (exists("ddc"))
    retval$colDendrogram <- ddc
  retval$breaks <- breaks
  retval$col <- col
  if (!invalid(na.color) & any(is.na(x))) {
    # load library(gplots)
    mmat <- ifelse(is.na(x), 1, NA)
    image(
      1:nc,
      1:nr,
      mmat,
      axes = FALSE,
      xlab = "",
      ylab = "",
      col = na.color,
      add = TRUE
    )
  }
  axis(
    1,
    1:nc,
    labels = labCol,
    las = 2,
    line = -0.5,
    tick = 0,
    cex.axis = cexCol
  )
  if (!is.null(xlab))
    mtext(xlab, side = 1, line = margins[1] - 1.25)
  axis(
    4,
    iy,
    labels = labRow,
    las = 2,
    line = -0.5,
    tick = 0,
    cex.axis = cexRow
  )
  if (!is.null(ylab))
    mtext(ylab, side = 4, line = margins[2] - 1.25)
  if (!missing(add.expr))
    eval(substitute(add.expr))
  if (!missing(colsep))
    for (csep in colsep)
      rect(
        xleft = csep + 0.5,
        ybottom = rep(0, length(csep)),
        xright = csep + 0.5 + sepwidth[1],
        ytop = rep(ncol(x) + 1, csep),
        lty = 1,
        lwd = 1,
        col = sepcolor,
        border = sepcolor
      )
  if (!missing(rowsep))
    for (rsep in rowsep)
      rect(
        xleft = 0,
        ybottom = (ncol(x) + 1 - rsep) - 0.5,
        xright = nrow(x) + 1,
        ytop = (ncol(x) + 1 - rsep) - 0.5 - sepwidth[2],
        lty = 1,
        lwd = 1,
        col = sepcolor,
        border = sepcolor
      )
  min.scale <- min(breaks)
  max.scale <- max(breaks)
  x.scaled <- scale01(t(x), min.scale, max.scale)
  if (trace %in% c("both", "column")) {
    retval$vline <- vline
    vline.vals <- scale01(vline, min.scale, max.scale)
    for (i in colInd) {
      if (!is.null(vline)) {
        abline(v = i - 0.5 + vline.vals,
               col = linecol,
               lty = 2)
      }
      xv <- rep(i, nrow(x.scaled)) + x.scaled[, i] - 0.5
      xv <- c(xv[1], xv)
      yv <- 1:length(xv) - 0.5
      lines(
        x = xv,
        y = yv,
        lwd = 1,
        col = tracecol,
        type = "s"
      )
    }
  }
  if (trace %in% c("both", "row")) {
    retval$hline <- hline
    hline.vals <- scale01(hline, min.scale, max.scale)
    for (i in rowInd) {
      if (!is.null(hline)) {
        abline(h = i + hline,
               col = linecol,
               lty = 2)
      }
      yv <- rep(i, ncol(x.scaled)) + x.scaled[i,] - 0.5
      yv <- c(yv[1], yv)
      xv <- length(yv):1 - 0.5
      lines(
        x = xv,
        y = yv,
        lwd = 1,
        col = tracecol,
        type = "s"
      )
    }
  }
  if (!missing(cellnote))
    text(
      x = c(row(cellnote)),
      y = c(col(cellnote)),
      labels = c(cellnote),
      col = notecol,
      cex = notecex
    )
  par(mar = c(margins[1], 0, 0, 0))
  if (dendrogram %in% c("both", "row")) {
    plot(
      ddr,
      horiz = TRUE,
      axes = FALSE,
      yaxs = "i",
      leaflab = "none"
    )
  }
  else
    plot.new()
  par(mar = c(0, 0, if (!is.null(main))
    5
    else
      0, margins[2]))
  if (dendrogram %in% c("both", "column")) {
    plot(ddc,
         axes = FALSE,
         xaxs = "i",
         leaflab = "none")
  }
  else
    plot.new()
  if (!is.null(main))
    title(main, cex.main = 1.5 * op[["cex.main"]])
  if (key) {
    par(mar = c(5, 4, 2, 1), cex = 0.75)
    tmpbreaks <- breaks
    if (symkey) {
      max.raw <- max(abs(c(x, breaks)), na.rm = TRUE)
      min.raw <- -max.raw
      tmpbreaks[1] <- -max(abs(x), na.rm = TRUE)
      tmpbreaks[length(tmpbreaks)] <-
        max(abs(x), na.rm = TRUE)
    }
    else {
      min.raw <- min(x, na.rm = TRUE)
      max.raw <- max(x, na.rm = TRUE)
    }
    
    z <- seq(min.raw, max.raw, length = length(col))
    image(
      z = matrix(z, ncol = 1),
      col = col,
      breaks = tmpbreaks,
      xaxt = "n",
      yaxt = "n"
    )
    par(usr = c(0, 1, 0, 1))
    lv <- pretty(breaks)
    xv <- scale01(as.numeric(lv), min.raw, max.raw)
    axis(1, at = xv, labels = lv)
    if (scale == "row")
      mtext(side = 1, "Row Z-Score", line = 2)
    else if (scale == "column")
      mtext(side = 1, "Column Z-Score", line = 2)
    else
      mtext(side = 1, KeyValueName, line = 2)
    if (density.info == "density") {
      dens <- density(x, adjust = densadj, na.rm = TRUE)
      omit <- dens$x < min(breaks) | dens$x > max(breaks)
      dens$x <- dens$x[-omit]
      dens$y <- dens$y[-omit]
      dens$x <- scale01(dens$x, min.raw, max.raw)
      lines(dens$x,
            dens$y / max(dens$y) * 0.95,
            col = denscol,
            lwd = 1)
      axis(2, at = pretty(dens$y) / max(dens$y) * 0.95, pretty(dens$y))
      title("Color Key\nand Density Plot")
      par(cex = 0.5)
      mtext(side = 2, "Density", line = 2)
    }
    else if (density.info == "histogram") {
      h <- hist(x, plot = FALSE, breaks = breaks)
      hx <- scale01(breaks, min.raw, max.raw)
      hy <- c(h$counts, h$counts[length(h$counts)])
      lines(
        hx,
        hy / max(hy) * 0.95,
        lwd = 1,
        type = "s",
        col = denscol
      )
      axis(2, at = pretty(hy) / max(hy) * 0.95, pretty(hy))
      title("Color Key\nand Histogram")
      par(cex = 0.5)
      mtext(side = 2, "Count", line = 2)
    }
    else
      title("Color Key")
  }
  else
    plot.new()
  retval$colorTable <-
    data.frame(
      low = retval$breaks[-length(retval$breaks)],
      high = retval$breaks[-1],
      color = retval$col
    )
  invisible(retval)
}
```

## Simple heatmap

Here is an example for running the `heatmap.function` on a data set with few parameters. 
 
```{r garden_matrix}

garden_matrix <- matrix(1:100, byrow = T, nrow = 10)

column_annotation <-
  sample(c("red", "blue", "green"), 10, replace = T)
column_annotation <- as.matrix(column_annotation)

colnames(column_annotation) <- c("Variable X")

row_annotation <- sample(c("red", "blue", "green"), 10, replace = T)
row_annotation <- as.matrix(t(row_annotation))

rownames(row_annotation) <- c("Variable Y")
```

Run the `heatmap.function` on the simple data.

```{r simple.heatmap.function}

heatmap.function(garden_matrix,
                 RowSideColors = row_annotation,
                 ColSideColors = column_annotation)

```

## Complex heatmap

Here is a more complex example with many parameters. First we create a dataset for demonstration purposes.

```{r plant_abundance_data}
plant_abundance = replicate(100, rnorm(20))
Garden_names = paste("Garden", letters[1:20], sep = "_")
Plant_ids = paste("Plant", c(1:100), sep = "_")
rownames(plant_abundance) = Garden_names
colnames(plant_abundance) = Plant_ids
```

Create color side bars to represent other variables to compare.

```{r complex_colors}
owner_gendercolors = sample(
  c("darkorchid", "darkred"),
  length(Garden_names),
  replace = TRUE,
  prob = NULL
)
salescolors = sample(c("green", "darkgreen"),
                     length(Garden_names),
                     replace = TRUE,
                     prob = NULL)
subtypecolors = sample(
  c("red", "blue", "cyan", "pink", "yellow", "green"),
  length(Plant_ids),
  replace = TRUE,
  prob = NULL
)
qualitycolors = sample(c("black", "white", "grey"),
                       length(Plant_ids),
                       replace = TRUE,
                       prob = NULL)
yieldscolors = sample(c("black", "white", "grey"),
                      length(Plant_ids),
                      replace = TRUE,
                      prob = NULL)
Productioncolors = sample(c("black", "white", "grey"),
                          length(Plant_ids),
                          replace = TRUE,
                          prob = NULL)
Usefulnesscolors = sample(c("black", "white", "grey"),
                          length(Plant_ids),
                          replace = TRUE,
                          prob = NULL)
Importancecolors = sample(c("black", "white", "grey"),
                          length(Plant_ids),
                          replace = TRUE,
                          prob = NULL)
NutritionContributioncolors = sample(c("black", "white", "grey"),
                                     length(Plant_ids),
                                     replace = TRUE,
                                     prob = NULL)
```

Create labels

```{r complex_labels}
rlab = t(cbind(owner_gendercolors, salescolors))

clab = cbind(
  subtypecolors,
  qualitycolors,
  yieldscolors,
  Productioncolors,
  Usefulnesscolors,
  Importancecolors,
  NutritionContributioncolors
)
```

Create row and column names.

```{r complex_col_names}

rownames(rlab) = c("Gender", "Economy")

colnames(clab) = c(
  "Subtype",
  "Quality",
  "Yield",
  "Production",
  "Usefulness",
  "Importance",
  "NutritionContribution"
)
```

Define `dist` and `hclust` functions for the heatmap using `stats` [@R-base]. 

```{r mydist_myclust}
mydist=function(c) {dist(c,method="euclidian")}
myclust=function(c) {hclust(c,method="average")}
```

Create heatmap using the `heatmap.function`. Colors for the map are from `grDevices` `terrain.colors` [@R-base].

```{r complex_heat_plot, out.width = "100%", dpi=300}

main_title = "Garden Plant Abundance"

par(cex.main = 1)
heatmap.function(
  plant_abundance,
  hclustfun = myclust,
  distfun = mydist,
  na.rm = TRUE,
  scale = "none",
  dendrogram = "both",
  margins = c(6, 12),
  row_values = TRUE,
  cols_values = TRUE,
  ColSideColors = clab,
  RowSideColors = rlab,
  symbreaks = FALSE,
  key = TRUE,
  symkey = FALSE,
  density.info = "none",
  trace = "none",
  main = main_title,
  labCol = FALSE,
  labRow = Garden_names,
  cexRow = 1,
  col = terrain.colors(n=10, alpha=1),
  ColSideColorsSize = 7,
  RowSideColorsSize = 2,
  KeyValueName = "Plant density"
)
legend(
  "topright",
  legend = c(
    "Ornamental",
    "Perennial",
    "Annual",
    "Shrub",
    "Tree",
    "Grass",
    "",
    "Poor (0)",
    "Fair (1)",
    "Good (2)",
    "",
    "Female",
    "Male",
    "",
    "Subsistence",
    "Sales"
  ),
  fill = c(
    "red",
    "blue",
    "cyan",
    "pink",
    "yellow",
    "green",
    "white",
    "pink",
    "grey",
    "black",
    "white",
    "darkorchid",
    "darkred",
    "white",
    "green",
    "darkgreen"
  ),
  border = FALSE,
  bty = "n",
  y.intersp = 0.7,
  cex = 0.7
)

```

# Make an image file

Create a png of the heatmap using `grDevices` [@R-base] and `heatmap.function`.

```{r offscreen}

png(file="garden_heatmap.png")
main_title="Garden Plant Abundance"

par(cex.main=1)

heatmap.function(
  plant_abundance,
  hclustfun = myclust,
  distfun = mydist,
  na.rm = TRUE,
  scale = "none",
  dendrogram = "both",
  margins = c(6, 12),
  row_values = TRUE,
  cols_values = TRUE,
  ColSideColors = clab,
  RowSideColors = rlab,
  symbreaks = FALSE,
  key = TRUE,
  symkey = FALSE,
  density.info = "none",
  trace = "none",
  main = main_title,
  labCol = FALSE,
  labRow = Garden_names,
  cexRow = 1,
  col = terrain.colors(75),
  ColSideColorsSize = 7,
  RowSideColorsSize = 2,
  KeyValueName = "Prob. Response"
)
legend(
  "topright",
  legend = c(
    "Ornamental",
    "Perennial",
    "Annual",
    "Shrub",
    "Tree",
    "Grass",
    "",
    "Poor (0)",
    "Fair (1)",
    "Good (2)",
    "",
    "Female",
    "Male",
    "",
    "Subsistence",
    "Sales"
  ),
  fill = c(
    "red",
    "blue",
    "cyan",
    "pink",
    "yellow",
    "green",
    "white",
    "white",
    "grey",
    "black",
    "white",
    "darkorchid",
    "darkred",
    "white",
    "green",
    "darkgreen"
  ),
  border = FALSE,
  bty = "n",
  y.intersp = 0.7,
  cex = 0.7
)
dev.off()
```

## Ethnobotany heatmap

Following the simple example above we run the `heatmap.function` for visualizing data typical collected in ethnobotany surveys. We use the `ethnobotanydata` data set from `ethnobotanyR` [@R-ethnobotanyR].
 
```{r ethno_matrix}

ethno_matrix <- as.matrix(dplyr::select(ethnobotanydata, -sp_name, -informant))

```

Run the `heatmap.function` on the simple data.

```{r ethno_heatmap}

heatmap.function(ethno_matrix)

```


# References