7.2.1 Visualize the distribution

As usual, we begin by visualizing the actual data. As a reminder, we will assume for the time being that these data constitute the entire population of interest. The histogram below shows the population distribution of nightly hours of sleep:

nhanes %>%
    ggplot(aes(x = SleepHrsNight)) +
    geom_histogram(binwidth = 1)

7.2.2 Is the normal distribution a good model?

Next, we want to know whether the normal distribution will be a good model of the population distribution we just visualized. To do this, we will draw a curve representing the normal distribution on top of the histogram we just made.

But remember that the normal distribution is just a shape. Its center and spread are determined by two population parameters: the mean (\(\mu\)) and standard deviation (\(\sigma\)). What should those parameters be?

Let’s first take a wild guess and say they should be \(\mu = 8\) and \(\sigma = 1\). We can draw this distribution using the chunk of code below:

nhanes %>%
    ggplot(aes(x = SleepHrsNight)) +
    geom_histogram(aes(y = ..density..), binwidth = 1) +
    stat_function(fun = dnorm, args = list(mean = 8, sd = 1), color = 'darkred')

Notice that we needed to add aes(y = ..density..) to our histogram line, for reasons which will be clear in the next paragraph. The final line lets us draw a function on our graph. The function is called dnorm for “density of a normal distribution”. We had to say what the arguments of that function should be, which are a list in which we say what the mean and sd of the normal distribution should be.

The whole “density” thing comes from this: A histogram shows absolute counts. But a normal distribution only specifies the relative frequency of different values. The relative frequency is called the “density”. This is what the dnorm function gives us for the normal distribution. So adding aes(y = ..density..) to the histogram line told R to show the density (relative frequency) rather than the absolute counts.

Anyway, it is clear that a mean of 8 hours and standard deviation of 1 hour do not make a good fit.

nhanes %>%
    ggplot(aes(x = SleepHrsNight)) +
    geom_histogram(aes(y = ..density..), binwidth = 1) +
    stat_function(fun = dnorm, args = list(mean = ___, sd = ___), color = 'darkred')

7.2.3 Bad days, bad models?

As part of the same survey, people were asked how many days out of the last 30 would they consider their mental health to have been poor. This is recorded under the variable named DaysMentHlthBad.

___ %>%
    ggplot(aes(x = ___)) +
    geom_histogram(binwidth = 1)

7.3.1 The “true” proportion

Because we are assuming that the nhanes data represent the entire population we are interested in, we can find the “true” value of the population parameter \(\pi_{\text{Marijuana}}\), that is, the proportion of people who have ever tried marijuana. We can use the code below:

nhanes %>%
    specify(response = Marijuana, success = "Yes") %>%
    calculate(stat = "prop")

## Response: Marijuana (factor)
## # A tibble: 1 × 1
##    stat
##   <dbl>
## 1 0.585

7.3.2 Simulating many samples

Last session, we used bootstrapping to simulate what would happen if we collected many samples from a population. We did that by using our sample to estimate the whole population. Now, we have access to that whole population. We can therefore draw as many samples of a given size as we want!

The following chunk of code randomly samples 5 people from the population and gets R to remember that sample under the name sample_size5.

sample_size5 <- nhanes %>%
    slice_sample(n = 5)

This is what that sample looks like:

Notice that we still have all of those people’s responses to each question on the NHANES survey. That means we can get a summary statistic for the Marijuana variable from this sample, just like we did from the population.

sample_size5 <- nhanes %>%
    slice_sample(n = 5)

___ %>%
    specify(response = ___, success = "___") %>%
    calculate(stat = "___")

As we’ve seen, the great thing about computers is that they can do boring repetitive things for us very quickly. Now, we will draw 1000 samples of size 5 from the population and get R to remember them under the name samples_size5.

samples_size5 <- nhanes %>%
    rep_slice_sample(n = 5, reps = 1000)

Here are the first three samples in our new samples_size5:

7.3.3 Summary statistics for each sample

Notice that the different samples are indexed using the new variable replicate. We can use this to get the proportion of marijuana triers in each sample using the old group_by routine. We will tell R to remember these sample proportions under the name sample_props_mari_size5:

sample_props_mari_size5 <- samples_size5 %>%
    group_by(replicate) %>%
    summarize(stat = mean(Marijuana == "Yes"))

Notice that the proportions calculated from each sample will tend to vary. We can visualize this using a histogram:

sample_props_mari_size5 %>%
    ggplot(aes(x = stat)) +
    geom_histogram(aes(y = ..density..), binwidth = 0.2)

7.3.4 Approximating with a normal distribution

To see whether we can approximate this sampling distribution with a normal distribution, we need to find the mean and standard deviation. Recall, that the standard deviation of a sampling distribution has a special name: the standard error (SE).

samples_size5 <- nhanes %>%
    rep_slice_sample(n = 5, reps = 1000)

sample_props_mari_size5 <- samples_size5 %>%
    group_by(replicate) %>%
    summarize(stat = mean(Marijuana == "Yes"))

sample_props_mari_size5 %>%
    ggplot(aes(x = stat)) +
    geom_histogram(aes(y = ..density..), binwidth = 0.2) +
    stat_function(fun = dnorm, args = list(mean = ___, sd = ___), color = 'darkred')

7.3.5 Larger samples and the central limit theorem

Recall that the central limit theorem says that the normal distribution will probably be a good model if the samples that go into our sampling distribution are sufficiently large (and if each observation in our sample is independent of the others, which is guaranteed by random sampling). If so, then the mean of the sampling distribution should be close to the population proportion \(\pi\) and the standard error should be close to

\[ SE = \sqrt{\frac{\pi \left(1 - \pi \right)}{n}} \]

Let’s see what happens if we use those when we fit our normal distribution to the sampling distribution. Remember that we found the population proportion earlier (it was about \(0.585\)).

samples_size50 <- nhanes %>%
    rep_sample_n(size = 50, reps = 1000)

sample_props_mari_size50 <- samples_size50 %>%
    group_by(replicate) %>%
    summarize(stat = mean(Marijuana == "Yes"))

sample_props_mari_size50 %>%
    ggplot(aes(x = stat)) +
    geom_histogram(aes(y = ..density..), binwidth = 0.02) +
    stat_function(fun = dnorm, args = list(mean = ___, sd = ___), color = 'darkred')

7.3.6 What if we don’t know the “true” population?

Unlike in the previous example, in practice we usually have a single sample and don’t know what the true value of the population parameter is. But we can still use our sample to learn something about the population in general. We’ve seen how to use bootstrapping to construct a 95% confidence interval for a population proportion based on a sample. Let’s do the same with the normal distribution.

To do this, we will do the following:

1. Draw a random sample of size 50 from the complete NHANES dataset.
2. Use R to calculate the proportion of people in the same who have tried marijuana; this is our sample statistic \(\hat{p}\) or p_hat.
3. Use \(\hat{p}\) from our sample to calculate the standard error according to the same formula as above: \[ SE = \sqrt{\frac{\hat{p} \left(1 - \hat{p} \right)}{n}} \]
4. Make a 95% confidence interval using our knowledge that 95% of a normal distribution is within 1.96 standard deviations of the mean; the 95% CI is \(\left( \hat{p} - 1.96 \times SE, \hat{p} + 1.96 \times SE \right)\).

my_sample <- nhanes %>%
    slice_sample(n = 50)

p_hat <- my_sample %>%
    specify(response = Marijuana, success = "Yes") %>%
    calculate(stat = "prop") %>%
    pull()

standard_error <- sqrt((___ * (1 - ___)) / ___)

c(___ - 1.96 * ___, ___ + 1.96 * ___)