<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/exploring-ggplot2-boxplots-defining-limits-and-adjusting-style/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://owi.usgs.gov/blog/boxplots/"> The USGS OWI blog </a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>Boxplots are often used to show data distributions, and <code>ggplot2</code> is often used to visualize data. A question that comes up is what exactly do the box plots represent? The <code>ggplot2</code> box plots follow standard Tukey representations, and there are many references of this online and in standard statistical text books. The base R function to calculate the box plot limits is <code>boxplot.stats</code>. The help file for this function is very informative, but it’s often non-R users asking what exactly the plot means. Therefore, this blog post breaks down the calculations into (hopefully!) easy-to-follow chunks of code for you to make your own box plot legend if necessary. Some additional goals here are to create boxplots that come <em>close</em> to USGS style. Features in this blog post take advantage of enhancements to <code>ggplot2</code> in version 3.0.0 or later.</p>
<p>First, let’s get some data that might be typically plotted in a USGS report using a boxplot. Here we’ll use chloride data (parameter code “00940”) measured at a USGS station on the Fox River in Green Bay, WI (station ID “04085139”). We’ll use the package <code>dataRetrieval</code> to get the data (see <a rel="nofollow" href="https://owi.usgs.gov/R/dataRetrieval.html" target="_blank">this tutorial</a> for more information on <code>dataRetrieval</code>), and plot a simple boxplot by month using <code>ggplot2</code>:</p>
<pre><code>library(dataRetrieval)
library(ggplot2)

# Get chloride data using dataRetrieval:
chloride &lt;- readNWISqw(&quot;04085139&quot;, &quot;00940&quot;)
# Add a month column:
chloride$month &lt;- month.abb[as.numeric(format(chloride$sample_dt, &quot;%m&quot;))]
chloride$month &lt;- factor(chloride$month, labels = month.abb)
# Pull out the official parameter and site names for labels:
cl_name &lt;- attr(chloride, &quot;variableInfo&quot;)[[&quot;parameter_nm&quot;]]
cl_site &lt;- attr(chloride, &quot;siteInfo&quot;)[[&quot;station_nm&quot;]]
# Create basic ggplot graph:
ggplot(data = chloride, 
       aes(x = month, y = result_va)) +
  geom_boxplot() +
  xlab(&quot;Month&quot;) +
  ylab(cl_name) +
  labs(title = cl_site)
</code></pre>
<p><img alt="ggplot2 defaults for boxplots." title="ggplot2 defaults for boxplots." src="https://i1.wp.com/owi.usgs.gov/blog/static/boxplots/getChoride-1.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>Is that graph great? YES! And for presentations and/or journal publications, that graph might be appropriate. However, for an official USGS report, USGS employees need to get the graphics approved to assure they follow specific style guidelines. The approving officer would probably come back from the review with the following comments:</p>
<table>
<colgroup>
<col width="51%">
<col width="48%">
</colgroup>
<thead>
<tr>
<th>Reviewer’s Comments</th>
<th><code>ggplot2</code> element to modify</th>
</tr>
</thead>
<tbody>
<tr>
<td>Remove background color, grid lines</td>
<td>Adjust theme</td>
</tr>
<tr>
<td>Add horizontal bars to the upper and lower whiskers</td>
<td>Add <code>stat_boxplot</code></td>
</tr>
<tr>
<td>Have tick marks go inside the plot</td>
<td>Adjust theme</td>
</tr>
<tr>
<td>Tick marks should be on both sides of the y axis</td>
<td>Add <code>sec.axis</code> to <code>scale_y_continuous</code></td>
</tr>
<tr>
<td>Remove tick marks from discrete data</td>
<td>Adjust theme</td>
</tr>
<tr>
<td>y-axis needs to start exactly at 0</td>
<td>Add <code>expand_limits</code></td>
</tr>
<tr>
<td>y-axis labels need to be shown at 0 and at the upper scale</td>
<td>Add <code>breaks</code> and <code>limits</code> to <code>scale_y_continuous</code></td>
</tr>
<tr>
<td>Add very specific legend</td>
<td>Create function <code>ggplot_box_legend</code></td>
</tr>
<tr>
<td>Add the number of observations above each boxplot</td>
<td>Add custom <code>stat_summary</code></td>
</tr>
<tr>
<td>Change text size</td>
<td>Adjust <code>geom_text</code> defaults</td>
</tr>
<tr>
<td>Change font (we’ll use &quot;serif&quot; in this blog, although that is not the official USGS font)</td>
<td>Adjust <code>geom_text</code> defaults</td>
</tr>
</tbody>
</table>
<p>As you can see, it will not be as simple as creating a single custom ggplot theme to comply with the requirements. However, we can string together ggplot commands in a list for easy re-use. This blog is <em>not</em> going to get you perfect compliance with the USGS standards, but it will get much closer. Also, while these style adjustments are tailored to USGS requirements, the process described here may be useful for other graphic guidelines as well.</p>
<p>So, let’s skip to the exciting conclusion and use some code that will be described later (<code>boxplot_framework</code> and <code>ggplot_box_legend</code>) to create the same plot, now closer to those USGS style requirements:</p>
<pre><code>library(cowplot)

