Resampling-based methods for biologists

Understanding foundational concepts using simulation-based methods

Sampling distributions

The sampling distribution of a statistic tells us about the values we might expect for the statistic if we were to repeatedly collect data sets of the same size from the same population and using the same sampling protocols. Again, this is a difficult concept to grasp because we usually only get to collect one data set. However, it is easy to use a computer to generate multiple samples from a population, and hence calculate multiple statistics to approximate the sampling distribution. To illustrate, we explore here the sampling distribution of a sample mean using a data set containing the count of the number of polar bears within each of 164 roughly 3 × 3 km quadrats on Rowley Island in northern Foxe Basin, Nunavut (Fig. 1A; Stapleton et al., 2014). First, we read in the data set which contains two variables, an id for each quadrat (Quadrat) and a count of the number of bears in each quadrat (Num.Bears), making use of the here package to make it easy to read in data from a subdirectory (Müller, 2017). We use <- to assign the data to an object called bdat. We then use the head function to look at the first five rows of the data set (here and throughout the paper, we provide R code followed by its associated output).

bdat <- read.csv(here("data", "bears.csv"))

head(bdat, n=5)

Quadrat Num.Bears

1 1 2

2 2 1

3 3 0

4 4 0

5 5 1

We use the mean function to calculate the mean number of bears counted per quadrat in the population (i.e., the set of 164 quadrats). This function, like most in the mosaic package, follows a common syntax, which for statistical summaries of a single variable is: goal(∼variable name, data=Mydata) (Pruim, Kaplan & Horton, 2017).

mean(∼Num.Bears, data=bdat)

[1] 0.561

In this situation, we can calculate the population mean, $\mu =0.561$, exactly because we have polar bear counts for all sample units. Nonetheless, we will use these data to illustrate the concept of the sampling distribution. Specifically, we consider a hypothetical situation in which we have access to the counts of bears in each of 164 quadrats (i.e., the population), to look at what could have happened if we had only sampled a random subset of 75 quadrats. Our sample statistic, the sample mean ($\overline{x}$), provides an estimate of the population mean (μ) bear count per quadrat. Below, we repeatedly take simple random samples of 75 quadrats using the sample function and calculate the sample mean for each of the resulting data sets. We store the resulting sample means in an object named samp.mean. This feat can be accomplished with essentially one line of code using the do function in the mosaic package, where its argument (10,000 in this case) tells R how many times to repeatedly execute the code within the { } (To reproduce our results exactly, readers should set the seed of the random number generator using set.seed(03222007)):

samp.mean<-do(10000)*{

mean(∼Num.Bears, data=sample(bdat, size= 75))

}

The do function captures the 10,000 sample means and stores them in a variable named result within the samp.mean object. We can look at the first five sample means using the head function in R:

head(samp.mean, n=5)

result

1 0.507

2 0.600

3 0.453

4 0.480

5 0.520

We can easily visualize the sampling distribution using a histogram or density plot (Fig. 1C), and we can calculate the standard error using the sd function in the mosaic library (the standard error is just the standard deviation of the sampling distribution; Box 1):

(se.mean<-sd(∼result, data=samp.mean))

[1] 0.0772

When statistics have bell-shaped (i.e., approximately Normal) sampling distributions, we expect 95% of sample statistics to be within roughly 2 standard deviations of the mean of the sampling distribution (1.96 standard deviations to be more exact); we will refer to this result as the 2 SE rule (Box 1). We can use the above simulation results (i.e., the 10,000 sample means contained in samp.mean), to verify this claim. First, we find this interval by taking the mean of the sampling distribution ± 2 times the SE of the mean (we stored this SE in an object called se.mean and use the c operator to create a vector containing the values −2 and 2).

(rule.2SE<-mean(∼result, data=samp.mean)+c(−2, 2)*se.mean)

[1] 0.407 0.715

We then determine the proportion of sample means that fall in this interval (i.e., the proportion of sample means that are greater than the lower threshold and less than the upper threshold) using the tally function. We use I(result>= 0.407 & result <=0.715) to create an indicator variable that is 1 when result is in the interval and zero otherwise.

tally(∼I(result >= 0.407 & result <=0.715),

data=samp.mean, format="proportion")

I(result >= 0.407 & result <= 0.715)

