STAC32 Assignment 4

Packages

library(tidyverse)

You are expected to complete this assignment on your own: that is, you may discuss general ideas with others, but the writeup of the work must be entirely your own. If your assignment is unreasonably similar to that of another student, you can expect to be asked to explain yourself.

If you run into problems on this assignment, it is up to you to figure out what to do. The only exception is if it is impossible for you to complete this assignment, for example a data file cannot be read. (There is a difference between you not knowing how to do something, which you have to figure out, and something being impossible, which you are allowed to contact me about.)

You must hand in a rendered document that shows your code, the output that the code produces, and your answers to the questions. This should be a file with .html on the end of its name. There is no credit for handing in your unrendered document (ending in .qmd), because the grader cannot then see whether the code in it runs properly. After you have handed in your file, you should be able to see (in Attempts) what file you handed in, and you should make a habit of checking that you did indeed hand in what you intended to, and that it displays as you expect.

Hint: render your document frequently, and solve any problems as they come up, rather than trying to do so at the end (when you may be close to the due date). If your document will not successfully render, it is because of an error in your code that you will have to find and fix. The error message will tell you where the problem is, but it is up to you to sort out what the problem is.

Livestock and fluoride

Livestock such as cows “graze” (that is, eat grass). Do they absorb pollution from the environment? In a study, some livestock were grazed in an area known to be polluted with fluorine, and some other livestock were grazed in an area not polluted with fluorine. At the end of the study, the fluorine content of the livestock’s urine was measured. The sample sizes, however, are very small. The data are in http://datafiles.ritsokiguess.site/livestock-fluorine.csv. The column pollution indicates whether the animal was grazing in a Polluted or an Unpolluted area, and the column fluoride indicates the amount of fluorine (in parts per million).

  1. (1 point) Read in and display (most of) the data.

The usual:

my_url <- "http://datafiles.ritsokiguess.site/livestock-fluorine.csv"
livestock <- read_csv(my_url)
Rows: 12 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): pollution
dbl (1): fluoride

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
livestock

I think the name livestock describes what the data set is, without duplicating any column names. There are only 12 animals altogether.

Extra: these are the data set xmp15.04 from the Devore7 package. These came in the wrong shape:

Devore7::xmp15.04

so I rearranged things into the right shape using pivot_longer (which we will see later in the course):

Devore7::xmp15.04 %>% 
  pivot_longer(everything(), names_to = "pollution", values_to = "fluoride") %>% 
  drop_na()

The original data has some “structural missings” caused by there being seven animals grazing in polluted areas and only five in unpolluted areas. With the original shape, the two columns have to be the same length, so one of the columns has acquired some missing values that are not actually missing data but are a consequence of the layout. I no longer need these in the new layout, hence the drop_na at the end.

  1. (3 points) Make a suitable plot to assess the normality of the fluoride measurements for each level of pollution. What do you see?

Assessing the normality means a normal quantile plot, using facets to do each level of pollution separately:

ggplot(livestock, aes(sample = fluoride)) +
  stat_qq() + stat_qq_line() +
  facet_wrap(~ pollution)

As far as we can tell from such small samples, the Polluted values are close to normal, and the Unpolluted values have a low outlier (and maybe a high one too).

My take on the Polluted values is that the highest and lowest values are close to the line, and there is no apparent curve overall. Non-normality generally shows up in the tails rather than the middle of the distribution. On the Unpolluted data, the extreme observations at both ends seem far from the line.

Extra: the problem is that with samples as small as these, the values can look rather non-normal even if they actually came from a normal distribution. One way to assess this is to make what is called a “lineup” of several normal quantile plots that are actually normal data with the same sample size (say, 5 like the Unpolluted data):

tibble(sim = 1:9) %>% 
  rowwise() %>% 
  mutate(my_sample = list(rnorm(5))) %>% 
  unnest_longer(my_sample) %>% 
  ggplot(aes(sample = my_sample)) +
    stat_qq() + stat_qq_line() +
    facet_wrap(~ sim, scales = "free")

These nine plots are all genuinely normal data, and then the question becomes “does my actual data look definitely different from these?” For example, plot #8 here looks quite a lot like our Unpolluted data, so maybe our Unpolluted data are not as far as you might think from what normal data might look like.