# NOTE! This is a preview of the FUTURE!
# We'll create the functions ggplot_box_legend and boxplot_framework
# later in this blog post. 
# So....by the end of this post, you will be able to:
legend_plot &lt;- ggplot_box_legend()

chloride_plot &lt;- ggplot(data = chloride, 
       aes(x = month, y = result_va)) +
  boxplot_framework(upper_limit = 70) + 
  xlab(label = &quot;Month&quot;) +
  ylab(label = cl_name) +
  labs(title = cl_site)

plot_grid(chloride_plot, 
          legend_plot,
          nrow = 1, rel_widths = c(.6,.4))
</code></pre>
<p><img alt="Chloride by month, styled." title="Chloride by month styled." src="https://i1.wp.com/owi.usgs.gov/blog/static/boxplots/chlorideWithLegend-1.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>As can be seen in the code chunk, we are now using a function <code>ggplot_box_legend</code> to make a legend, <code>boxplot_framework</code> to accommodate all of the style requirements, and the <code>cowplot</code> package to plot them together.</p>
<h1><code>boxplot_framework</code>: Styling the boxplot</h1>
<p>Let’s get our style requirements figured out. First, we can set some basic plot elements for a theme. We can start with the <code>theme_bw</code> and add to that. Here we remove the grid, set the size of the title, bring the y-ticks inside the plotting area, and remove the x-ticks:</p>
<pre><code>theme_USGS_box &lt;- function(base_family = &quot;serif&quot;, ...){
  theme_bw(base_family = base_family, ...) +
  theme(
    panel.grid = element_blank(),
    plot.title = element_text(size = 8),
    axis.ticks.length = unit(-0.05, &quot;in&quot;),
    axis.text.y = element_text(margin=unit(c(0.3,0.3,0.3,0.3), &quot;cm&quot;)), 
    axis.text.x = element_text(margin=unit(c(0.3,0.3,0.3,0.3), &quot;cm&quot;)),
    axis.ticks.x = element_blank(),
    aspect.ratio = 1,
    legend.background = element_rect(color = &quot;black&quot;, fill = &quot;white&quot;)
  )
}
</code></pre>
<p>Next, we can change the defaults of the geom_text to a smaller size and font.</p>
<pre><code>update_geom_defaults(&quot;text&quot;, 
                   list(size = 3, 
                        family = &quot;serif&quot;))
</code></pre>
<p>We also need to figure out what other <code>ggplot2</code> functions need to be added. The basic ggplot code for the chloride plot would be:</p>
<pre><code>n_fun &lt;- function(x){
  return(data.frame(y = 0.95*70,
                    label = length(x)))
}

ggplot(data = chloride, 
       aes(x = month, y = result_va)) +
  stat_boxplot(geom ='errorbar', width = 0.6) +
  geom_boxplot(width = 0.6, fill = &quot;lightgrey&quot;) +
  stat_summary(fun.data = n_fun, geom = &quot;text&quot;, hjust = 0.5) +
  expand_limits(y = 0) +
  theme_USGS_box() +
  scale_y_continuous(sec.axis = dup_axis(label = NULL, 
                                         name = NULL),
                     expand = expand_scale(mult = c(0, 0)),
                     breaks = pretty(c(0,70), n = 5), 
                     limits = c(0,70))