TRUE FALSE

0.9548 0.0452

We see that roughly 95% of the sample means fall within 2 standard errors of the mean.

Central limit theorem

The sampling distribution of means (or sums) approaches that of a Normal distribution as the sample size increases. This result, given by the Central Limit Theorem, forms the basis of many formula-based methods of statistical inference and can be illustrated through simulation (Kwak & Kim, 2017). For example, the sampling distribution in our bear example is right-skewed for small sample sizes (e.g., n = 10; Fig. 1B); however, for samples of size 75, the sampling distribution is bell-shaped and well-approximated by a Normal distribution. Thus, rather than blindly trusting that the Central Limit Theorem applies, we can demonstrate it first hand by sampling repeatedly using different sample sizes. Additionally, the Central Limit Theorem guarantees for sufficiently large samples that the sampling distribution of $\overline{x}$ will be approximately Normal with 95% of the sample means falling within 2 standard errors of μ in the polar bear example. This result holds even though the population distribution is right-skewed and highly discrete, with individual bear counts taking on only integer values. Further, many estimators (including all maximum likelihood estimators) are calculated via a sum of independent measurements, which ensures that their sampling distributions are asymptotically Normal. Understanding this powerful result helps uncover why many back-of-the-envelope calculations use estimate ± 2SE to form approximate 95% confidence intervals.

Bootstrap confidence intervals

In real applications, we do not have access to data from the full population. Instead, we have but a single sample. When making inference to the population, via traditional parametric methods or with resampling-based methods, we must assume our sample data are representative of the population, as is assumed when a large number of observations are selected by simple random sampling. Then, we can use the distribution of values in our sample to approximate the distribution of values in the population. For example, we can make many, many copies of our sample data and use the resulting data set as an estimate of the whole population. With this estimated population in place, we could repeatedly sample from it, forming many data sets that are the same size as our original data set, and calculate statistics for each sample to approximate sampling distributions (Fig. 2). In practice, we do not actually need to make multiple copies of our sample data to estimate the population; instead, we form new data sets that are the same size as our original data set by sampling our original data with replacement, which effectively does the same thing. Sampling with replacement means that we select cases one at a time, and after each selection, we put the selected case back in the population so it can be chosen again. Thus, each observation in the original data set can occur zero, one, two, or more times in the generated data set, whereas it occurred exactly once in the original data set.

This process allows us to create a bootstrap sampling distribution to determine how much our estimates vary among repeated samples. We quantify this variability using the standard deviation of our estimates across repeated samples and refer to this standard deviation as our bootstrap standard error (BSE; Fig. 2). If the bootstrap sampling distribution is bell-shaped and centered on our sample statistic (indicating there is no sampling bias), then we can use this estimated standard error and the 2 SE rule for Normally distributed data to calculate a 95% confidence interval for our parameter of interest: estimate ± 2BSE. We can also think about repeating this process many times (i.e., collect data, estimate a parameter, use the bootstrap to calculate a SE and confidence interval), in which case we expect 95% of our confidence intervals to include the true population parameter of interest. That is, the bootstrap method for calculating a 95% confidence interval is what statisticians call well calibrated. In practice, we only have one data set, leading to a single confidence interval that either does or does not contain the true value of the population parameter, but knowing the method is well calibrated allows us to state that we are 95% confident that our interval contains the population parameter.

We illustrate this approach by considering data collected by the Minnesota Department of Natural Resources (MN DNR) to explore the potential impact of changing fishing regulations on the size distribution of northern pike (Esox lucius) in Medicine Lake, an approximately 460-acre lake in Beltrami County, MN. In 1989, the MN DNR instituted a slot limit of 22–30 inches in this lake (i.e., all caught fish within this size interval had to be released). We consider length data from 73 and 81 fish collected in trap nets in 1988 and 1993, respectively (before and after the fishing regulation was put in place). Importantly, these data come from only one lake, and many other factors may have changed between 1988 and 1993. Therefore, we must be cautious when interpreting any changes in the distribution of fish sizes (i.e., attributing the cause of length changes to the management regulation or generalizing results to other lakes). Nonetheless, we can ask, “How much did fish length, on average, change between 1988 and 1993 in Medicine Lake?” To address this question, we estimate the difference in the mean length of fish in the two samples. We also quantify our uncertainty in this estimated difference in means, recognizing that we would get a different estimate if we could go back and collect other samples of fish in those 2 years.