The real way to do a lineup is to make eight fabricated normal quantile plots from actual normal data, add your data to it, and then shuffle the results. Then you ask “which one of these looks different from the others?”, and you get your friends to do the same, and if more of you than chance can pick out your actual data, then your actual data really is different.

  1. (3 points) Obtain and make a normal quantile plot of the bootstrap sampling distribution of the sample mean for the observations from Polluted areas. Describe what you see. (Here, and in the next question, do not use the bstrap package.)

Obtain just the fluorine values for the Polluted areas first:

livestock %>% filter(pollution == "Polluted") -> polluteds
polluteds

Then use these to obtain the bootstrap sampling distribution. At this point you might like to do a set.seed first (before any code that has random number generation in it):

set.seed(457299)

and then

tibble(sim = 1:10000) %>% 
  rowwise() %>% 
  mutate(my_sample = list(sample(polluteds$fluoride, replace = TRUE))) %>% 
  mutate(my_mean = mean(my_sample)) %>% 
  ggplot(aes(sample = my_mean)) + stat_qq() + stat_qq_line()

This looks very close to normal to me (not surprising since the original distribution was close to normal). The lower tail is, if anything, a little short.

  1. (3 points) Obtain and make a normal quantile plot of the bootstrap sampling distribution of the sample mean for the observations from Unpolluted areas. Describe what you see.

(I copied and pasted all my code from the previous question and carefully edited it.)

Obtain just the fluorine values for the Unpolluted areas first:

livestock %>% filter(pollution == "Unpolluted") -> unpolluteds
unpolluteds

Then use these to obtain the bootstrap sampling distribution:

tibble(sim = 1:10000) %>% 
  rowwise() %>% 
  mutate(my_sample = list(sample(unpolluteds$fluoride, replace = TRUE))) %>% 
  mutate(my_mean = mean(my_sample)) %>% 
  ggplot(aes(sample = my_mean)) + stat_qq() + stat_qq_line()

This is skewed to the left, but only slightly (and, given the low outlier in the original distribution, you might have been expecting things to come out much worse). I think there is a definite curved shape to the plot: the middle points are on or above the line, while the ones in the tails, both ends, are almost all below the line.

My guess is that the actual sample is not that different from what might have come from a normal distribution with equally small sample size, so the sampling distribution is more like a normal than you would expect.

  1. (2 points) Comment briefly on the appropriateness of using a two-sample \(t\)-test to compare the mean fluorine levels for the two groups of livestock.

Make a call about the normality of the bootstrap sampling distributions in each of the previous two questions. Remember that these two do need to be close to normal, in an absolute sense (because the bootstrap sampling distribution is accounting for sample size). If you are going to say that the two-sample \(t\) is appropriate, you need to make the case that both distributions are close enough to normal.

So: if you think that both distributions are normal enough, \(t\) is good. If you can find a problem with one of them, that’s enough to invalidate the \(t\)-test.

My take is that the bootstrap sampling distribution of the sample mean of the Unpolluted values is definitely skewed left and not close enough to normal, and therefore that the two-sample \(t\) is no good.

Extra: you might be wondering whether we can get both bootstrap sampling distributions at once, rather than one at a time. The problem is in sampling the right thing at the right time, but something like this makes it easier:

livestock %>% 
  nest_by(pollution)

This seems to have made things worse, but now the list-column data contains the fluoride values for the Polluted areas (first row) and the Unpolluted areas (second row). With a little care, we can mimic what we did for one group at a time for both groups at once. I’ll illustrate with three simulations first, so you can see what’s going on, and scale up later:

livestock %>% 
  nest_by(pollution) %>% 
  mutate(sim = list(1:3)) %>% 
  unnest_longer(sim)

Even though the column data looks the same all the way down, it actually contains the fluoride values for the right pollution value (on the left), so now, when we do the bootstrap sampling, it will sample from the right thing:

livestock %>% 
  nest_by(pollution) %>% 
  mutate(sim = list(1:3)) %>% 
  unnest_longer(sim) %>% 
  rowwise() %>% 
  mutate(my_sample = list(sample(data$fluoride, replace = TRUE)))

The first three bootstrap samples have 7 observations in them, because there were seven observations from Polluted areas, and the others have 5 (five observations from Unpolluted areas). So this looks as if it has done the right thing.

Just to check that this is working, let’s finish the procedure before we scale it up: work out the mean of each bootstrap sample, and then make normal quantile plots of the means, facetted by the value in pollution:1

