Motivating Example: You have data from two groups and want to take a look before doing any formal statistics.
Learning Goals: By the end of this section, you will be able to:
For this section we will compare the average number of pollinator visits over a 15 minute interval to white and pink Clarkia RILs at site SR. I’ll load and prep the data to get started.
library(readr)
library(ggplot2)
library(dplyr)
ril_link <- "https://raw.githubusercontent.com/ybrandvain/datasets/refs/heads/master/clarkia_rils.csv"
SR_visits <- readr::read_csv(ril_link) |>
dplyr::rename(visits = mean_visits) |>
dplyr::filter(location == "SR",
!is.na(petal_color),
!is.na(visits)) |>
dplyr::select(petal_color, visits ) The first step in any data analysis is to look at the data. Here we’ll reintroduce some of the best ways to do this when comparing two groups. See our introduction to ggplot for our initial presentation of these ideas.
Parametric statistical approaches - like the two sample t-test, which we focus on in this chapter, make assumptions about the distribution of data. Specifically, linear models assume a normal distribution of residuals. We will consider these assumptions later in this chapter. But here I introduce two kinds of plots – histograms and density plots – that highlight the distribution of the data.
I have introduced the idea of visualizing a continuous response variable as a function of a categorical explanatory variable previously. See the sections on visualizing one continuous variable or a continuous response to a categorical predictor in our introduction to ggplot.
Histograms show the distribution of our data by putting each data point in a “bin” within a range. That bin-range is shown on the x-axis. The y-axis – labelled “count” by default in ggplot – shows the number of data points in each bin. Of course, there is rarely a “natural” or “correct” size of a bin, rather we must choose which size is appropriate. This choice involves a trade-off:
If you do not make an explicit choice, you are choosing 30 bins (R’s default).
Fortunately, you are not bound by an initial choice. In fact, I suggest experimenting with a few bin sizes to help you best understand your data, before choosing one that best communicates your results. There are two ways to specify the smoothing of a histogram in ggplot:
binwidth argument in geom_histogram().bins argument in geom_histogram().There are a few ways to use a histogram to compare the distributions between groups. The code below shows that you can do this by noting group by color and piling observations on top of one another. But, I find this difficult to interpret, and prefer using small multiples by adding:
facet_wrap(~petal_color, ncol = 1) to show the data in one column (“commented out” below).facet_wrap(~petal_color, nrow = 1) to show the data in one row (change ncol = 1 to nrow = 1 in the “commented out” section below).Use the webR environment, below, to compare the distribution of the pollinator visits of white- and pink-flowered Clarkia RILs. Then answer the challenge questions. In doing so:
ncol and nrow arguments in the facet_wrap function.
Before starting these questions, change the number or width of bins to find a reasonable “smoothing”.
Q1. Which of the following is FALSE?Q2. Because on average pink-flowered genotypes receive more pollinator visits than white-flowered genotypes, the data appear bimodal before considering petal color.
Q3. Visually approximate the binwidth for the data shown in Figure 1. .
Rather than counting the observations in a bin, we can fit a mathematical function to the shape of the distribution. Such “density plots” often allow readers to better compare distributions than histograms. A density plot displays the same message as a histogram, but now the y-axis is labelled “density” and represents the relative proportion of the data estimated to lie at a specific x-value from the model of counts. Thus the area under the curve integrates to one.
As in a histogram, we can control the smoothness of a density plot, and again we have the same trade-off between smoothing over noise (raising geom_density’s adjust argument above one), and revealing subtle patterns (lowering geom_density’s adjust argument below one).
Finding the right transparency.
A major advantage of density plots over histograms is that we can visually compare overlaid distributions. However to do so we must be able to see both distributions. So we must make the plots semi-transparent. In ggplot, this is controlled by the alpha argument.
Use the webR environment, below, to compare the distribution of the pollinator visits of white- and pink-flowered Clarkia RILs. Then answer the challenge questions. In doing so:
alpha to 0.8, 0.5, 0.2, and 0.adjust to a stupidly large number. Start by trying adjust values of 2, 1, and 0.5.
Use your R efforts and Figure 2 below to answer the following questions:
Q4. geom_density(adjust =0.3, alpha= 1) makes .
Q5.geom_density(adjust =1, alpha= 1) makes .
Q6. geom_density(adjust= 3, alpha= 0.1) makes .
Density plots do not show the actual data. Rather density plots show a function that is meant to approximate the data. Always look at the actual data to be sure that the density plot is not misleading you.
Histograms and density plots are the first step in exploratory data analysis. As such, these plots are for the person analyzing the data (you), and others deeply interested in the data analysis (your collaborators and mentors etc). Histograms and density plots are rarely best for communicating differences in means for a broad audience. Rather than highlighting the distribution, a good presentation of a between-group comparison should show:
Let’s get started with our first attempt at such a plot. To show our data, we map our explanatory variable onto the x-axis, and our response variable onto the y-axis. To show the estimated mean and uncertainty about it, we use the stat_summary() function. It is best to have our plots match our stats, so in stat_summary() we set fun.data = "mean_cl_normal". When we specify mean_cl_normal for the fun.data argument, ggplot uses a function from the Hmisc package, so you may need to install it by entering install.packages("Hmisc") and load it by entering library(Hmisc).
library(Hmisc)
SR_visits |>
ggplot(aes(x = petal_color, y = visits)) +
geom_point()+
stat_summary(fun.data = "mean_cl_normal")
Say what error bars represent! Most scientists know to add “error bars” to their plots. But unfortunately, many different “error bars” accompany an estimate. The most common are 95% confidence intervals (like those above), or some number of standard errors. Sometimes people use bars, not to show uncertainty but rather to summarize variability (e.g. as standard deviations).
The best practice is to show 95% confidence intervals and tell the reader that in the figure caption.
Showing all the data plus the estimated mean and uncertainty is not as easy as it sounds. As you can see from Figure 3, simply using geom_point() hides a bunch of data. For example, we know from our histograms that white-flowered genotypes have many “zero” visits, but we can only see one because they are all on top of each other. This phenomenon, known as “overplotting”, ironically hides your data in trying to show it all.
The solution to overplotting is to spread the data out a bit on the x-axis. I often accompany this with slightly larger and more transparent points than I would use otherwise. In the introduction to ggplot, we saw how to use geom_jitter() to overcome overplotting (for example geom_jitter(size = 3, alpha = 0.7, height = 0, width = 0.1)). Here, I introduce the sina plot from the ggforce package as an alternative to a traditional jitter.
(almost) always set height = 0 in geom_jitter() – We do not want to mislead by giving the impression of excess variation in our response variable. Similarly, you can experiment with width but be sure that a reader can unambiguously know which data point is associated with which value of the categorical predictor.
As you can see from Figure 4, the sina plot spreads data out - much like a jitter plot - but in doing so, it makes the local spread of the data on the x-axis proportional to the density of these data. I also made the points larger and more transparent to better see each data point, and showed estimated means and uncertainty in red to ensure these summaries do not get confused for raw data.
library(Hmisc)
library(ggforce)
SR_visits |>
ggplot(aes(x = petal_color, y = visits)) +
geom_sina(size = 3,alpha = .5)+
stat_summary(fun.data = "mean_cl_normal", color = "red",
size = .7, linewidth = 1)
From our visualizations we see that very few white- and pink-flowered RILs receive many visits while most receive few visits - i.e. the data are right-skewed. In particular many white-flowered RILs received no visitors. We also see that pink-flowered RILs receive more pollinator visits on average than do white-flowered RILs, and that the 95% confidence intervals hardly overlap.
The rest of this chapter is dedicated to estimating the difference in means, the uncertainty about this difference and how to test the null that the observed difference represents sampling error from a single population with the same mean.
To do this, we will build off of our understanding of the t-distribution from the previous chapter. But doing this requires that we evaluate whether data fit the model assumptions.