We begin by calculating the sample size, mean length in the 2 years, and the difference in sample means using the tally, mean, and diffmean functions in the mosaic library. Again, note that the functions in the mosaic library have a common syntax for data summaries involving two variables: goal(y∼x, data=Mydata). For example, mean(y∼x, data=Mydata) will calculate the mean of y for each level of a categorical variable x, and diffmean(y∼x, data=MyData) will calculate the difference in sample means when x is a categorical variable taking on only two levels. Here, we store the difference in means in an object named effect_size:

pikedat<-read.csv(here("data", "Pikedata.csv"))

tally(∼year, data=pikedat)

year

1988 1993

73 81

mean(length.inches∼year, data=pikedat)

1988 1993

18.6 21.2

(effect_size <- diffmean(length.inches∼year, data=pikedat))

diffmean

2.58

We estimate that the mean size of pike increased by roughly 2.6 inches between 1988 and 1993. To evaluate uncertainty in our estimated effect size, we explore the variability in the difference in means across bootstrapped samples. Here, we use the sample data from 1988 and 1993 as our estimate of the distribution of lengths in the population in 1988 and 1993, respectively. We repeatedly sample (10,000 times) from these estimated populations and calculate our sample statistic (the difference in means) for each of these bootstrapped data sets, making use of the do and resample functions in the mosaic library. Whereas the sample function, by default, samples without replacement, the resample function samples with replacement; the group argument is used to ensure that bootstrapping is conducted separately for each year.

bootdist<-do(10000)*{

diffmean(length.inches∼year, data=resample(pikedat, group=year))

}

Above, we stored the 10,000 differences in means as a variable named diffmean contained in the bootdist object that we created. We can look at the first five differences in means using the head function.

head(bootdist, n=5)

diffmean

1 3.27

2 1.19

3 2.32

4 3.11

5 2.92

Because the bootstrap distribution is bell-shaped and centered on our sample statistic (Fig. 3), we can estimate a confidence interval for the difference in mean length in 1993 relative to 1988 using the 2 SE rule (i.e., based on an unstated Normal distributional assumption):

(SE<-sd(∼diffmean, data=bootdist))

[1] 0.673

(CIa<-effect_size + c(−2, 2)*SE)

[1] 1.24 3.93

Alternatively, we could use quantiles of the bootstrap distribution, which does not require Normality (Fig. 2). The 0.025 quantile (or, equivalently, 2.5th percentile) refers to the value of x, such that 2.5% of the distribution falls below x. Here, we use the qdata function in the mosaic library to determine the 0.025 and 0.975 quantiles. These quantiles capture the middle 95% of the bootstrap distribution for the difference in means:

(CIb<-qdata(∼diffmean, data=bootdist, p=c(0.025, 0.975)))

quantile p

2.5% 1.26 0.025

97.5% 3.90 0.975

For bell-shaped bootstrap distributions, Normal-based and percentile-based bootstrap confidence intervals work well and will usually give similar results. For small samples or skewed distributions, better methods exist (Davison & Hinkley, 1997; Hesterberg, 2015; Puth, Neuhäuser & Ruxton, 2015). Nonetheless, simple examples like this can facilitate understanding confidence intervals and other measures of uncertainty.

Null distributions and p-values

Although null hypothesis testing has largely fallen out of favor in many disciplines, including ecology (Johnson, 1999; Hobbs & Hilborn, 2006), p-values, when fully understood, can play an important role in assessing the plausibility of statistical hypotheses (De Valpine, 2014; Murtaugh, 2014; Dushoff, Kain & Bolker, 2019). Further, it is critically important for biologists to be able to interpret the p-values they see in papers and in output generated by statistical software. To understand p-values, we must consider the sampling distribution of our statistic across repeated samples in the case where the null hypothesis is true. Simulation-based methods are also useful for understanding this concept, provided we can generate data sets consistent with the null hypothesis. To illustrate, we will consider data from an experiment used to test various hypotheses about the mating preferences of female sagebrush crickets (Cyphoderris strepitans) (Johnson, Ivy & Sakaluk, 1999). Males of this species will often allow a female to eat part of their hind wings when mating, which decreases the male’s attractiveness to future potential mates. To explore how diet might influence mating behaviors of females, Johnson, Ivy & Sakaluk (1999) randomly assigned 24 female crickets to either a low-nutrient (starved) or high-nutrient (fed) diet before putting them in a cage with a male cricket to mate. They then measured the waiting time before mating occurred. These data are contained in the abd package and were used by Whitlock & Schluter (2009) to introduce the Mann–Whitney U-test.