livestock %>% 
  nest_by(pollution) %>% 
  mutate(sim = list(1:3)) %>% 
  unnest_longer(sim) %>% 
  rowwise() %>% 
  mutate(my_sample = list(sample(data$fluoride, replace = TRUE))) %>% 
  mutate(my_mean = mean(my_sample)) %>% 
  ggplot(aes(sample = my_mean)) +
    stat_qq() + stat_qq_line() +
    facet_wrap(~ pollution, scales = "free")

The graphs look kind of dopey with only three points each, but we have three different bootstrap sample means from each sample, so it looks as if we are good.2

Now, hold your breath and scale it up. This is literally a matter of copying and pasting the above code and changing the 3 to 10,000 without making any other changes:

livestock %>% 
  nest_by(pollution) %>% 
  mutate(sim = list(1:10000)) %>% 
  unnest_longer(sim) %>% 
  rowwise() %>% 
  mutate(my_sample = list(sample(data$fluoride, replace = TRUE))) %>% 
  mutate(my_mean = mean(my_sample)) %>% 
  ggplot(aes(sample = my_mean)) +
    stat_qq() + stat_qq_line() +
    facet_wrap(~ pollution, scales = "free")

and you get more or less the same message as before: the Polluted distribution has a short lower tail but is otherwise close to normal, and the Unpolluted distribution is slightly left-skewed.

I don’t know how much you think the extra work here is worth it. After all, doing the two groups separately is no more than selecting the data for each group, then doing an ordinary one-sample bootstrap. Then, for the second one, copy-paste the code and change some variable names. There is a bit of care needed in the editing, but not that much. If you had three or four groups (as in an ANOVA), then you might start thinking (by the time you get to about the third copy-paste) “I wonder if I can do these all at once?”, and then you can use the procedure above, which works for any number of groups.

This, by the way, is what bstrap does behind the scenes for more than one group. Once again, you can use this package to check your work, but don’t hand it in:

library(bstrap)
with(livestock, plot_several(fluoride, pollution))

Power when data are uniformly distributed

A (continuous) random variable on the interval from min to max has a uniform distribution if it has equal probability of taking any value in that interval. The mean of the distribution is the mean of min and max. In R, the function runif draws a random sample from a uniform distribution. It has three inputs, in order: (i) the number of random observations to draw, (ii) the value of min, (iii) the value of max.

  1. (3 points) Draw a random sample of 100 observations from a uniform distribution on the interval from 0 to 1. Demonstrate that your observations have approximately the correct mean and their distribution has approximately the correct shape.

Here, min is 0 and max is 1, and the sample size is 100. I’m going to put these in a dataframe because I want to do some plotting and some calculations later. Use any name you like for the dataframe and the column in it. I’m using “generic” ones:

d <- tibble(x = runif(100, 0, 1))
d

The distribution should have mean \((0+1)/2 = 0.5\) and a flat-topped shape because any value is equally likely. The sample mean is

d %>% summarize(mean_x = mean(x))

which is close to 0.5 (to within randomness), and the distribution shape can be assessed from a histogram with not too many bins:

ggplot(d, aes(x = x)) + geom_histogram(bins = 6)

The best answer explains why this is not completely flat-topped: this uniform distribution only has values between 0 and 1, and the outermost bins go beyond these limits, so they have fewer observations in them. You don’t need to fix this, but you do need to explain what happened, because your reader is expecting a distribution that has a flat top, and they will want to know why this one does not.

Graphs that do not demonstrate that the distribution has the right shape:

  • a one-sample boxplot (does not show shape beyond symmetry)
  • a normal quantile plot (we are not assessing normality).

Extra: you’ll have to do some research (cite your sources) to find out how to improve this histogram. I looked in the help for geom_histogram (type ?geom_histogram in the console) to find that you actually have (at least) three options for specifying bins on a histogram:

  • bins to specify the number of bins
  • binwidth to specify the width of the bins
  • breaks to specify exactly where the bins should be.3

The last one of these, breaks, does exactly the job we want: we can specify bins that start exactly at 0 and end exactly at 1. For example, to put the bin boundaries at 0, 0.2, 0.4, …, 1, we can do this:

bin_breaks <- seq(0, 1, 0.2)
bin_breaks
[1] 0.0 0.2 0.4 0.6 0.8 1.0