</code></pre>
<p>Breaking that code down:</p>
<table>
<colgroup>
<col width="41%">
<col width="58%">
</colgroup>
<thead>
<tr>
<th>Function</th>
<th>What’s happening?</th>
</tr>
</thead>
<tbody>
<tr>
<td>stat_boxplot(geom =’errorbar’)</td>
<td>The &quot;errorbars&quot; are used to make the horizontal lines on the upper and lower whiskers. This needs to happen first so it is in the back of the plot.</td>
</tr>
<tr>
<td>geom_boxplot</td>
<td>Regular boxplot</td>
</tr>
<tr>
<td>stat_summary(fun.data = n_fun, geom = &quot;text&quot;, hjust = 0.5)</td>
<td>The <code>stat_summary</code> function is very powerful for adding specific summary statistics to the plot. In this case, we are adding a <code>geom_text</code> that is calculated with our custom <code>n_fun</code>. That function comes back with the count of the boxplot, and puts it at 95% of the hard-coded upper limit.</td>
</tr>
<tr>
<td>expand_limits</td>
<td>This forces the plot to include 0.</td>
</tr>
<tr>
<td>theme_USGS_box</td>
<td>Theme created above to help with grid lines, tick marks, axis size/fonts, etc.</td>
</tr>
<tr>
<td>scale_y_continuous</td>
<td>A tricky part of the USGS requirements involve 4 parts: Add ticks to the right side, have at least 4 &quot;pretty&quot; labels on the left axis, remove padding, and have the labels start and end at the beginning and end of the plot. Breaking that down further:</td>
</tr>
<tr>
<td>scale_y_continuous(sec.axis = dup_axis</td>
<td>Handy function to add tick marks to the right side of the graph.</td>
</tr>
<tr>
<td>scale_y_continuous(expand = expand_scale(mult = c(0, 0))</td>
<td>Remove padding</td>
</tr>
<tr>
<td>scale_y_continuous(breaks = pretty(c(0,70), n = 5))</td>
<td>Make pretty label breaks, assuring 5 pretty labels if the graph went from 0 to 70</td>
</tr>
<tr>
<td>scale_y_continuous(limits = c(0,70))</td>
<td>Assure the graph goes from 0 to 70.</td>
</tr>
</tbody>
</table>
<p>Let’s look at a few other common boxplots to see if there are other ggplot2 elements that would be useful in a common <code>boxplot_framework</code> function.</p>
<h2>Logrithmic boxplot</h2>
<p>For another example, we might need to make a boxplot with a logarithm scale. This data is for phosphorus measurements on the Pheasant Branch Creek in Middleton, WI.</p>
<pre><code>site &lt;- &quot;05427948&quot;
pCode &lt;- &quot;00665&quot;

# Get phosphorus data using dataRetrieval:
phos_data &lt;- readNWISqw(site, pCode)
# Create a month column:
phos_data$month &lt;- month.abb[as.numeric(format(phos_data$sample_dt, &quot;%m&quot;))]
phos_data$month &lt;- factor(phos_data$month, labels = month.abb)
# Get site name and paramter name for labels:
phos_name &lt;- attr(phos_data, &quot;variableInfo&quot;)[[&quot;parameter_nm&quot;]]
phos_site &lt;- attr(phos_data, &quot;siteInfo&quot;)[[&quot;station_nm&quot;]]

n_fun &lt;- function(x){
    return(data.frame(y = 0.95*log10(50),
                      label = length(x)))
}

prettyLogs &lt;- function(x){
  pretty_range &lt;- range(x[x > 0])
  pretty_logs &lt;- 10^(-10:10)
  log_index &lt;- which(pretty_logs &lt; pretty_range[2] &amp; 
                       pretty_logs > pretty_range[1])
  log_index &lt;- c(log_index[1]-1,log_index, log_index[length(log_index)]+1)
  pretty_logs_new &lt;-  pretty_logs[log_index] 
  return(pretty_logs_new)
}

fancyNumbers &lt;- function(n){
  nNoNA &lt;- n[!is.na(n)]
  x &lt;-gsub(pattern = &quot;1e&quot;,replacement = &quot;10^&quot;,
           x = format(nNoNA, scientific = TRUE))
  exponents &lt;- as.numeric(sapply(strsplit(x, &quot;\\^&quot;), function(j) j[2]))

  base &lt;- ifelse(exponents == 0, &quot;1&quot;, ifelse(exponents == 1, &quot;10&quot;,&quot;10^&quot;))
  exponents[base == &quot;1&quot; | base == &quot;10&quot;] &lt;- &quot;&quot;
  textNums &lt;- rep(NA, length(n))  
  textNums[!is.na(n)] &lt;- paste0(base,exponents)
  
  textReturn &lt;- parse(text=textNums)
  return(textReturn)
}


phos_plot &lt;- ggplot(data = phos_data, 
       aes(x = month, y = result_va)) +
  stat_boxplot(geom ='errorbar', width = 0.6) +
  geom_boxplot(width = 0.6, fill = &quot;lightgrey&quot;) +
  stat_summary(fun.data = n_fun, geom = &quot;text&quot;, hjust = 0.5) +
  theme_USGS_box() +
  scale_y_log10(limits = c(0.001, 50),
                expand = expand_scale(mult = c(0, 0)),
                labels=fancyNumbers,
                breaks=prettyLogs) +
  annotation_logticks(sides = c(&quot;rl&quot;)) +
  xlab(&quot;Month&quot;) +
  ylab(phos_name) +
  labs(title = phos_site) 

phos_plot
</code></pre>
<p><img alt="Phosphorus distribution by month." title="Phosphorus distribution by month." src="https://i0.wp.com/owi.usgs.gov/blog/static/boxplots/phosDistribution-1.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>What are the new features we have to consider for log scales?</p>
<table>
<colgroup>
<col width="41%">
<col width="58%">
</colgroup>
<thead>
<tr>
<th>Function</th>
<th>What’s happening?</th>
</tr>
</thead>
<tbody>
<tr>
<td>stat_boxplot</td>
<td>The <code>stat_boxplot</code> function is the same, but our custom function to calculate counts need to be adjusted so the position would be in log units.</td>
</tr>
<tr>
<td>scale_y_log10</td>
<td>This is used instead of <code>scale_y_continuous</code>.</td>
</tr>
<tr>
<td>annotation_logticks(sides = c(&quot;rl&quot;))</td>
<td>Adds nice log ticks to the right (&quot;r&quot;) and left (&quot;l&quot;) side.</td>
</tr>
<tr>
<td>prettyLogs</td>
<td>This function forces the y-axis breaks to be on every 10^x. This could be adjusted if a finer scale was needed.</td>
</tr>
<tr>
<td>fancyNumbers</td>
<td>This is a custom formatting function for the log axis. This function could be adjusted if other formatting was needed.</td>
</tr>
</tbody>
</table>
<h2>Grouped boxplots</h2>
<p>We might also want to make grouped boxplots. In ggplot, it’s pretty easy to add a “fill” to the <code>aes</code> argument. Here we’ll plot temperature distributions at 4 USGS stations. We’ll group the measurements by a “daytime” and “nighttime” factor. Temperature might be a parameter that would not be required to start at 0.</p>
<pre><code>library(dplyr)

# Get water temperature data for a variety of USGS stations
temp_q_data &lt;- readNWISuv(siteNumbers = c(&quot;04026561&quot;, &quot;04063700&quot;,
                                          &quot;04082400&quot;, &quot;05427927&quot;),
                          parameterCd = '00010', 
                          startDate = &quot;2018-06-01&quot;, 
                          endDate = &quot;2018-06-03&quot;)

temperature_name &lt;- attr(temp_q_data, &quot;variableInfo&quot;)[[&quot;variableName&quot;]]

# add an hour of day to create groups (daytime or nighttime)
temp_q_data &lt;- temp_q_data %>% 
  renameNWISColumns() %>%
  mutate(hourOfDay = as.numeric(format(dateTime, &quot;%H&quot;)),
         timeOfDay = case_when(hourOfDay &lt; 20 &amp; hourOfDay > 6 ~ &quot;daytime&quot;,
                               TRUE ~ &quot;nighttime&quot; # catchall 
  ))

n_fun &lt;- function(x){
  return(data.frame(y = 0.95*30,
                    label = length(x)))
}

temperature_plot &lt;- ggplot(data = temp_q_data, 
       aes(x=site_no, y=Wtemp_Inst, fill=timeOfDay)) +
  stat_boxplot(geom ='errorbar', width = 0.6) +
  geom_boxplot(width = 0.6) +
  stat_summary(fun.data = n_fun, geom = &quot;text&quot;, 
               aes(group=timeOfDay),
               hjust = 0.5, position = position_dodge(0.6)) +
  expand_limits(y = 0) +
  scale_y_continuous(sec.axis = dup_axis(label = NULL, 
                                         name = NULL),
                     expand = expand_scale(mult = c(0, 0)),
                     breaks = pretty(c(10,30), n = 5), 
                     limits = c(10,30)) +
  theme_USGS_box() +
  xlab(&quot;Station ID&quot;) +
  ylab(temperature_name) +
  scale_fill_discrete(name = &quot;EXPLANATION&quot;) +
  theme(legend.position = c(0.175, 0.78))

temperature_plot
</code></pre>
<p><img alt="Grouped boxplots" title="Grouped boxplots." src="https://i1.wp.com/owi.usgs.gov/blog/static/boxplots/groupPlot1-1.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>What are the new features we have to consider for log scales?</p>
<table>
<colgroup>
<col width="41%">
<col width="58%">
</colgroup>
<thead>
<tr>
<th>Function</th>
<th>What’s happening?</th>
</tr>
</thead>
<tbody>
<tr>
<td>stat_summary(position)</td>
<td>We need to move the counts to above the boxplots. This is done by shifting them the same amount as the width.</td>
</tr>
<tr>
<td>stat_summary(aes(group=timeOfDay))</td>
<td>We need to include how the boxplots are grouped.</td>
</tr>
<tr>
<td>scale_fill_discrete</td>
<td>Need include a fill legend.</td>
</tr>
</tbody>
</table>
<p>Additionally, the parameter name that comes back from <code>dataRetrieval</code> could use some formatting. The following function can fix that for both <code>ggplot2</code> and base R graphics:</p>
<pre><code>unescape_html &lt;- function(str){
  fancy_chars &lt;- regmatches(str, gregexpr(&quot;&amp;#\\d{3};&quot;,str)) 

  unescaped &lt;- xml2::xml_text(xml2::read_html(paste0(&quot;&quot;)))

  fancy_chars &lt;- gsub(pattern = &quot;&amp;#\\d{3};&quot;, 
                      replacement = unescaped, x = str)

  fancy_chars &lt;- gsub(&quot;Â&quot;,&quot;&quot;, fancy_chars)
  return(fancy_chars)
}
</code></pre>
<p>We’ll use this function in the next section.</p>
<h2>Framework function</h2>
<p>Finally, we can bring all of those elements together into a single list for <code>ggplot2</code> to use. While we’re at it, we can create a function that is flexible for both linear and logarithmic scales, as well as grouped boxplots. It’s a bit clunky because you need to specify the upper and lower limits of the plot. Much of the USGS style requirements depend on specific upper and lower limits, so I decided this was an acceptable solution for this blog post. There’s almost certainly a slicker way to do that, but for now, it works:</p>
<pre><code>boxplot_framework &lt;- function(upper_limit, 
                              family_font = &quot;serif&quot;,
                              lower_limit = 0, 
                              logY = FALSE, 
                              fill_var = NA,
                              fill = &quot;lightgrey&quot;, width = 0.6){
  
  update_geom_defaults(&quot;text&quot;, 
                     list(size = 3, 
                          family = family_font))
  
  n_fun &lt;- function(x, lY = logY){
    return(data.frame(y = ifelse(logY, 0.95*log10(upper_limit), 0.95*upper_limit),
                      label = length(x)))
  }
  
  prettyLogs &lt;- function(x){
    pretty_range &lt;- range(x[x > 0])
    pretty_logs &lt;- 10^(-10:10)
    log_index &lt;- which(pretty_logs &lt; pretty_range[2] &amp; 
                         pretty_logs > pretty_range[1])
    log_index &lt;- c(log_index[1]-1,log_index,
                   log_index[length(log_index)]+1)
    pretty_logs_new &lt;-  pretty_logs[log_index] 
    return(pretty_logs_new)
  }
  
  fancyNumbers &lt;- function(n){
    nNoNA &lt;- n[!is.na(n)]
    x &lt;-gsub(pattern = &quot;1e&quot;,replacement = &quot;10^&quot;,
             x = format(nNoNA, scientific = TRUE))
    exponents &lt;- as.numeric(sapply(strsplit(x, &quot;\\^&quot;), function(j) j[2]))
  
    base &lt;- ifelse(exponents == 0, &quot;1&quot;, ifelse(exponents == 1, &quot;10&quot;,&quot;10^&quot;))
    exponents[base == &quot;1&quot; | base == &quot;10&quot;] &lt;- &quot;&quot;
    textNums &lt;- rep(NA, length(n))  
    textNums[!is.na(n)] &lt;- paste0(base,exponents)
    
    textReturn &lt;- parse(text=textNums)
    return(textReturn)
  }
  
  if(!is.na(fill_var)){
      basic_elements &lt;- list(stat_boxplot(geom ='errorbar', width = width),
                            geom_boxplot(width = width),
                            stat_summary(fun.data = n_fun, 
                                         geom = &quot;text&quot;, 
                                         position = position_dodge(width),
                                         hjust =0.5,
                                         aes_string(group=fill_var)),
                            expand_limits(y = lower_limit),
                            theme_USGS_box())
  } else {
      basic_elements &lt;- list(stat_boxplot(geom ='errorbar', width = width),
                            geom_boxplot(width = width, fill = fill),
                            stat_summary(fun.data = n_fun, 
                                         geom = &quot;text&quot;, hjust =0.5),
                            expand_limits(y = lower_limit),
                            theme_USGS_box())
  }
  
  if(logY){
    return(c(basic_elements,
              scale_y_log10(limits = c(lower_limit, upper_limit),
                  expand = expand_scale(mult = c(0, 0)),
                  labels=fancyNumbers,
                  breaks=prettyLogs),
              annotation_logticks(sides = c(&quot;rl&quot;))))      
  } else {
    return(c(basic_elements,
              scale_y_continuous(sec.axis = dup_axis(label = NULL, 
                                                     name = NULL),
                                 expand = expand_scale(mult = c(0, 0)),
                                 breaks = pretty(c(lower_limit,upper_limit), n = 5), 
                                 limits = c(lower_limit,upper_limit))))    
  }
}
</code></pre>
<h1>Examples with our framework</h1>
<p>Let’s see if it works! Let’s build the last set of example figures using our new function <code>boxplot_framework</code>. I’m also going to use the ‘cowplot’ package to print them all together. I’ll also include the <code>ggplot_box_legend</code> which will be described in the next section.</p>
<pre><code>legend_plot &lt;- ggplot_box_legend()

chloride_plot &lt;- ggplot(data = chloride, 
       aes(x = month, y = result_va)) +
  boxplot_framework(upper_limit = 70) + 
  xlab(label = &quot;Month&quot;) +
  ylab(label = cl_name) +
  labs(title = cl_site)

phos_plot &lt;- ggplot(data = phos_data, 
       aes(x = month, y = result_va)) +
  boxplot_framework(upper_limit = 50,
                    lower_limit = 0.001,
                    logY = TRUE) +
  xlab(label = &quot;Month&quot;) +
  #Shortened label since the graph area is smaller
  ylab(label = &quot;Phosphorus in milligraphs per liter&quot;) + 
  labs(title = phos_site) 

temperature_plot &lt;- ggplot(data = temp_q_data, 
       aes(x=site_no, y=Wtemp_Inst, fill=timeOfDay)) +
  boxplot_framework(upper_limit = 30, 
                    lower_limit = 10, 
                    fill_var = &quot;timeOfDay&quot;) + 
  xlab(label = &quot;Station ID&quot;) +
  ylab(label = unescape_html(temperature_name)) +
  labs(title = &quot;Daytime vs Nighttime Temperature Distribution&quot;) +
  scale_fill_discrete(name = &quot;EXPLANATION&quot;) +
  theme(legend.position = c(0.225, 0.72),
        legend.title =  element_text(size = 7))

plot_grid(chloride_plot, 
          phos_plot,
          temperature_plot,
          legend_plot,
          nrow = 2)
</code></pre>
<p><img src="https://i0.wp.com/owi.usgs.gov/blog/static/boxplots/comboGraph-1.png?w=456&ssl=1" alt="Combining boxplots using framework function and cowplot" title="Combining boxplots using framework function and cowplot" data-recalc-dims="1"></p>
<h1><code>ggplot_box_legend</code>: What is a boxplot?</h1>
<p>A non-trivial requirement to the USGS boxplot style guidelines is to make a detailed, prescribed legend. In this section we’ll first verify that <code>ggplot2</code> boxplots use the same definitions for the lines and dots, and then we’ll make a function that creates the prescribed legend. To start, let’s set up random data using the R function <code>sample</code> and then create a function to calculate each value.</p>
<pre><code>set.seed(100)

sample_df &lt;- data.frame(parameter = &quot;test&quot;,
                        values = sample(500))

# Extend the top whisker a bit:
sample_df$values[1:100] &lt;- 701:800
# Make sure there's only 1 lower outlier:
sample_df$values[1] &lt;- -350
</code></pre>
<p>Next, we’ll create a function that calculates the necessary values for the boxplots:</p>
<pre><code>ggplot2_boxplot &lt;- function(x){
  
  quartiles &lt;- as.numeric(quantile(x, 
                                   probs = c(0.25, 0.5, 0.75)))
  
  names(quartiles) &lt;- c(&quot;25th percentile&quot;, 
                        &quot;50th percentile\n(median)&quot;,
                        &quot;75th percentile&quot;)
  
  IQR &lt;- diff(quartiles[c(1,3)])

  upper_whisker &lt;- max(x[x &lt; (quartiles[3] + 1.5 * IQR)])
  lower_whisker &lt;- min(x[x > (quartiles[1] - 1.5 * IQR)])
    
  upper_dots &lt;- x[x > (quartiles[3] + 1.5*IQR)]
  lower_dots &lt;- x[x &lt; (quartiles[1] - 1.5*IQR)]

  return(list(&quot;quartiles&quot; = quartiles,
              &quot;25th percentile&quot; = as.numeric(quartiles[1]),
              &quot;50th percentile\n(median)&quot; = as.numeric(quartiles[2]),
              &quot;75th percentile&quot; = as.numeric(quartiles[3]),
              &quot;IQR&quot; = IQR,
              &quot;upper_whisker&quot; = upper_whisker,
              &quot;lower_whisker&quot; = lower_whisker,
              &quot;upper_dots&quot; = upper_dots,
              &quot;lower_dots&quot; = lower_dots))
}

ggplot_output &lt;- ggplot2_boxplot(sample_df$values)
</code></pre>
<p>What are those calculations?</p>
<ul>
<li>Quartiles (25, 50, 75 percentiles), 50% is the median</li>
<li>Interquartile range is the difference between the 75th and 25th percentiles</li>
<li>The upper whisker is the maximum value of the data that is within 1.5 times the interquartile range over the 75th percentile.</li>
<li>The lower whisker is the minimum value of the data that is within 1.5 times the interquartile range under the 25th percentile.</li>
<li>Outlier values are considered any values over 1.5 times the interquartile range over the 75th percentile or any values under 1.5 times the interquartile range under the 25th percentile.</li>
</ul>
<p>Let’s check that the output matches <code>boxplot.stats</code>:</p>
<pre><code># Using base R:
base_R_output &lt;- boxplot.stats(sample_df$values)

# Some checks:

# Outliers:
all(c(ggplot_output[[&quot;upper_dots&quot;]], 
      ggplot_output[[&quot;lowerdots&quot;]]) %in%
    c(base_R_output[[&quot;out&quot;]]))
</code></pre>
<pre><code>## [1] TRUE
</code></pre>
<pre><code># whiskers:
ggplot_output[[&quot;upper_whisker&quot;]] == base_R_output[[&quot;stats&quot;]][5]
</code></pre>
<pre><code>## [1] TRUE
</code></pre>
<pre><code>ggplot_output[[&quot;lower_whisker&quot;]] == base_R_output[[&quot;stats&quot;]][1]
</code></pre>
<pre><code>## [1] TRUE
</code></pre>
<h2>Boxplot Legend</h2>
<p>Let’s use this information to generate a legend, and make the code reusable by creating a standalone function that we used in earlier code (<code>ggplot_box_legend</code>). There is a <em>lot</em> of <code>ggplot2</code> code to digest here. Most of it is style adjustments to approximate the USGS style guidelines for a boxplot legend.</p>
<p><p>
Show/Hide Code
</p>
<div>
<pre><code>ggplot_box_legend &lt;- function(family = &quot;serif&quot;){
  
  # Create data to use in the boxplot legend:
  set.seed(100)

  sample_df &lt;- data.frame(parameter = &quot;test&quot;,
                        values = sample(500))

  # Extend the top whisker a bit:
  sample_df$values[1:100] &lt;- 701:800
  # Make sure there's only 1 lower outlier:
  sample_df$values[1] &lt;- -350
  
  # Function to calculate important values:
  ggplot2_boxplot &lt;- function(x){
  
    quartiles &lt;- as.numeric(quantile(x, 
                                     probs = c(0.25, 0.5, 0.75)))
    
    names(quartiles) &lt;- c(&quot;25th percentile&quot;, 
                          &quot;50th percentile\n(median)&quot;,
                          &quot;75th percentile&quot;)
    
    IQR &lt;- diff(quartiles[c(1,3)])
  
    upper_whisker &lt;- max(x[x &lt; (quartiles[3] + 1.5 * IQR)])
    lower_whisker &lt;- min(x[x > (quartiles[1] - 1.5 * IQR)])
      
    upper_dots &lt;- x[x > (quartiles[3] + 1.5*IQR)]
    lower_dots &lt;- x[x &lt; (quartiles[1] - 1.5*IQR)]
  
    return(list(&quot;quartiles&quot; = quartiles,
                &quot;25th percentile&quot; = as.numeric(quartiles[1]),
                &quot;50th percentile\n(median)&quot; = as.numeric(quartiles[2]),
                &quot;75th percentile&quot; = as.numeric(quartiles[3]),
                &quot;IQR&quot; = IQR,
                &quot;upper_whisker&quot; = upper_whisker,
                &quot;lower_whisker&quot; = lower_whisker,
                &quot;upper_dots&quot; = upper_dots,
                &quot;lower_dots&quot; = lower_dots))
  }
  
  # Get those values:
  ggplot_output &lt;- ggplot2_boxplot(sample_df$values)
  
  # Lots of text in the legend, make it smaller and consistent font:
  update_geom_defaults(&quot;text&quot;, 
                     list(size = 3, 
                          hjust = 0,
                          family = family))
  # Labels don't inherit text:
  update_geom_defaults(&quot;label&quot;, 
                     list(size = 3, 
                          hjust = 0,
                          family = family))
  
  # Create the legend:
  # The main elements of the plot (the boxplot, error bars, and count)
  # are the easy part.
  # The text describing each of those takes a lot of fiddling to
  # get the location and style just right:
  explain_plot &lt;- ggplot() +
    stat_boxplot(data = sample_df,
                 aes(x = parameter, y=values),
                 geom ='errorbar', width = 0.3) +
    geom_boxplot(data = sample_df,
                 aes(x = parameter, y=values), 
                 width = 0.3, fill = &quot;lightgrey&quot;) +
    geom_text(aes(x = 1, y = 950, label = &quot;500&quot;), hjust = 0.5) +
    geom_text(aes(x = 1.17, y = 950,
                  label = &quot;Number of values&quot;),
              fontface = &quot;bold&quot;, vjust = 0.4) +
    theme_minimal(base_size = 5, base_family = family) +
    geom_segment(aes(x = 2.3, xend = 2.3, 
                     y = ggplot_output[[&quot;25th percentile&quot;]], 
                     yend = ggplot_output[[&quot;75th percentile&quot;]])) +
    geom_segment(aes(x = 1.2, xend = 2.3, 
                     y = ggplot_output[[&quot;25th percentile&quot;]], 
                     yend = ggplot_output[[&quot;25th percentile&quot;]])) +
    geom_segment(aes(x = 1.2, xend = 2.3, 
                     y = ggplot_output[[&quot;75th percentile&quot;]], 
                     yend = ggplot_output[[&quot;75th percentile&quot;]])) +
    geom_text(aes(x = 2.4, y = ggplot_output[[&quot;50th percentile\n(median)&quot;]]), 
              label = &quot;Interquartile\nrange&quot;, fontface = &quot;bold&quot;,
              vjust = 0.4) +
    geom_text(aes(x = c(1.17,1.17), 
                  y = c(ggplot_output[[&quot;upper_whisker&quot;]],
                        ggplot_output[[&quot;lower_whisker&quot;]]), 
                  label = c(&quot;Largest value within 1.5 times\ninterquartile range above\n75th percentile&quot;,
                            &quot;Smallest value within 1.5 times\ninterquartile range below\n25th percentile&quot;)),
                  fontface = &quot;bold&quot;, vjust = 0.9) +
    geom_text(aes(x = c(1.17), 
                  y =  ggplot_output[[&quot;lower_dots&quot;]], 
                  label = &quot;Outside value&quot;), 
              vjust = 0.5, fontface = &quot;bold&quot;) +
    geom_text(aes(x = c(1.9), 
                  y =  ggplot_output[[&quot;lower_dots&quot;]], 
                  label = &quot;-Value is >1.5 times and&quot;), 
              vjust = 0.5) +
    geom_text(aes(x = 1.17, 
                  y = ggplot_output[[&quot;lower_dots&quot;]], 
                  label = &quot;&lt;3 times the interquartile range\nbeyond either end of the box&quot;), 
              vjust = 1.5) +
    geom_label(aes(x = 1.17, y = ggplot_output[[&quot;quartiles&quot;]], 
                  label = names(ggplot_output[[&quot;quartiles&quot;]])),
              vjust = c(0.4,0.85,0.4), 
              fill = &quot;white&quot;, label.size = 0) +
    ylab(&quot;&quot;) + xlab(&quot;&quot;) +
    theme(axis.text = element_blank(),
          axis.ticks = element_blank(),
          panel.grid = element_blank(),
          aspect.ratio = 4/3,
          plot.title = element_text(hjust = 0.5, size = 10)) +
    coord_cartesian(xlim = c(1.4,3.1), ylim = c(-600, 900)) +
    labs(title = &quot;EXPLANATION&quot;)

  return(explain_plot) 
  
}

ggplot_box_legend()
</code></pre>
</div>
<p><img alt="ggplot2 box plot with explanation." title="ggplot2 box plot with explanation." src="https://i1.wp.com/owi.usgs.gov/blog/static/boxplots/visualizeBox-1.png?w=456&ssl=1" data-recalc-dims="1"></p>
<h1>Bring it together</h1>
<p>What’s nice about leaving this in the world of <code>ggplot2</code> is that it is still possible to use other <code>ggplot2</code> elements on the plot. For example, let’s add a reporting limit as horizontal lines to the phosphorous graph:</p>
<pre><code>phos_plot_with_DL &lt;- phos_plot +
  geom_hline(linetype = &quot;dashed&quot;,
             yintercept = 0.01)

explain_plot_DL &lt;- ggplot_box_legend() +
  geom_segment(aes(y = -650, yend = -650,
               x = 0.6, xend = 1.6),
           linetype=&quot;dashed&quot;) +
  geom_text(aes(y = -650, x = 1.8, label = &quot;Reporting Limit&quot;))

plot_grid(phos_plot_with_DL, 
          explain_plot_DL,
          nrow = 1, rel_widths = c(.6,.4))
</code></pre>
<p><img alt="Boxplot with additional ggplot2 features." title="Boxplot with additional ggplot2 features." src="https://i0.wp.com/owi.usgs.gov/blog/static/boxplots/picsWithDL-1.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>I hoped you like my “deep dive” into <code>ggplot2</code> boxplots. Many of the techniques here can be used to modify other <code>ggplot2</code> plots.</p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://owi.usgs.gov/blog/boxplots/"> The USGS OWI blog </a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/YV7dnd9VUPQ" height="1"></a></div>