We can access the data using the data function after first loading the abd package:

library(abd)

data("SagebrushCrickets")

We then use the str function to explore the structure of the data set, finding that there are two variables: treatment, a factor variable with values fed or starved, and time.to.mating, a numerical variable containing the mating times for each cricket.

str(SagebrushCrickets)

‘data.frame’: 24 obs. of 2 variables:

$ treatment : Factor w/ 2 levels "fed", "starved": 2 2 2 2 2 2 2 2 2 2 ...

$ time.to.mating: num 1.9 2.1 3.8 9 9.6 13 14.7 17.9 21.7 29 ...

Our objective is to test whether the mean waiting times depend on the diet of the female: $H_0:{\mu }_{\mathrm{f}\mathrm{e}\mathrm{d}}={\mu }_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{v}\mathrm{e}\mathrm{d}}$ versus $H_A:{\mu }_{\mathrm{f}\mathrm{e}\mathrm{d}}\ne {\mu }_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{v}\mathrm{e}\mathrm{d}}$ where $H_0$ and $H_A$ stand for the null and alternative hypotheses, and ${\mu }_{\mathrm{f}\mathrm{e}\mathrm{d}}$ and ${\mu }_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{v}\mathrm{e}\mathrm{d}}$ represent the mean waiting times of crickets fed the high- and low-nutrient diets, respectively. It is natural to consider using the difference in sample means, ${\overline{x}}_{\mathrm{f}\mathrm{e}\mathrm{d}}-{\overline{x}}_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{v}\mathrm{e}\mathrm{d}}$ to conduct our hypothesis test. Given the small sample sizes (we have 13 and 11 cases in the fed and starved treatment groups, respectively), we cannot safely assume the Central Limit Theorem ensures that the sampling distribution of ${\overline{x}}_{\mathrm{f}\mathrm{e}\mathrm{d}}-{\overline{x}}_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{v}\mathrm{e}\mathrm{d}}$ is Normal. Further, the distribution of waiting times within each treatment group is far from Normal, and the measurements do not appear to be equally variable within each treatment group (Fig. 4). Thus, the assumptions of a parametric t-test are likely not met.

Although we could use a rank-based non-parametric test, such as the Mann–Whitney U-test (Whitlock & Schluter, 2009), we can gain additional insights into key statistical concepts (null distributions and p-values) if we construct our own randomization test, which also allows us to relax the Normality assumption.

We begin by estimating the mean waiting times for the two treatment groups (fed and starved) as well as the difference in sample means.

mean(time.to.mating∼treatment, data=SagebrushCrickets)

fed starved

36.0 17.7

(effect_size2<-diffmean(time.to.mating∼treatment, data=SagebrushCrickets))

diffmean

−18.3

We see that the mean time to waiting is 18.3 h longer for the crickets in the fed treatment group. A skeptic might point out that each time you repeat this experiment you will get a different value for this statistic (i.e., the difference in sample means). Perhaps the longer mean waiting time was just a fluke. Could that be the case? To find out, we need to determine what values we might expect to see for the difference in means when the null hypothesis is true. If the null hypothesis is true, then the labels (fed and starved) should not matter. So, to simulate a sample statistic (difference in means) that we might get if the null hypothesis were true, we can: (1) randomly shuffle the treatment variable among cases to form a new data set, and (2) calculate the difference in sample means. Using the do, diffmean, and shuffle functions in the mosaic library, we can repeat this process many times:

randomization_dist<-do(10000)*{

diffmean(time.to.mating∼shuffle(treatment), data=SagebrushCrickets)

}