I already knew about seq, but you probably don’t. You can find it by throwing something like “r how to specify sequence” into your favourite search engine. The first link I got back was https://www.statology.org/seq-function-in-r/, and this tells you that seq is what to use, and how to use it. This site often comes up (for me) in web searches, and I have found its results to be generally reliable. Usually you use the first three inputs, which are the start, the end, and the step from one value to the next.

Having found out about that, you then feed this into breaks in geom_histogram:

ggplot(d, aes(x = x)) + geom_histogram(breaks = bin_breaks)

This does indeed (to within randomness) have a flat top, showing that all values are about equally likely. Your graph won’t have the tallest bars in the same place as mine, but they should all be about the same height.

I said before that a normal quantile plot is not the right thing, but in fact you can make a quantile plot for any distribution. To find out how, look in the help for stat_qq to find that both it and stat_qq_line have an input distribution. The default distribution is the normal4 which is specified as stats::qnorm. To disentangle that, look in the help for qnorm to see that there are actually all of:

  • dnorm: normal probability density function (the thing that’s shaped like a bell)
  • pnorm: probability that a normal random variable is less than something (as in the tables you used in your first course)
  • qnorm: the value (of \(z\)) that has a certain probability of being less than it (using your normal tables backwards to find the \(z\) that goes with a given probability)
  • rnorm: random values from a normal distribution

In R, many different distributions are available: they each have a “base name” like norm for normal, binom for binomial and so on, onto which you add a d, p, q, or r as needed. We used runif to generate our random uniformly-distributed values, so the “base name” is unif and we must need qunif to draw a uniform quantile plot. In that, stat_qq and stat_qq_line work independently, so you specify the distribution in both places:

ggplot(d, aes(sample = x)) + 
  stat_qq(distribution = stats::qunif) + 
  stat_qq_line(distribution = stats::qunif)

These points follow the line well, so they have close to a uniform distribution, as they should.

In case you are curious (I was), here’s what a normal quantile plot for these data looks like:

ggplot(d, aes(sample = x)) + 
  stat_qq() + stat_qq_line()

Compared to a normal, a uniform distribution has short tails: this uniform distribution literally cannot go beyond 0 or 1, and you see that a normal distribution (with the same quartiles) would go up and down much further.

  1. (3 points) Using data from uniform distributions with min = 0, estimate (by simulation) the power of a \(t\)-test to reject a mean of 0.5 in favour of a two-sided alternative, when the mean is actually 0.6, using a sample size of \(n = 20\). Hint: if the mean is 0.6, what is max?

Sample from the truth, and test whether the mean is 0.5. “The truth” here means a uniform distribution with min of 0 and a mean of 0.6, which means a max of 1.2 (since \((0 + 1.2) / 2 = 0.6\)):

tibble(sim = 1:1000) %>% 
  rowwise() %>% 
  mutate(my_sample = list(runif(20, 0, 1.2))) %>% 
  mutate(my_test = list(t.test(my_sample, mu = 0.5))) %>% 
  mutate(p_value = my_test$p.value) %>% 
  count(p_value <= 0.05)

My estimated power is \(222/1000 = 0.222\).

This is not especially high (I think because the null and the truth are not very far apart, and so it will take a large sample to distinguish them with any reliability).

  1. (3 points) By experimentation, find the sample size that is required for this test to have power 0.60. Demonstrate a well-motivated process. Hint: it is enough to get power between 0.55 and 0.65.

This means reusing your code of the previous part with different sample sizes until you find one that hits a power between 0.55 and 0.65. By “well-motivated” I mean that you describe what you are trying next and why.

So far, we know that a sample size of 20 goes with a power of 0.222 (to within the accuracy of the simulation). We want a bigger power, so we need a bigger sample size, but we have at this point no idea how much bigger. All we can do is to try something and see what happens. How about \(n = 30\)? Copy and paste your previous code, but change the sample size from 20 to 30:

tibble(sim = 1:1000) %>% 
  rowwise() %>% 
  mutate(my_sample = list(runif(30, 0, 1.2))) %>% 
  mutate(my_test = list(t.test(my_sample, mu = 0.5))) %>% 
  mutate(p_value = my_test$p.value) %>% 
  count(p_value <= 0.05)

That improved the power by some, to 0.345, but not enough, so the sample size needs to be bigger still. \(n = 50\)?

