mod_data-viz.qmd

---
title: "Data Visualization & Exploration"
code-annotations: hover
---

## Overview

Data visualization is a fundamental part of working with data. Visualization can be only used in the final stages of a project to make figures for publication but it can also be hugely valuable for quality control and hypothesis development processes. This module focuses on the fundamentals of graph creation in an effort to empower you to apply those methods in the various contexts where you might find visualization to be helpful.

## Learning Objectives

After completing this module you will be able to: 

- <u>Explain</u> how data visualization can be used to explore data
- <u>Define</u> fundamental `ggplot2` vocabulary
- <u>Identify</u> appropriate graph types for given data type/distribution
- <u>Discuss</u> differences between presentation- and publication-quality graphs
- <u>Explain</u> how your graphs can be made more accessible

## Preparation

1. Each Synthesis fellow should download one data file identified for your group's project
2. _If you are a Mac user_, install [XQuartz](https://www.xquartz.org/)
3. _If you are an R user_, run the following code:

```{r install-pkgs}
#| eval: false

install.packages("librarian")
librarian::shelf(tidyverse, summarytools, datacleanr, lterdatasampler, supportR, cowplot)
```

## Networking Session

:::{.panel-tabset}
### 2024 Guests

We'll have two guests to kick off today's class. Each has been involved in synthesis as an early career researcher and each uses visualization in different ways to assess, clarify, and communicate their data and analyses.

- [Tim Ohlert](https://www.researchgate.net/scientific-contributions/Timothy-Ohlert-2172949124), Postdoctoral Researcher, Colorado State University; DroughtNet Coordinator
- [Kyle Cavanaugh](https://www.ioes.ucla.edu/person/kyle-cavanaugh/), Associate Professor, UCLA Institute of the Environment and Sustainability and the UCLA Geography Department

:::

## Data Visualization & The Synthesis Workflow

As shown in the graphic below, visualization can be valuable throughout the lifecycle of a synthesis project, albeit in different ways at different phases of a project.

<p align="center">
<img src="images/image_data-stages-with-viz.png" alt="Diagram depicting how raw data is transformed to cleaned data, then standardized data, and finally to published data products by a set of scripts between each 'type' of data" width="90%"/>
<figcaption>Diagram of data stages from raw data to published products. Credit: Margaret O'Brian & Li Kui & Sarah Elmendorf </figcaption>
</p>

## Visualization for Exploration

Exploratory data visualization is an important part of any scientific project. Before launching into analysis it is valuable to make some simple plots to scan the contents. These plots may reveal any number of issues, such as typos, sensor calibration problems or differences in the protocol over time.

These "fitness for use" visualizations are even more critical for synthesis projects. In synthesis, we are often repurposing publicly available datasets to answer questions that differ from the original motivations for data collection. As a result, the metadata included with a published dataset may be insufficient to assess whether the data are useful for your group's question. Datasets may not have been carefully quality-controlled prior to publication and could include any number of 'warts' that can complicate analyses or bias results. Some of these idiosyncrasies you may be able to anticipate in advance (e.g. spelling errors in taxonomy) and we encourage you to explicitly test for those and rectify them during the data harmonization process (see the [Data Wrangling module](https://lter.github.io/ssecr/mod_wrangle.html)). Others may come as a surprise.

During the early stages of a synthesis project, you will want to gain skill to quickly scan through large volumes of data. The figures you make will typically be for internal use only, and therefore have low emphasis on aesthetics.

### Exploratory Visualization Applications

Specific applications of exploratory data visualization include identifying:

1. Dataset coverage (temporal, spatial, taxonomic)
    - For example, the metadata might indicate a dataset covers the period 2010-2020. That could mean one data point in 2010 and one in 2020! This may not be useful for a time-series analysis.
2. Errors in metadata 
    - Do the units "make sense" with the figure? Typos in metadata do occur, so if you find yourself with elephants weighing only a few grams, it may be necessary to reach out to the dataset contact.
3. Differences in methodology
    - Do the data from sequential years, replicate sites, different providers generally fall into the same ranges or is there sensor drift or changes in protocols that need to be addressed?
    - A risk of synthesis projects is that you may find you are comparing apples to oranges across datasets, as the individual datasets included in your project were likely not collected in a coordinated fashion.
    - A benefit of synthesis projects is you will typically have large volumes of data, collected from many locations or timepoints. This data volume can be leveraged to give you a good idea of how your response variable looks at a 'typical' location as well as inform your gestalt sense of how much site-to-site, study-to-study, or year-to-year variability is expected. In our experience, where one particular dataset, or time period, strongly differs from the others, the most common root cause is differences in methodology that need to be addressed in the data harmonization process. 

In the data exploration stage you may find:

- Harmonization issues
    - Are all your datasets measured in units that can be converted to the same units?
    - If not, can you envision metrics (relative abundance? Effect size?) that would make datasets intercomparable?
- Some entire datasets cannot be used
- Parts of some datasets cannot be used
- Additional quality control is needed (e.g. filtering large outliers)

These steps are an important precursor to the data harmonization stage, where you will process the datasets you have selected into an analysis-ready format.

:::{.callout-note icon="false"}
#### Activity: Data Sleuth

In this activity, you'll play the role of data detective. You will have many potential datasets to look through. It is important to do it correctly, but you likely won't need or want to develop boutique code to examine each dataset, especially since some may be discarded after an initial pass.

As a project team, discuss the following points:

1. Decide on a structure for tracking results of exploratory data checks
    - Git issues? Additional columns in your team-data-inventory google sheet? Something else?
    - Draft a list of 'generic checks' you would want to apply to each dataset before inclusion in your synthesis
2. Use the `summarytools` and/or `datacleanr` packages to explore one exemplar dataset that you intend to include in your project
    - Discuss any issues you discover
    - Create a "to do" list for the exemplar dataset that details additional steps needed to make that dataset analysis ready (e.g. remove 1993 due to incomplete sampling, convert concentrations from mmols to mg/L, contact dataset providers to ask about anomalous values in April 2021)
    - Note we will work on skills to **implement** these steps in the [Data Wrangling module](https://lter.github.io/ssecr/mod_wrangle.html) in a few weeks.
    - Revise the list of 'generic checks' for remaining datasets as necessary
3. If you choose to save any images and/or code you used in your exploratory data visualization, decide on a naming convention and storage location
    - Will you add these files to your `.gitignore` or do you plan on committing them?
4. What additional plots would you ideally make that are not available through these generic tools?

::::{.panel-tabset}
##### `summarytools` Package

```{r data-sleuth_summarytools}
#| eval: false

# Load the library
library(summarytools)

# Load data
dataset_1 <- read_csv("your_file_name_here.csv")

# View the data in your Rstudio environment
summarytools::view(summarytools::dfSummary(dataset_1), footnote = NA) # <1>

# Alternatively,save the results for viewing later, or to share with your team
print(summarytools::dfSummary(dataset_1), footnote = NA,
      file = 'dataset_01_summary.html')
```
1. Careful! Use lowercase 'v' in the `view` function of the `summarytools` package

##### `datacleanr` Package

```{r data-sleuth_datacleanr}
#| eval: false

# Load the library
library(datacleanr)

# Load data
dataset_1 <- read_csv("your_file_name_here.csv")

# Launch the shiny app and view the data interactively
datacleanr::dcr_app(dataset_1)
```

::::

<br>

Both of these packages have extensive vignettes and online instructional materials. See [here](https://cran.r-project.org/web/packages/summarytools/vignettes/introduction.html) for one from `summarytools` and [here](https://the-hull.github.io/datacleanr/) for one from `datacleanr`.

:::

## Graphing with `ggplot2`

You may already be familiar with the `ggplot2` package in R but if you are not, it is a popular graphing library based on [The Grammar of Graphics](https://bookshop.org/p/books/the-grammar-of-graphics-leland-wilkinson/1518348?ean=9780387245447). Every ggplot is composed of four elements:

1. A 'core' `ggplot` function call
2. Aesthetics
3. Geometries
4. Theme

Note that the theme component may be implicit in some graphs because there is a suite of default theme elements that applies unless otherwise specified. 

This module will use example data to demonstrate these tools but as we work through these topics you should <u>feel free to substitute a dataset of your choosing</u>! If you don't have one in mind, you can use the example dataset shown in the code chunks throughout this module. This dataset comes from the [`lterdatasampler` R package](https://lter.github.io/lterdatasampler/) and the data are about fiddler crabs (_Minuca pugnax_) at the [Plum Island Ecosystems (PIE) LTER](https://pie-lter.mbl.edu/) site.

```{r data-prep}
#| message: false
#| warning: false

# Load needed libraries
library(tidyverse); library(lterdatasampler)

# Load the fiddler crab dataset
data(pie_crab)

# Check its structure
str(pie_crab)
```

With this dataset in hand, let's make a series of increasingly customized graphs to demonstrate some of the tools in `ggplot2`.

::::{.panel-tabset}
### 1. Starter Graph

Let's begin with a scatterplot of crab size on the Y-axis with latitude on the X. We'll forgo doing anything to the theme elements at this point to focus on the other three elements.

```{r gg-1}
#| fig-align: center
#| fig-width: 9
#| fig-height: 4

ggplot(data = pie_crab, mapping = aes(x = latitude, y = size, fill = site)) + # <1>
  geom_point(pch = 21, size = 2, alpha = 0.5) # <2>
```
1. We're defining both the data and the X/Y aesthetics in this top-level bit of the plot. Also, note that each line ends with a plus sign
2. Because we defined the data and aesthetics in the `ggplot()` function call above, this geometry can assume those mappings without re-specificying

### 2. Custom Theme

We can improve on this graph by tweaking theme elements to make it use fewer of the default settings.

```{r gg-2}
#| fig-align: center
#| fig-width: 9
#| fig-height: 4

ggplot(data = pie_crab, mapping = aes(x = latitude, y = size, fill = site)) +
  geom_point(pch = 21, size = 2, alpha = 0.5) +
  theme(legend.title = element_blank(), # <1>
        panel.background = element_blank(),
        axis.line = element_line(color = "black"))
```
1. All theme elements require these `element_...` helper functions. `element_blank` removes theme elements but otherwise you'll need to use the helper function that corresponds to the type of theme element (e.g., `element_text` for theme elements affecting graph text)

### 3. Multiple Geometries

We can further modify `ggplot2` graphs by adding _multiple_ geometries if you find it valuable to do so. Note however that <u>geometry order matters</u>! Geometries added later will be "in front of" those added earlier. Also, adding too much data to a plot will begin to make it difficult for others to understand the central take-away of the graph so you may want to be careful about the level of information density in each graph. Let's add boxplots behind the points to characterize the distribution of points more quantitatively.

```{r gg-3}
#| fig-align: center
#| fig-width: 9
#| fig-height: 4

ggplot(data = pie_crab, mapping = aes(x = latitude, y = size, fill = site)) +
  geom_boxplot(pch = 21) + # <1>
  geom_point(pch = 21, size = 2, alpha = 0.5) +
  theme(legend.title = element_blank(), 
        panel.background = element_blank(),
        axis.line = element_line(color = "black"))
```
1. By putting the boxplot geometry first we ensure that it doesn't cover up the points that overlap with the 'box' part of each boxplot

### 4. Multiple Datasets

`ggplot2` also supports adding more than one data object to the same graph! While this module doesn't cover map creation, maps are a common example of a graph with more than one data object. Another common use would be to include both the full dataset and some summarized facet of it in the same plot.

Let's calculate some summary statistics of crab size to include that in our plot.

```{r gg-4-prep}
#| message: false

# Load the supportR library
library(supportR)

# Summarize crab size within latitude groups
crab_summary <- supportR::summary_table(data = pie_crab, groups = c("site", "latitude"),
                                        response = "size", drop_na = TRUE)

# Check the structure
str(crab_summary)
```

With this data object in-hand, we can make a graph that includes both this and the original, unsummarized crab data. To better focus on the 'multiple data objects' bit of this example we'll pare down on the actual graph code.

```{r gg-4}
#| fig-align: center
#| fig-width: 9
#| fig-height: 4

ggplot() + # <1>
  geom_point(pie_crab, mapping = aes(x = latitude, y = size, fill = site),
             pch = 21, size = 2, alpha = 0.2) + 
  geom_errorbar(crab_summary, mapping = aes(x = latitude, # <2>
                                            ymax = mean + std_error,
                                            ymin = mean - std_error),
                width = 0.2) +
  geom_point(crab_summary, mapping = aes(x = latitude, y = mean, fill = site),
             pch = 23, size = 3) + 
  theme(legend.title = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(color = "black"))
```
1. If you want multiple data objects in the same `ggplot2` graph you need to leave this top level `ggplot()` call _empty!_ Otherwise you'll get weird errors with aesthetics later in the graph
2. This geometry adds the error bars and it's important that we add it before the summarized data points themselves if we want the error bars to be 'behind' their respective points

::::

:::{.callout-note icon="false"}
#### Activity: Graph Creation (P1)

In a script, attempt the following with one of either yours or your group's datasets:

- Make a graph using `ggplot2`
    - Include at least one geometry
    - Include at least one aesthetic (beyond X/Y axes)
    - Modify at least one theme element from the default

:::

## Streamlining Graph Aesthetics

Synthesis projects often generate an entire network of inter-related papers. Ensuring that all graphs across papers from a given team have a similar "feel" is a nice way of implying a certain standard of robustness for all of your group's projects. However, copy/pasting the theme elements of your graphs can (A) be cumbersome to do even once and (B) needs to be re-done every time you make a change anywhere. Fortunately, there is a better way!

`ggplot2` supports adding theme elements to an object that can then be reused as needed elsewhere. This is the same theory behind wrapping repeated operations into custom functions.

```{r std-theme}
#| fig-align: center
#| fig-width: 9
#| fig-height: 4

# Define core theme elements
theme_synthesis <- theme(legend.position = "none",
                         panel.background = element_blank(),
                         axis.line = element_line(color = "black"),
                         axis.text = element_text(size = 13)) # <1>

# Create a graph
ggplot(pie_crab, aes(y = water_temp, x = air_temp, color = size, size = size)) +
  geom_point() +
  theme_synthesis +
  theme(legend.position = "right") # <2>
```
1. This theme element controls the text on the tick marks. `axis.title` controls the text in the _labels_ of the axes
2. As a bonus, subsequent uses of `theme()` will replace defaults defined in your earlier theme object. So, you can design a set of theme elements that are _usually_ appropriate and then easily change just some of them as needed

:::{.callout-note icon="false"}
#### Activity: Graph Creation (P2)

In a script, attempt the following:

- Remove all theme edits from the graph you made in the preceding activity and assign them to a separate object
    - Then add that object to your graph
- Make a second (different) graph and add your consolidated theme object to that graph as well

:::

## Multi-Panel Graphs

It is sometimes the case that you want to make a single graph file that has multiple panels. For many of us, we might default to creating the separate graphs that we want, exporting them, and then using software like Microsoft PowerPoint to stitch those panels into the single image we had in mind from the start. However, as all of us who have used this method know, this is hugely cumbersome when your advisor/committee/reviewers ask for edits and you now have to redo all of the manual work behind your multi-panel graph. 

Fortunately, there are two nice entirely scripted alternatives that you might consider: **Faceted graphs** and **Plot grids**. See below for more information on both.

:::{.panel-tabset}
### Facets

In a faceted graph, <u>every panel of the graph has the same aesthetics</u>. These are often used when you want to show the relationship between two (or more) variables but separated by some other variable. In synthesis work, you might show the relationship between your core response and explanatory variables but facet by the original study. This would leave you with one panel per study where each would show the relationship only at that particular study.

Let's check out an example.

```{r facet-1}
#| fig-align: center
#| fig-width: 6
#| fig-height: 6

ggplot(pie_crab, aes(x = date, y = size, color = site))+
  geom_point(size = 2) +
  facet_wrap(. ~ site) + # <1>
  theme_bw() +
  theme(legend.position = "none") # <2>
```
1. This is a `ggplot2` function that assumes you want panels laid out in a regular grid. There are other `facet_...` alternatives that let you specify row versus column arrangement. You could also facet by multiple variables by putting something to the left of the tilde
2. We can remove the legend because the site names are in the facet titles in the gray boxes

### Plot Grids

In a plot grid, <u>each panel is completely independent of all others</u>. These are often used in publications where you want to highlight several _different_ relationships that have some thematic connection. In synthesis work, your hypotheses may be more complicated than in primary research and such a plot grid would then be necessary to put all visual evidence for a hypothesis in the same location. On a practical note, plot grids are also a common way of circumventing figure number limits enforced by journals.

Let's check out an example that relies on the `cowplot` library.

```{r grid-1}
#| message: false
#| fig-align: center
#| fig-width: 9
#| fig-height: 5

# Load a needed library
library(cowplot)

# Create the first graph
crab_p1 <- ggplot(pie_crab, aes(x = site, y = size, fill = site)) + # <1>
  geom_violin() +
  coord_flip() + # <2>
  theme_bw() +
  theme(legend.position = "none")

# Create the second
crab_p2 <- ggplot(pie_crab, aes(x = air_temp, y = water_temp)) +
  geom_errorbar(aes(ymax = water_temp + water_temp_sd, ymin = water_temp - water_temp_sd),
                 width = 0.1) +
  geom_errorbarh(aes(xmax = air_temp + air_temp_sd, xmin = air_temp - air_temp_sd), # <3>
                 width = 0.1) +
  geom_point(aes(fill = site), pch = 23, size = 3) +
  theme_bw()

# Assemble into a plot grid
cowplot::plot_grid(crab_p1, crab_p2, labels = "AUTO", nrow = 1) # <4>
```
1. Note that we're assigning these graphs to objects!
2. This is a handy function for flipping X and Y axes without re-mapping the aesthetics
3. This geometry is responsible for _horizontal_ error bars (note the "h" at the end of the function name)
4. The `labels = "AUTO"` argument means that each panel of the plot grid gets the next sequential capital letter. You could also substitute that for a vector with labels of your choosing
:::

:::{.callout-note icon="false"}
#### Activity: Graph Creation (P3)

In a script, attempt the following:

- Assemble the two graphs you made in the preceding two activities into the appropriate type of multi-panel graph

:::

## Accessibility Considerations

After you've made the graphs you need, it is good practice to revisit them with to ensure that they are as accessible as possible. You can of course also do this during the graph construction process but it is sometimes less onerous to tackle as a penultimate step in the figure creation process. There are many facets to accessibility and we've tried to cover just a few of them below.

### Color Choice

One of the more well-known facets of accessibility in data visualization is choosing colors that are "colorblind safe". Such palettes still create distinctive colors for those with various forms of color blindness (e.g., deuteranomoly, protanomaly, etc.). The classic red-green heatmap for instance is very colorblind unsafe in that people with some forms of colorblindness cannot distinguish between those colors (hence the rise of the yellow-blue heatmap in recent years). Unforunately, the `ggplot2` default rainbow palette--while nice for exploratory purposes--_is not_ colorlbind sfae.

Some websites (such as [colorbewer2.org](https://colorbrewer2.org/#type=sequential&scheme=YlGnBu&n=9)) include a simple checkbox for colorblindness safety which automatically limits the listed options to those that are colorblind safe. Alternately, you could use a browser plug-in (such as [Let's get color blind](https://chromewebstore.google.com/detail/lets-get-color-blind/bkdgdianpkfahpkmphgehigalpighjck) on Google Chrome) to simulate colorblindness on a particular page.

One extreme approach you could take is to dodge this issue entirely and format your graphs such that color either isn't used at all or only conveys information that is also conveyed in another graph aesthetic. We don't necessarily recommend this as color--when the palette is chosen correctly--can be a really nice way of making information-dense graphs more informative and easily-navigable by viewers.

### Multiple Modalities

Related to the color conversation is the value of mapping multiple aesthetics to the same variable. By presenting information in multiple ways--even if that seems redundant--you enable a wider audience to gain an intuitive sense of what you're trying to display.

```{r multi-modal}
#| fig-align: center
#| fig-width: 9
#| fig-height: 4

ggplot(data = pie_crab, mapping = aes(x = latitude, y = size, 
                                      fill = site, shape = site)) + # <1>
  geom_jitter(size = 2, width = 0.1, alpha = 0.6) +  
  scale_shape_manual(values = c(21:25, 21:25, 21:23)) + # <2>
  theme_bw() +
  theme(legend.title = element_blank())
```
1. In this graph we're mapping both the fill and shape aesthetics to site
2. This is a little cumbersome but there are only five 'fill-able' shapes in R so we need to reuse some of them to have a unique one for each site. Using fill-able shapes is nice because you get a crisp black border around each point. See `?pch` for all available shapes

In the above graph, even though the rainbow palette is not ideal for reasons mentioned earlier, it is now much easier to tell the difference between sites with similar colors. For instance, "NB", "NIB", and "PIE" are all shades of light blue/teal. Now that they have unique shapes it is dramatically easier to look at the graph and identify which points correspond to which site.


:::{.callout-warning icon="false"}
#### Discussion: Graph Accessibility

With a group discuss (some of) the following questions:

- What are other facets of accessibility that you think are important to consider when making data visualizations?
- What changes do you make to your graphs to increase accessibility?
    - What changes _could_ you make going forward?

:::


### Presentation vs. Publication

One final element of accessibility to consider is the difference between a '_presentation_-quality' graph and a '_publication_-quality' one. While it may be tempting to create a single version of a given graph and use it in both contexts that is likely to be less effective in helping you to get your point across than making small tweaks to two separate versions of what is otherwise the same graph.

:::{.panel-tabset}
### Presentation-Focused

**Do:**

- Increase size of text/points **greatly**
    - If possible, sit in the back row of the room where you'll present and look at your graphs from there
- _Consider_ adding graph elements that highlight certain graph regions
- Present summarized data (increases focus on big-picture trends and avoids discussion of minutiae)
- Map multiple aesthetics to the same variables

**Don't:**

- Use technical language / jargon
- Include _unnecessary_ background elements
- Use multi-panel graphs (either faceted or plot grid)
    - If you have multiple graph panels, put each on its own slide!

```{r talk-graph}
#| fig-align: center
#| fig-width: 9
#| fig-height: 4

ggplot(crab_summary, aes(x = latitude, y = mean, 
                         shape = reorder(site, latitude), # <1>
                         fill = reorder(site, latitude))) +
  geom_vline(xintercept = 36.5, color = "black", linetype = 1) +
  geom_vline(xintercept = 41.5, color = "black", linetype = 2) + # <2>
  geom_errorbar(mapping = aes(ymax = mean + std_error, ymin = mean - std_error),
                width = 0.2) +
  geom_point(size = 4) + 
  scale_shape_manual(values = c(21:25, 21:25, 21:23)) +
  labs(x = "Latitude", y = "Mean Crab Size (mm)") + # <3>
  theme(legend.title = element_blank(),
        axis.line = element_line(color = "black"),
        panel.background = element_blank(),
        axis.title = element_text(size = 17),
        axis.text = element_text(size = 15))
```
1. We can use the `reorder` function to make the order of sites in the legend (from top to bottom) match the order of sites in the graph (from left to right)
2. Adding vertical lines at particular parts in the graph can make comparisons within the same graph easier
3. `labs` lets us customize the title and label text of a graph

### Publication-Focused

**Do:**

- Increase size of text/points **slightly**
    - You want to be legible but you can more safely assume that many readers will be able to increase the zoom of their browser window if needed
- Present un-summarized data (with or without summarized points included)
    - Many reviewers will want to get a sense for the "real" data so you should include unsummarized values wherever possible
- Use multi-panel graphs
    - If multiple graphs "tell a story" together, then they should be included in the same file!
- Map multiple aesthetics to the same variables
- If publishing in a journal available in print, check to make sure your graph still makes sense in grayscale
    - There are nice browser plug-ins (like [Grayscale the Web](https://chromewebstore.google.com/detail/grayscale-the-web-save-si/mblmpdpfppogibmoobibfannckeeleag) for Google Chrome) for this too

**Don't:**

- Include _unnecessary_ background elements
- Add graph elements that highlight certain graph regions
    - You can--and should--lean more heavily on the text of your publication to discuss particular areas of a graph

```{r pub-graph}
#| fig-align: center
#| fig-width: 9
#| fig-height: 4

ggplot() +
  geom_point(pie_crab, mapping = aes(x = latitude, y = size,
                                     color = reorder(site, latitude)),
             pch = 19, size = 1, alpha = 0.3) +
  geom_errorbar(crab_summary, mapping = aes(x = latitude, y = mean, 
                                  ymax = mean + std_error, 
                                  ymin = mean - std_error),
                width = 0.2) +
  geom_point(crab_summary, mapping = aes(x = latitude, y = mean, 
                           shape = reorder(site, latitude),
                           fill = reorder(site, latitude)),
            size = 4) +
  scale_shape_manual(values = c(21:25, 21:25, 21:23)) +
  labs(x = "Latitude", y = "Mean Crab Carapace Width (mm)") + # <1>
  theme(legend.title = element_blank(),
        axis.line = element_line(color = "black"),
        panel.background = element_blank(),
        axis.title = element_text(size = 15),
        axis.text = element_text(size = 13))
```
1. Here we are using a reasonable amount of technical language

### Other Considerations

Some other factors you might consider _regardless of where the graphs will be embedded_ include:

- **White Background**. Ensure figures have a plain, white background for clarity and compatibility with journal formats.
- **High Resolution**. Use a resolution of at least 300 dpi for print quality. Journals often specify the minimum dpi required.
- **Bounding Box and Borders**. Add a bounding box or border if it enhances clarity, but avoid excessive framing unless necessary to separate elements clearly.
- **Clear Axis Labels**. Label axes with clear, concise descriptions, including units of measurement (e.g., "Temperature (°C)"). Use readable font sizes that remain legible when scaled.
- **Consistent Font Style and Size**. Use a uniform font style (e.g., Arial, Helvetica) across all figures and a size that is readable but not overwhelming (typically 8–12 points).
- **Color Scheme**. Choose a color palette that remains clear in both color and grayscale. Use distinct colors for different categories or groups, and avoid colors that may be difficult for colorblind readers to differentiate (e.g., red-green combinations).
- **Legend Placement**. Place legends within the figure space if possible, ensuring they don't overlap data or distract from the main content. Keep legends concise.
- **Minimal Gridlines**. Use minimal and subtle gridlines for reference, but avoid heavy or cluttered lines that may distract from the data.
- **Error Bars and Statistical Indicators**. Add error bars, confidence intervals, or statistical significance markers as needed to represent variability and support interpretation.
- **Descriptive Figure Caption**. Include a detailed caption that summarizes the figure's purpose, data source, and any essential methods or abbreviations. Captions should be self-contained to ensure figures are understandable independently.

:::

## Code Demo: Post-Harmonization Visualization

After harmonizing your data, you'll want to generate one last set of 'sanity check' plots to make sure (1) you have interpreted the metadata correctly (2) you haven't made any obvious errors in the harmonization and (3) your data are ready for analysis. Nothing is less fun than finding out your analytical results are due to an error in the underlying data.

The following is a multi-part code demonstration of three common post-harmonization uses of visualization. In addition to being useful graphs, there is also example code on how to export multiple panels of graphs into separate pages of a PDF which can be really helpful when reviewing exploratory visualizations as a group (without needing to scroll through a ton of separate graph files).

### Additional Needed Packages

If you'd like to follow along with the code chunks included throughout this demo, you'll need to install the following packages:

```{r load-pkgs}
#| message: false

## install.packages("librarian")
librarian::shelf(tidyverse, scales, ggforce, slider)
```

The three sets of plots below encompass many of the most common data structures
we have encountered types in ecological synthesis projects. These include 
quantitative measurements collected over many sites, taxonomic data collected
over many sites, and seasonal time series data.

::: panel-tabset
### Graph _All_ Numeric Variables

It can be helpful to visualize all numeric variables in your dataset, grouped by site (or dataset source) to check that the data have been homogenized correctly. As an example, we'll use a 2019 dataset on lake water quality, chemistry, and zooplankton community composition near the [Niwot Ridge](https://nwt.lternet.edu/) LTER. The dataset is a survey of 16 high alpine lakes and has structure similar to one that might be included in a multi-site synthesis. For more information on these data, check out [the data package on EDI](https://portal.edirepository.org/nis/mapbrowse?packageid=knb-lter-nwt.12.1). 

```{r demo_all-num-vars_data}
# Read in data
green_biochem <- read.csv(file = file.path("data", "green-lakes_water-chem-zooplank.csv")) %>%  # <1>
  dplyr::mutate(date = as.Date(date))

# Check structure
str(green_biochem)
```
1. Note that you could also read in this data directly from EDI. See ~line 31 of [this script](https://github.com/lter/ssecr/blob/main/scripts/prep-data_data-viz-demo.R) for a syntax example

Once we have the data file, we can programmatically identify all columns that R knows to be numeric.

```{r demo_all-num-vars_numcols}
# determine which columns are numeric in green_biochem
numcols <- green_biochem %>%
  dplyr::select(dplyr::where(~ is.numeric(.x) == TRUE)) %>% # <1>
  names(.) %>% 
  sort(.)

# Check that out
numcols # <2>
```
1. The tilde (`~`) is allowing us to evaluate each column against this conditional
2. You may notice that these columns all have `"num"` next to them in their structure check. The scripted method is _dramatically_ faster and more reproducible than writing these names down by hand

Now that we have our data and a vector of numeric column names, we can generate a multi-page PDF of scatterplots where each page is specific to a numeric variable and each graph panel within a given page reflects the time series at each site.

```{r demo_all-num-vars_viz-code-fake}
#| eval: false

# Open PDF 'device'
grDevices::pdf(file = file.path("qc_all_numeric.pdf")) # <1>

# Loop across numeric variables
for (var in numcols) {
  
  # Create a set of graphs for one variable
  myplot <- ggplot(green_biochem, aes(x = date, y = .data[[var]])) +
    geom_point(alpha = 0.5) + # <2>
    facet_wrap(. ~ local_site)
  
  # Print that variable
  print(myplot)
}

# Close the device
dev.off() # <3>
```
1. This function tells R that the following code should be saved as a PDF
2. A scatterplot may not be the best tool for your data; adjust appropriately
3. This function (when used after a 'device' function like `grDevices::pdf`) tells R when to stop adding things to the PDF and actually save it

The first page of the resulting plot should look something like the following, with each page having the same content but a different variable on the Y axis.

```{r demo_all-num-vars_viz-code-real}
#| fig-align: center
#| fig-width: 8
#| fig-height: 8
#| echo: false
#| warning: false

# Grab just one variable
var <- sort(numcols)[1]

# Graph it
ggplot(green_biochem, aes(x = date, y = .data[[var]])) +
  geom_point(alpha = 0.5) +
  facet_wrap(~local_site, ncol = 3)
```

### Taxonomic Consistency

Taxonomic time series can be tricky to work with due to inconsistencies in nomenclature and/or sampling effort. In particular, 'pseudoturnover' where one species 'disappears' with or without the simultaneous 'appearance' of another taxa can be indicative of either true extinctions, or changes in species names, or changes in methodology that cause particular taxa not to be detected. A second complication is that taxonomic data are often archived as 'presence-only' so it is necessary to _infer_ the absences based on sampling methodology and add them to your dataset before analysis.

While there are doubtless many field-collected datasets that have this issue, we've elected to simulate data so that we can emphasize the visualization elements of this problem while avoiding the "noise" typical of real data. This simulation is not necessarily vital to the visualization so we've left it out of the following demo. _However_, if that is of interest to you, see [this script](https://github.com/lter/ssecr/blob/main/scripts/prep-data_data-viz-demo.R)--in particular \~line 41 through \~80.

A workflow for establishing taxonomic consistency and plotting the results is included below.

```{r demo_tax-consist_data}
# Read in data
taxa_df <- read.csv(file.path("data", "simulated-taxa-df.csv"))

# Check structure
str(taxa_df)
```

First, we'll define units of sampling (year, plot and taxon) and 'pad out' the zeros. In this example, we have only added zeroes for taxa-plot-year combinations where that taxa is present in at least one year at a given plot. Again, this zero-padding is prerequisite to the visualization but not necessarily part of it so see \~lines 84-117 of the [prep script](https://github.com/lter/ssecr/blob/main/scripts/prep-data_data-viz-demo.R) if that process is of interest.

```{r  demo_tax-consist_data-2}
# Read in data
withzeros <- read.csv(file.path("data", "simulated-taxa-df_with-zeros.csv")) %>% 
  dplyr::mutate(plot = factor(plot))

# Check structure
str(withzeros) # <1>
```
1. Notice how there are more rows than the preceding data object and several new zeroes in the first few rows?

Now that we have the data in the format we need, we'll create a plot of species counts over time with zeros filled in. Because there are many plots and it is difficult to see so many panels on the same page, we'll use the `facet_wrap_paginate` function from the `ggforce` package to create a multi-page PDF output.

```{r demo_tax-consist_viz-code-fake}
#| eval: false

# Create the plot of species counts over time (with zeros filled in)
myplot <- ggplot(withzeros, aes(x = year, y = n, group = plot, color = plot)) +
  geom_line() +
  scale_x_continuous(breaks = scales::pretty_breaks()) +
  ggforce::facet_wrap_paginate(~ taxon, nrow = 2, ncol = 2) # <1>

# Start the PDF output
grDevices::pdf(file.path("counts_by_taxon_with_zeros.pdf"),
               width = 9, height = 5)

# Loop across pages (defined by `ggforce::facet_wrap_paginate`)
for (i in seq_along(ggforce::n_pages(myplot))) {
  
  page_plot <- myplot + 
      ggforce::facet_wrap_paginate(~taxon, page = i, 
                                   nrow = 2, ncol = 2)
  
  print(page_plot)
}

# Close the PDF output
dev.off()
```
1. This allows a faceted graph to spread across more than one page. See `?ggforce::facet_wrap_paginate` for details

The first page of the resulting plot should look something like this:

```{r demo_tax-consist_viz-code-real}
#| fig-align: center
#| fig-width: 7
#| fig-height: 5
#| message: false
#| echo: false

# Make multi-page graph
myplot <- ggplot(withzeros, aes(x = year, y = n, group = plot, color = plot)) +
  geom_line() +
  scale_x_continuous(breaks = scales::pretty_breaks()) +
  ggforce::facet_wrap_paginate(~ taxon, nrow = 2, ncol = 2)

# Check out just the first page
myplot +
  ggforce::facet_wrap_paginate(~ taxon, page = 1,
                               nrow = 2, ncol = 2)
```

Notice how "Taxon_A" is absent from all plots in 2014 whereas "Taxon_B" has extremely high counts in the same year. Often this can signify inconsistent use of taxonomic names over time.

### Seasonal Time Series

For time series, intra-annual variation can often make data issues difficult to spot. In these cases, it can be helpful to plot each year onto the same figure and compare trends across study years.

As an example, we'll use a 2024 dataset on streamflow near the [Niwot Ridge](https://nwt.lternet.edu/) LTER. The dataset is a 22 year time-series of daily streamflow. For more information on these data, check out [the data package on EDI](https://portal.edirepository.org/nis/mapbrowse?packageid=knb-lter-nwt.105.18). 

```{r demo_seasons_data}
# Read data
green_streamflow <- read.csv(file.path("data", "green-lakes_streamflow.csv")) # <1>

# Check structure
str(green_streamflow)
```
1. Note again that you could also read in this data directly from EDI. See ~line 129 of [this script](https://github.com/lter/ssecr/blob/main/scripts/prep-data_data-viz-demo.R) for a syntax example

Let's now calculate a moving average encompassing the 5 values before and after each focal value.

```{r demo_seasons_wrangle}
#| message: false

# Do necessary wrangling
stream_data <- green_streamflow %>%
  # Calculate moving average for each numeric variable
  dplyr::mutate(dplyr::across(.cols = dplyr::all_of(c("discharge", "temperature")),
                              .fns = ~ slider::slide_dbl(.x = .x, .f = mean,
                                                         .before = 5, .after = 5),
                              .names = "{.col}_move.avg" )) %>%
  # Handle date format issues
  dplyr::mutate(yday = lubridate::yday(date),
                year = lubridate::year(date))

# Check the structure of that
str(stream_data)
```

Plot seasonal timeseries of each numeric variable as points with the moving
average included as lines

```{r demo_seasons_viz-code-fake}
#| eval: false

# Start PDF output
grDevices::pdf(file = file.path("qc_all_numeric_seasonal.pdf"))

# Loop across variables
for (var in c("discharge", "temperature")) {
  
  # Make the graph
  myplot <- ggplot(stream_data, aes(x = yday, group = year, color = year)) +
    geom_point(aes(y = .data[[var]])) + # <1>
    geom_line(aes(y = .data[[paste0(var, "_move.avg")]])) + # <2>
    viridis::scale_color_viridis()
  
  # Print it
  print(myplot)
}

# End PDF creation
dev.off()
```
1. Add points based on the year
2. Adding lines based on the average

The first page of the resulting figure should look something like this:

```{r demo_seasons_viz-code-real}
#| fig-align: center
#| fig-width: 8
#| fig-height: 8
#| echo: false
#| warning: false

# Print the first graph
ggplot(stream_data, aes(x = yday, group = year, color = year)) +
  geom_point(aes(y = .data[["discharge"]])) +
  geom_line(aes(y = .data[["discharge_move.avg"]])) +
  viridis::scale_color_viridis()
```

One of these years is not like the others...
:::

## Multivariate Visualization

If you are working with multivariate data (i.e., data where multiple columns are all response variables collectively) you may need to use visualization methods unique to that data structure. For more information, check out the [bonus multivariate visualization module](https://lter.github.io/ssecr/mod_multivar-viz.html).

## Maps

You may find it valuable to create a map as an additional way of visualizing data. Many synthesis groups do this--particularly when there is a strong spatial component to the research questions and/or hypotheses. Check out the [bonus spatial data module](https://lter.github.io/ssecr/mod_spatial.html) for more information on map-making if this is of interest!

## Additional Resources

### Papers & Documents

- Chang, W. _et al._, [`ggplot2`: Elegant Graphics for Data Analysis](https://ggplot2-book.org/). 3^rd^ edition. **2023**.
- National Center for Ecological Analysis and Synthesis (NCEAS). [Colorblind Safe Color Schemes](https://www.nceas.ucsb.edu/sites/default/files/2022-06/Colorblind%20Safe%20Color%20Schemes.pdf). **2022**.
- Wilke, C.O. [Fundamentals of Data Visualization](https://clauswilke.com/dataviz/). **2020**.

### Workshops & Courses

- The Carpentries. [Data Analysis and Visualization in R for Ecologists: Data Visualization with `ggplot2`](https://datacarpentry.org/R-ecology-lesson/visualizing-ggplot.html). **2024**.
- The Carpentries. [Data Analysis and Visualization in Python for Ecologists: Making Plots with `plotnine`](https://datacarpentry.org/python-ecology-lesson/07-visualization-ggplot-python.html). **2024**.
- LTER Scientific Computing Team. [Coding in the Tidyverse: 'Visualize' Module](https://lter.github.io/workshop-tidyverse/visualize.html). **2023**.

### Websites

- [The R Graph Gallery](https://r-graph-gallery.com/)