Here, the shuffle function permutes the treatment variable, ensuring that, on average, the difference in sample means between the two treatment groups is equal to zero. Again, the do function stores the 10,000 differences in means in a variable named diffmean in the randomization_dist object that we created. Plotting this randomization distribution (Fig. 5), we find the range of statistics (difference in sample means) that we would expect to see if the null hypothesis were true. We can then use this information to determine the probability of observing an absolute difference in sample means as large or larger than we saw in our experiment, $|{\overline{x}}_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{v}\mathrm{e}\mathrm{d}}-{\overline{x}}_{\mathrm{f}\mathrm{e}\mathrm{d}}=18.3|$, if the null hypothesis is true. This probability is our two-sided p-value. We can approximate it using our randomization distribution by counting the number of times the difference in sample means was ≤−18.3 or ≥18.3 (equivalent to the number of times the absolute value of the difference was ≥ 18.3).

(p.value<- tally(∼abs(diffmean)>=18.3, data=randomization_dist,

format="proportion"))

abs(diffmean) >= 18.3

TRUE FALSE

0.126 0.874

Given a p-value of 0.126, we expect that roughly 13% of the time we would get an absolute difference in mean waiting times as large or larger than 18.3 if the null hypothesis were true. Thus, our result is perhaps not all that surprising. Still, our sample sizes were small, leading to differences in sample means that were highly variable from sample to sample. It would behoove us to also report a confidence interval to go with our result, so the reader has an understanding of the range of effect sizes that are plausibly supported by the data. Using a bootstrap, we find that we are 95% confident that the difference in population means is between −38.9 and 3.6.

boot_dist<-do(10000)*{

diffmean(time.to.mating∼treatment, data=resample(SagebrushCrickets))

}

qdata(∼diffmean, data=boot_dist, p=c(0.025, 0.975))

quantile p

2.5% −38.91 0.025

97.5% 3.61 0.975

Generalizations to other sample statistics

The methods illustrated here easily generalize to other statistics and data sets, provided cases can be assumed to be independent (e.g., there are no repeated measures on the same sample unit). For example, we could calculate a bootstrap distribution, and hence confidence interval, for a correlation coefficient using:

do(10000)*{cor(y∼x, data=resample(Mydata))}

We could determine an appropriate null distribution for testing whether a set of regression coefficients are all zero versus an alternative hypothesis that at least one is non-zero by shuffling the response variable:

do(10000)*{lm(shuffle(y)∼x1 + x2 + x3, data=Mydata)}

Shuffling the response variable breaks any association between y and the predictor variables, thus meeting the assumption of the null hypothesis. It is tempting to think that we could also determine an appropriate null distribution for testing whether a regression coefficient for x₂ is zero, while adjusting for two other variables (x₁ and x₃) using:

do(10000)*{lm(y∼x1 + shuffle(x2) + x3, data=Mydata)}

Although shuffling the values of x₂ for each case will break any relationship between x₂ and y, this process also breaks any association between x₂ and the other two predictor variables, x₁ and x₃. Thus, shuffling x₂ imposes additional unintended restrictions on the data beyond those required to ensure the null hypothesis is true. Alternative solutions have been proposed, which involve permuting residuals, either residuals of the response (Freedman & Lane, 1983; Anderson & Robinson, 2001), or the predictor at stake (Collins, 1987; Dekker, Krackhardt & Snijders, 2007). The latter approach is perhaps simplest, has good statistical properties, and is described below (see Appendix A for a simple coded example):

1. Fit a linear regression model relating $x_2$ to the other predictor variables in the model: lm(x2∼x1+x3, data=Mydata). The residuals from this model, which we label x2r, capture the part of $x_2$ that is not explained by the linear effect of other variables in the model (here, $x_1$ and $x_3$). If we replace x2 in our original model with x2r, fitting the following model: lm(y∼x1+x2r+x3, data=Mydata), the coefficient, standard error, t-statistic, and p-value for x2r will be identical to the coefficient for x2 in our original model.

2. Create an appropriate null distribution of t-statistics associated with the coefficient for x2r by shuffling x2r: do(10000)*{lm(y∼x1 + shuffle(x2r) + x3, data=Mydata)}. The test is based on the t-statistic from the above randomization distribution rather than the coefficient for x2r itself because the t-statistic ensures the permutation distribution does not depend on additional unknown parameters, whereas the regression coefficient does not have this property.