tibble(sim = 1:1000) %>% 
  rowwise() %>% 
  mutate(my_sample = list(runif(50, 0, 1.2))) %>% 
  mutate(my_test = list(t.test(my_sample, mu = 0.5))) %>% 
  mutate(p_value = my_test$p.value) %>% 
  count(p_value <= 0.05)

0.495, so still not quite there yet. What about \(n = 75\)?

tibble(sim = 1:1000) %>% 
  rowwise() %>% 
  mutate(my_sample = list(runif(75, 0, 1.2))) %>% 
  mutate(my_test = list(t.test(my_sample, mu = 0.5))) %>% 
  mutate(p_value = my_test$p.value) %>% 
  count(p_value <= 0.05)

Now we have a power, 0.692, that is too high, so we need a sample size that is a little smaller.5 So we should go down a bit. 65 seems like a good guess:

tibble(sim = 1:1000) %>% 
  rowwise() %>% 
  mutate(my_sample = list(runif(65, 0, 1.2))) %>% 
  mutate(my_test = list(t.test(my_sample, mu = 0.5))) %>% 
  mutate(p_value = my_test$p.value) %>% 
  count(p_value <= 0.05)

and that just hits the target with a power of 0.643.

You can stop there, but I’m guessing that 60 is actually closer to the mark, so just to confirm that:

tibble(sim = 1:1000) %>% 
  rowwise() %>% 
  mutate(my_sample = list(runif(60, 0, 1.2))) %>% 
  mutate(my_test = list(t.test(my_sample, mu = 0.5))) %>% 
  mutate(p_value = my_test$p.value) %>% 
  count(p_value <= 0.05)

and so it is, with estimated power 0.611.

Your process will be different from mine, but by “well-motivated” I mean that you sensibly use previous guesses at the sample size to inform your next guess, without just leaping to \(n = 60\): that is, you try a sample size bigger than 20 (anything) at the first step, and use the resulting power to suggest your next sample size, repeating as necessary. This kind of open-ended exploration is something I can reasonably give you on an assignment but not so much on an exam. I think it is hard to justify going straight to \(n = 60\) on your second step (unless you knew that was already a reasonable answer); something between about 30 and 50 seems like a better next guess at the sample size, or you might try something like \(n = 100\) expecting to overshoot. If the grader sees a sensible process with three or more steps (including your initial simulation for \(n = 20\)) with a properly motivated sample size each time, finally reaching a power between 0.55 and 0.65, you’ll be good.

Extra:

I gave you a range to hit because power by simulation has randomness in it and you can only hope to get so close. On Worksheet 5, in the Extra to question 8, I showed you how you can get a confidence interval for the true power based on the randomness in your simulations, using prop.test. With 1000 simulations, this is typically close to your estimate plus or minus 3 percentage points (with 95% confidence). Here’s \(n = 65\):

prop.test(643, 1000)

    1-sample proportions test with continuity correction

data:  643 out of 1000, null probability 0.5
X-squared = 81.225, df = 1, p-value < 2.2e-16
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
 0.6123020 0.6725851
sample estimates:
    p 
0.643

The target 0.60 is not quite in this interval, so this sample size is a bit too big. \(n = 60\) gives

prop.test(611, 1000)

    1-sample proportions test with continuity correction

data:  611 out of 1000, null probability 0.5
X-squared = 48.841, df = 1, p-value = 2.776e-12
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
 0.5799086 0.6412279
sample estimates:
    p 
0.611

and now 0.60 is squarely in this interval. What about \(n = 57\)?

tibble(sim = 1:1000) %>% 
  rowwise() %>% 
  mutate(my_sample = list(runif(57, 0, 1.2))) %>% 
  mutate(my_test = list(t.test(my_sample, mu = 0.5))) %>% 
  mutate(p_value = my_test$p.value) %>% 
  count(p_value <= 0.05)
prop.test(580, 1000)

    1-sample proportions test with continuity correction

data:  580 out of 1000, null probability 0.5
X-squared = 25.281, df = 1, p-value = 4.956e-07
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
 0.5486573 0.6107204
sample estimates:
   p 
0.58

0.60 is in this one as well. So given the randomness present in the simulation, both 57 and 60 are supportable sample sizes, in that they both produced (for me) confidence intervals for the true power that contained the target power of 0.60; at this level of accuracy there is no reason to prefer one over the other.