As the above example illustrates, solutions to more complex applications may require custom-written code (and additional care), but the underlying concepts remain the same. For confidence intervals, we need to generate bootstrapped data sets that mimic how we obtained the original sample data, calculate our sample statistic for each of these data sets, and then use the variability in those statistics to form our confidence interval. For testing null hypotheses, we need to generate data sets for which we have assured compliance with the null hypothesis and then compare our observed statistic with those simulated under that null hypothesis. In the words of Cobb (2007), “I like to think of (the core logic of statistical inference) as three Rs: randomize, repeat, reject. Randomize data production; repeat by simulation to see what’s typical and what’s not; reject any model that puts your data in its tail.”

For computationally intensive applications, it may not be possible to generate as many as 10,000 random data sets as in our examples. In such cases, it is customary and important to include the original sample statistic in the permutation or bootstrap distribution (Davison & Hinkley, 1997; Phipson & Smyth, 2010). Doing so ensures statistical hypothesis tests are exact even for small numbers of permutations (Phipson & Smyth, 2010), that is, across repeated tests with a significance level of $\alpha =0.05$, we will reject at most 5% of the time when the null hypothesis is true.

Bootstrapping versus permutations

Our simple examples above demonstrate how to form confidence intervals by bootstrapping and how to test statistical hypotheses using permutations. In truth, one can use bootstrapping to perform hypothesis tests and permutations to form confidence intervals (e.g., by “inverting” a hypothesis test so that the confidence interval includes only parameter values for which we would fail to reject the null hypothesis). Ideally, the chosen approach (resampling or permutation) should reflect the original study design. For randomized experiments, it makes sense to use permutations that reflect alternative outcomes of the randomization procedure. By contrast, bootstrap resampling may be more appropriate for observational data where the “randomness” is reflected by the sampling procedure (Lock et al., 2013).

Consider again the fish length example used to demonstrate bootstrap confidence intervals. If we wanted to conduct a statistical hypothesis test that mean fish length was the same in both years, we could make the null hypothesis true (by adding the difference between the two sample means to the length of fish in the first sample), then bootstrap. This approach is illustrated below. We begin by creating a new variable, length.null, by adding the difference in sample means to all fish lengths collected in 1988. We can accomplish this task using the mutate function in the dplyr library (Wickham et al., 2019). The ifelse function evaluates the first expression (year=="1988") and executes the next argument when true and the last argument when false. The difference in the mean of this new variable is zero (reflecting our null hypothesis).

delta<-diffmean(length.inches∼year, data=pikedat)

pikedat<-mutate(pikedat,

length.null=ifelse(year=="1988", length.inches+delta,

length.inches))

diffmean(length.null∼year, data=pikedat)

diffmean

0

We can then use the bootstrap to quantify variability in difference in sample means when the null hypothesis is true, from which we can calculate the p-value for our hypothesis test.

randfish<-do(10000)*{

diffmean(length.null∼year, data=resample(pikedat, group=year))

}

(p.value2<-tally(∼abs(diffmean)>=delta, data=randfish, format="proportion"))

abs(diffmean) >= delta

TRUE FALSE

0 1

None of the bootstrap differences in sample means was as large or larger than the difference we observed using our original data, resulting in a p-value $\approx$ 0.

This example raises another important point of consideration when choosing an appropriate method for testing null hypotheses, namely differences between bootstrap- and permutation-based null distributions. In the sagebrush cricket example, we permuted treatment labels to create our null distribution for the difference in sample means. By permuting treatment labels, we ensure the distribution of the response variable, $y$, is the same for both treatment groups. Thus, means, variances, and higher-order moments will all, on average, be the same for the two groups. By contrast, the bootstrap-based test ensures only that the sample means of the two groups are, on average, the same under the null hypothesis. Both methods should reject the null hypothesis in cases where the population means differ, provided our sample size is large enough. However, the permutation-based test may also reject the null hypothesis (with probability > α) if the population means are the same in the two groups, but the population variances differ. Thus, the bootstrap-based test may be preferred if the goal is to detect only differences in population means, allowing for the possibility that the population variances differ.

Case Study I: relaxing the assumptions of linear regression

Zuur et al. (2009) lead off their popular book, Mixed Effects Models and Extensions in Ecology with R, with the statement (p. 19), “In Ecology, the data are seldom modeled adequately by linear regression models. If they are, you are lucky.” Here, we consider one of the data sets from their book, containing measurements of species richness at 45 sampling stations in the Netherlands (five stations located within each of nine beaches, Janssen & Mulder, 2004, 2005). When plotting species richness against an index of exposure of each beach to waves and surf, we see that mean species richness decreases linearly with exposure level (Fig. 6A), but the distribution of residuals is right skewed and, thus, non-Normal (Fig. 6B). In addition, observations were clustered in space so we might expect two observations from the same beach to be more alike than two observations from two different beaches.

The rest of Zuur et al.’s book is devoted to describing regression methods that allow one to relax these assumptions; for example, using random effects to account for correlation and different statistical distributions to relax the Normality and constant variance assumptions. Without knowledge of these tools, however, is there a way to estimate uncertainty in the regression line as a summary of the relationship between species richness and exposure (over the observed range of the data)? Yes! We just need to adapt the bootstrap to mimic the way data were collected (i.e., cluster sampling, with multiple observations from each of several beaches). In this case, we need to sample beaches with replacement, keeping all observations from a beach when selected (we will refer to this approach as a cluster-level bootstrap). Doing so allows us to evaluate how much our estimates of regression parameters would change if we were to collect another sample of observations from a different set of nine beaches. Our estimate of uncertainty would not require assumptions of Normality or constant variance. Further, we would only need to assume that data from different beaches are independent (not that observations from the same beach are independent). Although it would be nice to have a larger sample of beaches to represent the population, alternative analysis methods (e.g., models appropriate for clustered data) would face a similar challenge of having to estimate among-beach variability from nine beaches.

Conceptually, this example is no more difficult to understand than the fisheries slot limit example. We are trying to determine how much our regression parameter estimates would vary if we were able to collect another sample of observations using the same sampling effort (nine beaches, each with five stations). Implementation becomes slightly more challenging than the previous examples because we can no longer rely on the resample function in the mosaic package (it does not accommodate cluster-level bootstrapping, but see Deen & De Rooij, 2019, for an alternative using the ClusterBootstrap package). The bootstrap distributions of the regression coefficients in this example are not symmetric, so it is beneficial to consider alternative bootstrap confidence interval procedures. It is important to note that not all confidence intervals perform equally well. Normal- and percentile-based intervals are simple but less accurate on average than intervals that attempt to account for bias and skewness of the sampling distribution (e.g., the “BCa” interval illustrated in Appendix B); accuracy is measured in terms of an interval estimator’s coverage rate, defined as the frequency with which the interval estimator includes the true population parameter. We include, as Supplemental Material, code necessary to implement the cluster-level bootstrap, demonstrating a range of confidence interval procedures implemented using the boot package (Canty & Ripley, 2019; Appendix B).

Case Study II: combining data sources: estimating uncertainty when the sampling distribution is unknown

As mentioned above, the Central Limit Theorem guarantees that the sampling distribution of many estimators, including all maximum likelihood estimators, will be Normally distributed as sample sizes become large (Casella & Berger, 2002). In practice, however, it may be difficult to derive an analytical expression for the standard error to use in conjunction with a Normal approximation. This is particularly true for estimation problems that require fitting multiple models to different response variables collected from the same or different sets of sample units. As an illustrative example, we consider data from Zicus, Rave & Fieberg (2006), used to evaluate the relative cost-effectiveness of single- and double-cylinder nesting structures for mallard (Anas platyrhynchos) ducks. To measure duckling productivity, Zicus, Rave & Fieberg (2006) fit a linear regression model to duckling counts collected from 110 nesting structures over an 8-year period. To estimate costs, they fit an additional discrete-time survival model to quantify the probability that a structure would fall over and hence require a visit (with associated maintenance costs) for continued use. Both models included predictors capturing cylinder type (single or double) and size of the wetland holding the structure. Zicus, Rave & Fieberg (2006) combined the output from the two models and estimated, for a range of wetland sizes and for both cylinder types, the expected cost of maintaining the structure divided by the expected number of ducklings produced as a measure of cost-effectiveness. Estimating the standard error of this cost-effectiveness measure was complicated for several reasons: (1) the quantity of interest was a ratio and therefore not a linear function of model parameters, (2) there was a complex dependence structure resulting from repeated measures on the same set of structures, which were also shared across the two models; and (3) the residuals from the duckling model were not Normally distributed.