To nail things down more accurately, the idea is to run more simulations, like 10,000. Trying again with \(n = 57\):

tibble(sim = 1:10000) %>% 
  rowwise() %>% 
  mutate(my_sample = list(runif(57, 0, 1.2))) %>% 
  mutate(my_test = list(t.test(my_sample, mu = 0.5))) %>% 
  mutate(p_value = my_test$p.value) %>% 
  count(p_value <= 0.05)
prop.test(5566, 10000)

    1-sample proportions test with continuity correction

data:  5566 out of 10000, null probability 0.5
X-squared = 127.92, df = 1, p-value < 2.2e-16
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
 0.5467932 0.5663631
sample estimates:
     p 
0.5566

and at this point we see with 95% confidence that the power with a sample size of 57 is too small.

  1. (2 points) Re-do the previous two questions using power.t.test (find the power for a sample size of 20, and find the sample size for a power of 0.60.) Are the answers close to what you got by simulation for uniformly-distributed data? Hint: the standard deviation of a uniform distribution is (max minus min), all divided by the square root of 12.

The appropriate SD is the one for the true distribution, that is

u_sd <- (1.2 - 0) / sqrt(12)
u_sd
[1] 0.3464102

and hence for the power

power.t.test(n = 20, delta = 0.6 - 0.5, sd = u_sd, type = "one.sample")

     One-sample t test power calculation 

              n = 20
          delta = 0.1
             sd = 0.3464102
      sig.level = 0.05
          power = 0.2315831
    alternative = two.sided

and for the sample size

power.t.test(power = 0.6, delta = 0.6 - 0.5, sd = u_sd, type = "one.sample")

     One-sample t test power calculation 

              n = 60.7275
          delta = 0.1
             sd = 0.3464102
      sig.level = 0.05
          power = 0.6
    alternative = two.sided

The power of 0.231 is very close to my simulated 0.222, and the sample size of 61 is very close to my best answer of 60. Your answers in this question from power.t.test should be the same as mine, but your simulations from earlier will probably be different. Adjust accordingly, but I expect you to get answers that are close without necessarily being exactly the same.

Extra: This indicates that sampling from a uniform distribution is almost exactly the same as sampling from a normal distribution when it comes to power. This is because it actually takes a very small sample size from a uniform distribution to get a sampling distribution of the sample mean that is indistinguishable from normal. Let’s try samples of size 10. Note that we are sampling from a distribution now, rather than a sample as we would for a bootstrap:

tibble(sim = 1:10000) %>% 
  rowwise() %>% 
  mutate(my_sample = list(runif(10, 0, 1.2))) %>% 
  mutate(my_mean = mean(my_sample)) %>% 
  ggplot(aes(sample = my_mean)) + stat_qq() + stat_qq_line()

10 is not a big sample by any means, but the sampling distribution of the sample mean from a uniform distribution is indistinguishable from normal. That is to say, even as small a sample size as 10 is big enough to overcome the non-normality of the uniform distribution.

In an earlier Extra, we saw that the uniform distribution has short tails compared to the normal. The investigation here suggests that having short tails is much less of a problem than having long tails (or outliers) in that the Central Limit Theorem kicks in for much smaller sample sizes when the distribution has short tails, compared to when it has long tails.

Footnotes

  1. I used scales = "free" here because I want the best look at whether these are normal; I am not interested here in anything like comparing spreads.↩︎

  2. The first time I did this, the three bootstrap means were all the same, and I scratched my head a bit, drank some coffee, and realized my omission: bootstrap samples are with replacement, and I had forgotten the replace = TRUE, so that all I had done was to resample all of the original sample each time, and of course the mean was the same every time I did that.↩︎

  3. This one is actually part of stat_bin which does the counting of how many observations fall into each bin. You normally don’t call that directly yourself, but you let geom_histogram call it on your behalf when it draws the histogram. When we learn about writing functions (near the end of this course), we learn about ..., which means that any inputs to a function that it doesn’t recognize are passed on to something else. In the Usage section of the help, you see that geom_histogram has a ..., so that breaks is passed on to stat_bin.↩︎

  4. Which is why stat_qq() with no inputs gives you a normal quantile plot.↩︎

  5. I actually overshot on purpose to show you what happens if your sample size is too big and gets you more power than you want.↩︎