Fortunately, using a bootstrap to estimate a confidence interval for cost-effectiveness was relatively straightforward. We just generated new bootstrapped data sets, mimicking the original sampling design, and repeated the analysis using these new data sets. As illustrated in Appendix C, we: (1) resampled nest structures with replacement; (2) fit both models to the bootstrapped data sets; and (3) used the fitted models to estimate the expected number of ducklings, expected structure survival (and hence cost), and the ratio of expected cost to the expected number of ducklings for each cylinder type at different wetland sizes. This process allowed us to calculate confidence intervals for the cost effectiveness of both types of nest structures across a range of wetland sizes (Appendix C). We have used this same general approach elsewhere to quantify uncertainty in estimated growth rates from population models when parameters are estimated from multiple data sources (e.g., Ellner & Fieberg, 2003; Fieberg et al., 2010; Lenarz et al., 2010).

One important consideration in this example is that our estimator of cost-effectiveness, a ratio of expected means, may be biased. We can also use the bootstrap to estimate Bias = $E[\hat{\theta }]-\theta$, where $E[\hat{\theta }]$ refers to the mean of the sampling distribution of $\hat{\theta }$, our estimator of the unknown parameter $\theta$. To understand how this works, consider that we use the bootstrap to mimic sampling from our population: the bootstrap samples relate to the original sample in much the same way as the original sample relates to the population (Hall, 1988; Fig. 2). Thus, to estimate bias, we can compare the mean of our bootstrap statistics, ${\theta }_b^{\ast }$, to the same statistic calculated using our original sample, $\hat{\theta }$, as though the sample were the true population. In Appendix C, we estimate the bias of our cost-effectiveness measure for both deployment types and across a range of wetland sizes, showing that the bias is small relative to the standard error.

Case Study III: evaluating model-selection uncertainty

The bootstrap can be used to quantify model uncertainty in addition to sampling variation (Buckland, Burnham & Augustin, 1997; Fieberg & Johnson, 2015; Harrell, 2016). Consider, for example, the common situation in which a researcher desires to develop a predictive model by selecting a limited number of explanatory variables from some larger set (e.g., using AIC statistics associated with various sub-models). Several models may provide nearly equally good fits to the data, in which case we might expect a different model to come out as “best” each time we collect a new data set (where “best” is determined via backwards stepwise selection using AIC). We can explore whether this is the case by fitting the same set of candidate models to each of several bootstrapped data sets.

In this case study, we explore model-selection uncertainty arising from backwards stepwise selection using AIC (Appendix D). We begin by fitting a linear regression model to abundance data of longnose dace (Rhinichthys cataractae) collected from the Maryland Biological Stream Survey (downloaded from http://www.biostathandbook.com/multipleregression.html). Starting with six predictors (measurements of various stream attributes), we use the stepAIC function in the MASS package in R (Venables & Ripley, 2013) to sequentially eliminate predictors until we can no longer reduce AIC further by dropping any of the remaining predictors. Implementing this process resulted in a reduced model containing three variables: acreage drained by the stream, maximum depth, and NO₃. We can then use functions in the rms (regression modeling strategies) package in R (Harrell, 2019) to explore variability in the best model as determined by backwards stepwise selection. First, we refit the full model using ols (this function is equivalent to lm, but it stores additional information with the model fit). We then use the validate function to choose a best model by applying backwards selection to each of 100 bootstrapped data sets (i.e., validate fits a series of models, using backwards selection to choose the best model for each bootstrapped data set). Applying this process resulted in 31 different best models. The most frequently chosen model, selected in 24 of the 100 cases, included only acreage and NO₃. The model originally chosen using stepAIC rose to the top in only nine of the 100 bootstrapped data sets. One bootstrap sample resulted in a model with all six predictors, and five bootstrap samples led to a model with a single predictor (maximum depth in four instances and NO₃ in one instance). In addition to demonstrating the high degree of uncertainty associated with choosing a best model, these bootstrap results could be used to calculate model-averaged predictions (using the frequency with which a model was selected to determine model weights; Buckland, Burnham & Augustin, 1997).