Course textbook for F2026 under active development
đź“„ Download the R code from this chapter
đź“ť Download the class handout for this chapter
This chapter picks up two extensions of the inferential machinery from the last chapter. First, statistical power — the flip-side of \(\alpha\) that lets us reason about false negatives, understand what a given sample size can and cannot detect, and diagnose when a study “found nothing.” Second, the bootstrap — a computational trick that produces a sampling distribution for essentially any statistic we compute in a sample, even ones for which no analytical standard-error formula exists.
Let’s think again about the possible outcomes of a hypothesis test.
The thing we have control over here is \(\alpha\). It is what we explicitly set and therefore it is what we set it at. Keep in mind that 5% is…. pretty high! In my other job I work on the NBC news decision desk where our job is to explicitly reject the null hypothesis that the race is tied for the alternative hypothesis that a candidate has won. In 2022 we called around 600 races. If our standard was a 5% false positive rate we would get \(600*.05=30\) races wrong! We would absolutely lose our jobs. We say that we have an \(\alpha=.005\), but that would still be \(600*.005=3\) wrong calls. That would actually still be unacceptable.
\(\beta\), the probability of a false negative, is a bit more of a mystery. And we have to make some assumptions to get there.
Last time we ran a simulation to show that if we assume the true level of support for Kelly in Arizona is 52.4%, and have a sample of approximately 2600 people, we get a false negative (\(\beta\)) of 32%:
set.seed(19104)
false.negative <- rep(NA,10000)
n <- 2613
alpha <- .05
for(i in 1:10000){
samp <- rbinom(n,1, .524)
s <- sd(samp)
xbar <- mean(samp)
se <- s/sqrt(n)
t.val <- (xbar-.5)/se
#Two-tailed p-value via the absolute t so this works whether xbar happens
#to land above or below .5 (see the abs() explanation in the HypothesisTests
#chapter's CI/hypothesis equivalence loop).
p <- 2*pt(abs(t.val), df=n-1, lower.tail=F)
false.negative[i] <- p>=alpha
}
mean(false.negative)[1] 0.3215Where does that number come from? Can we visualize this given what we know about sampling distributions?
We are positing in the above simulation that the true level of support for Kelly is 52.4%. Given that, here is what the sampling distribution for sample means will look like:
x <- seq(.46,.57, .001)
n <- 2613
mu <- .524
se <- sqrt((mu*(1-mu))/n)
alpha <-.05
plot(x, dnorm(x, mu, se), type="l", col="dodgerblue", main="What is Beta?")
legend("topleft", "True Sampling Distribution", lty=1, col="dodgerblue")
Now, when we get samples in the real world we won’t know that this is the true sampling distribution, and we will perform a null hypothesis test. In most cases for election data like this we will perform the test to determine if the true population mean is 50%, or not.
To perform this hypothesis test, we place the sampling distribution on our null hypothesis, 50%:
no <- .5
plot(x, dnorm(x, mu, se), type="l", col="dodgerblue", main="What is Beta?")
points(x, dnorm(x, no, se), type="l", col="firebrick")
legend("topleft", c("True Sampling Distribution", "Null Hyp. Sampling Distribution"), lty=c(1,1), col=c("dodgerblue", "firebrick"))
In performing this hypothesis test, one of our steps would be to set up the \(\alpha\) region: the area in the tails of the null hypothesis distribution that contain 5% of the probability mass in total. This is the region where, if a mean falls, we will reject the null hypothesis.
plot(x, dnorm(x, mu, se), type="l", col="dodgerblue", main="What is Beta?")
points(x, dnorm(x, no, se), type="l", col="firebrick")
normal.shader(.46, .56, no, se, qnorm(alpha/2, no, se, lower.tail=T ), greater=F, col="firebrick")
normal.shader(.46, .56, no, se, qnorm(alpha/2, no, se, lower.tail=F ), greater=T, col="firebrick")
legend("topleft", c("True Sampling Distribution", "Null Hyp. Sampling Distribution"), lty=c(1,1), col=c("dodgerblue", "firebrick"))
If that’s true, then this is the region where if a mean falls we will fail to reject the null hypothesis:
plot(x, dnorm(x, mu, se), type="l", col="dodgerblue", main="What is Beta?")
points(x, dnorm(x, no, se), type="l", col="firebrick")
normal.shader(.46, .56, no, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="firebrick")
points(x, dnorm(x, mu, se), type="l", col="dodgerblue")
legend("topleft", c("True Sampling Distribution", "Null Hyp. Sampling Distribution"), lty=c(1,1), col=c("dodgerblue", "firebrick"))
So: the blue distribution is what is actually producing poll results (here the proportion of people who support Kelly, with the truth being 52.4% do), while the red region is the range where, if a mean is produced there, we will fail to reject the null hypothesis that the race is tied.
Here is what we are doing to determine \(\beta\): what is the area under the blue curve that corresponds to the shaded red area?
plot(x, dnorm(x, mu, se), type="l", col="dodgerblue", main="What is Beta?")
points(x, dnorm(x, no, se), type="l", col="firebrick")
normal.shader(.46, .56, no, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="firebrick")
points(x, dnorm(x, mu, se), type="l", col="dodgerblue")
legend("topleft", c("True Sampling Distribution", "Null Hyp. Sampling Distribution"), lty=c(1,1), col=c("dodgerblue", "firebrick"))
normal.shader(.46, .56, mu, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="dodgerblue")
What is this region? This blue region are the means that will get produced by the true sampling distribution that we will fail to reject the null hypothesis for. In other words, these are the false negatives!
What proportion of the blue curve is made up of the blue shaded region?
plot(x, dnorm(x, mu, se), type="l", col="dodgerblue", main="What is Beta?")
points(x, dnorm(x, no, se), type="l", col="firebrick")
normal.shader(.46, .56, no, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="firebrick")
points(x, dnorm(x, mu, se), type="l", col="dodgerblue")
legend("topleft", c("True Sampling Distribution", "Null Hyp. Sampling Distribution"), lty=c(1,1), col=c("dodgerblue", "firebrick"))
normal.shader(.46, .56, mu, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="dodgerblue")
#Use the null distribution to determine the rejection bounds
lower.bound <- qnorm(alpha/2, no, se, lower.tail=T )
upper.bound <- qnorm(alpha/2, no, se, lower.tail=F )
#Now use the *true* sampling distribution to determine the probability
pnorm(upper.bound, mu, se) - pnorm(lower.bound, mu, se)[1] 0.309763About 31%! Same as what we found through simulation! (Well, close, there is some rounding happening).
Put another way: 31% of the means that are produced by the true sampling distribution will lead us to fail to reject the null hypothesis. Our false-negative rate (\(\beta\)) is 31%.
What can we do to reduce the false negativity rate?
Looking at our charts, the first (and stupidest) thing we can do is to have the truth be further away from the null hypothesis:
n <- 2613
mu <- .53
se <- sqrt((mu*(1-mu))/n)
no <- .5
alpha <- .05
plot(x, dnorm(x, mu, se), type="l", col="dodgerblue", main="What is Beta?")
points(x, dnorm(x, no, se), type="l", col="firebrick")
normal.shader(.46, .56, no, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="firebrick")
points(x, dnorm(x, mu, se), type="l", col="dodgerblue")
legend("topleft", c("True Sampling Distribution", "Null Hyp. Sampling Distribution"), lty=c(1,1), col=c("dodgerblue", "firebrick"))
normal.shader(.46, .56, mu, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="dodgerblue")
#Use the null distribution to determine the rejection bounds
lower.bound <- qnorm(alpha/2, no, se, lower.tail=T )
upper.bound <- qnorm(alpha/2, no, se, lower.tail=F )
#Now use the *true* sampling distribution to determine the probability
pnorm(upper.bound, mu, se) - pnorm(lower.bound, mu, se)[1] 0.1329351The two distributions are the same shape as above, but now the “truth” is that Kelly is going to win 53% of the vote. Because this is simply further away from the null then we are much less likely to get a false negative (the blue region is much smaller). Taken to the extreme: if everyone was going to vote for Kelly there obviously would be no chance of getting a false negative! This is “stupid” because we can’t control what the truth is!
Given that the truth is fixed, we can also think about changing other parameters to lower our false positive rate.
The first thing that we can do is to raise \(\alpha\) the false positivity rate. If we are ok with getting more false positives, we will have less false negatives:
n <- 2613
mu <- .524
se <- sqrt((mu*(1-mu))/n)
no <- .5
alpha <- .1
plot(x, dnorm(x, mu, se), type="l", col="dodgerblue", main="What is Beta?")
points(x, dnorm(x, no, se), type="l", col="firebrick")
normal.shader(.46, .56, no, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="firebrick")
points(x, dnorm(x, mu, se), type="l", col="dodgerblue")
legend("topleft", c("True Sampling Distribution", "Null Hyp. Sampling Distribution"), lty=c(1,1), col=c("dodgerblue", "firebrick"))
normal.shader(.46, .56, mu, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="dodgerblue")
#Use the null distribution to determine the rejection bounds
lower.bound <- qnorm(alpha/2, no, se, lower.tail=T )
upper.bound <- qnorm(alpha/2, no, se, lower.tail=F )
#Now use the *true* sampling distribution to determine the probability
pnorm(upper.bound, mu, se) - pnorm(lower.bound, mu, se)[1] 0.2084848The area between the two lines is now somewhat smaller, and as such the probability that a mean from the red area is produced in that zone is smaller. We see that this drops the false negative rate to 20%. (Generally speaking we don’t want to make this trade-off! False positives are way worse than false negatives. We would much rather say a drug doesn’t work when it does, than say a drug does work when it doesn’t.)
The other thing that we can do is to reduce the width of the sampling distributions. Then, for a given level of alpha and a given difference between the truth and the null hypothesis there will be less overlap between the two distributions.
The easiest way to do this is to increase the sample size. Let’s increase the sample size to 5000 people:
n <- 5000
mu <- .524
se <- sqrt((mu*(1-mu))/n)
no <- .5
alpha <- .05
plot(x, dnorm(x, mu, se), type="l", col="dodgerblue", main="What is Beta?")
points(x, dnorm(x, no, se), type="l", col="firebrick")
normal.shader(.46, .56, no, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="firebrick")
points(x, dnorm(x, mu, se), type="l", col="dodgerblue")
legend("topleft", c("True Sampling Distribution", "Null Hyp. Sampling Distribution"), lty=c(1,1), col=c("dodgerblue", "firebrick"))
normal.shader(.46, .56, mu, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="dodgerblue")
#Use the null distribution to determine the rejection bounds
lower.bound <- qnorm(alpha/2, no, se, lower.tail=T )
upper.bound <- qnorm(alpha/2, no, se, lower.tail=F )
#Now use the *true* sampling distribution to determine the probability
pnorm(upper.bound, mu, se) - pnorm(lower.bound, mu, se)[1] 0.07520771Or to 10000 people:
n <- 10000
mu <- .524
se <- sqrt((mu*(1-mu))/n)
no <- .5
alpha <- .05
plot(x, dnorm(x, mu, se), type="l", col="dodgerblue", main="What is Beta?")
points(x, dnorm(x, no, se), type="l", col="firebrick")
normal.shader(.46, .56, no, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="firebrick")
points(x, dnorm(x, mu, se), type="l", col="dodgerblue")
legend("topleft", c("True Sampling Distribution", "Null Hyp. Sampling Distribution"), lty=c(1,1), col=c("dodgerblue", "firebrick"))
normal.shader(.46, .56, mu, se, qnorm(alpha/2, no, se, lower.tail=T ),qnorm(alpha/2, no, se, lower.tail=F ), between=T, col="dodgerblue")
#Use the null distribution to determine the rejection bounds
lower.bound <- qnorm(alpha/2, no, se, lower.tail=T )
upper.bound <- qnorm(alpha/2, no, se, lower.tail=F )
#Now use the *true* sampling distribution to determine the probability
pnorm(upper.bound, mu, se) - pnorm(lower.bound, mu, se)[1] 0.002216565The false negativity rate drops off substantially.
Given that the width of the sampling distribution is given by \(s/\sqrt{n}\) clearly we can either increase \(n\) or decrease \(s\) to make the sampling distribution smaller. In the case where we are sampling a bernoulli it’s very difficult to make \(s\) smaller, but it might be in other cases. For example a common survey question we ask are “Feeling Thermometers” where we ask people to rate how warmly they feel about a politician/group where 0 is cold and 100 is warm. As it turns out people are really bad at answering these questions and they have a lot of random variation. All things being equal, that high variance question is going to lead to a higher \(\beta\).
When we are performing research we often don’t think in terms of \(\beta\), but instead consider \(1-\beta\), which we call “statistical power”. If \(\beta\) is the probability of failing to reject the null hypothesis when we should, then power is the probability of rejecting the null hypothesis when we should. Put in other words, it’s our ability to say that something is there when something is in fact there.
Here’s an absurd example: Biden’s highest approval rating was 55%. Let’s say a survey firm completes a random survey of 15 people, generates an approval rating for Biden, but claims that they can’t conclude whether or not Biden is more popular than Trump was when he left office with a approval of 38%.
Let’s use the technique above to map out the true and null hypotheses for this firms test:
#Input
n <- 15
hA <- .55
hO <- .38
alpha <- .05
#SE of alternative
sd <- sqrt(hA*(1-hA))
seA <- sd/sqrt(n)
#SE of Null
sd <- sqrt(hO*(1-hO))
seO <- sd/sqrt(n)
#Calculate rejection zone for alpha
z.score <- qnorm(alpha/2)*-1
bounds <- c(hO+z.score*seO, hO-z.score*seO)
eval <- seq(.1,.7,.0001)
plot(eval, dnorm(eval, mean=hA, sd=seA), type="l", col="darkblue", main="What is Beta?")
legend("topleft","True Sampling Distribution", lty=1, col="darkblue")
points(eval, dnorm(eval, mean=hO, sd=seO), type="l", col="firebrick")
normal.shader(.1,.7, hO, seO, bounds[2], bounds[1], between=T, col="firebrick")
normal.shader(.1,.7, hA, seA, bounds[2], bounds[1], between=T, col="darkblue")
legend("topleft",c("True Sampling Distribution","Null Hyp. Sampling Distribution") , lty=c(1,1), col=c("darkblue", "firebrick"))
#Beta calculation:
pnorm(bounds[1], mean=hA, sd=seA) - pnorm(bounds[2], mean=hA, sd=seA) [1] 0.7214013With 15 people we have an enormously wide zone where we will fail to reject the null hypothesis. A wide band that is necessary to ensure our 5% false positive rate. We need extraordinary evidence to claim that a sample mean is different from 38% when we only have 15 people. Accordingly, we get a huge number of false negatives! Indeed, it’s more likely than not we will get a false negative.
Next week when we get to difference in means tests we will similarly be able to calculate the minimum detectable effect size given the \(n\) of an experiment, for example. An extraordinarily common critique of experiments is that authors failed to run a power analysis and, given the effect size that is plausible and their \(n\), they simply lacked the statistical power to find anything!
Similarly common in experiments is to what we saw above: authors may say: we ran an experiment and found nothing, so we conclude that nothing is there. But a lot of times power analysis will reveal that a 20 or 30% false negativity rate is expected. You can’t claim that nothing is happening when you get a result of “nothing” 30% of the time!
Can you have too much power? Not really, but the results of too much power can be somewhat misleading. Let’s say the true proportion of people who support Joe Biden is 38.5%, and a firm runs a poll with 1 million people in it. They get a result and conclude that Joe Biden’s support is significantly different than Donald Trump’s low point of 38%. What would that look like:
#Input
n <- 1000000
hA <- .385
hO <- .38
alpha <- .05
#SE of alternative
sd <- sqrt(hA*(1-hA))
seA <- sd/sqrt(n)
#SE of Null
sd <- sqrt(hO*(1-hO))
seO <- sd/sqrt(n)
#Calculate rejection zone for alpha
z.score <- qnorm(alpha/2)*-1
bounds <- c(hO+z.score*seO, hO-z.score*seO)
eval <- seq(.34,.44,.0001)
plot(eval, dnorm(eval, mean=hA, sd=seA), type="l", col="darkblue", main="What is Beta?")
legend("topleft","True Sampling Distribution", lty=1, col="darkblue")
points(eval, dnorm(eval, mean=hO, sd=seO), type="l", col="firebrick")
abline(v=bounds, lty=2, col="firebrick")
legend("topleft",c("True Sampling Distribution","Null Hyp. Sampling Distribution") , lty=c(1,1), col=c("darkblue", "firebrick"))
#Beta calculation:
pnorm(bounds[1], mean=hA, sd=seA) - pnorm(bounds[2], mean=hA, sd=seA) [1] 4.383925e-17It’s extraordinarily unlikely to get a false negative when you have an \(n\) that big, and as such very small differences from the null can be said to be “statistically significant”. But here is where we get into the art rather than the science of statistics, just because something is statistically significant doesn’t mean that it is substantively significant. I don’t think that anyone would claim that an approval rating of 38.5% is substantively different than an approval rating of 38%. You can get a very large sample and conclude that something is statistically significant, but that says nothing about the relative magnitude or importance of that effect.
My other favorite example of this is meat cancer. And by favorite I mean a disaster of science communication. Whole industries were launched, careers made, and documentaries written by the fact that the WHO labeled red meat a “Group 1 Carcinogen”, the same group that “arsenic, asbestos, and tobacco belong to”.
What was happening here? Is red meat as damaging to you as smoking? No! Not even close.
What happened was that there was a large accumulation of studies that, in total, showed that red meat had a small effect on the rates of colon cancer. The WHO’s “Group 1” is not a measure of how severe a carcinogen is, but is a measure of statistical significance. Enough studies had been done that there was sufficient evidence to overturn the null hypothesis that meat did not effect the rate of cancer.
Because people generally don’t think like this, people assumed what the WHO meant was that bacon was just as harmful as smoking.
The meat cancer scare was helped along by some truly weird choices by the epidemiologists about how to communicate the risk. All the headlines stated that “Eating 50g of processed meat a day increases the risk of colon cancer by 18%” A reasonable person would think that, if the baseline rate of colon cancer is 3%, then eating 50g of processed meat a day would lead to a rate of colon cancer of 3+18=21%. But that’s not at all what they were saying, they were saying that if the baseline rate of colon cancer was 3%, then eating 50g of processed meat a day would lead to a race of colon cancer of 3*1.18=3.54%. The 18% was a “relative risk” calculation, not a change in percentage points. That’s a lot less dramatic!
Let’s compare that to another “Group 1” carcinogen. The probability of lung cancer for never smokers is approximately .015%, that is, approximately 1 in 1000. The probability of lung cancer for smokers is approximately 25%, that is, 1 in 4. As such, compared to the “18%” relative risk of colon cancer from meat, there is a 1666% relative risk increase of lung cancer from smoking.
So if anyone tells you that eating meat is as dangerous as smoking, tell them Dr. Trussler said that they are a clown.
That being said we should probably all cut back on our red and processed meat consumption :).
If everyone misinterprets what you said as a scientist, it’s not their fault, it’s your fault.
Correlation, which we develop fully in the next chapter, is a statistic that measures how strongly two variables move together. R computes it with cor(x, y) and it returns a single number between -1 and 1. We use correlation here as our motivating example because it’s a statistic with no known analytical formula for the standard error — and that’s exactly the situation the bootstrap is built for.
There is a command for it in R, cor(), that we can use to see what it produces:
set.seed(19103)
sigma<-rbind(c(1,.2), c(.2,1))
d <- as.data.frame(mvrnorm(n=500, mu=c(0,0), Sigma=sigma))
names(d) <- c("x","y")
plot(d$x, d$y )
correlation <- cor(d$x, d$y)
correlation[1] 0.2569832One thing to note about the correlation function is that it won’t give you anything back if there missing values in either of the vectors:
z <- d$x
z[5] <- NA
cor(z, d$y)[1] NAOur general answer to how to deal with these things is to use na.rm=T, but that’s not what we do with a correlation. Because there are a few different ways that NA could be dealt with for a correaltion, here we need to specify the method which we want R to deal with NAs. The safe thing to use is:
cor(z, d$y, use="pairwise.complete")[1] 0.2596632That’s really great, but given what we have learned so far in this class about samples and population, what else would we like to know about this correlation coefficient?
Well we know that this sample of data is only one of an infinite number of samples we could get. We know that there is one “true” correlation between these two variables in the world, and we would like to make an inference about what that might be.
As with everything: we want to perform a hypothesis test, so we want a sampling distribution.
So let’s dig in to to the correlation function to find a standard error:
names(correlation)NULLUh-oh! It doesn’t give us one!
And if we go to to the textbook there is no mention of a standard error or sampling distribution for a correlation coefficient!
Well then what are we going to do??
I want to use the rest of class to introduce a method that can give you a sampling distribution for anything you calculate, called the bootstrap.
Why is it called the bootstrap? The phrase “picking yourself up by your bootstraps” specifically refers to doing something impossible. You can’t actually do it. Similarly, the bootstrap statistical method allows you to generate a sampling distribution for anything based off of a single sample, even in cases where there is no known equation for the standard error.
Before we get to the method (which is actually pretty easy), let’s generate the true sampling distribution for this correlation coefficient. In this case we can do so because we generated these data. So let’s repeatedly re-generate these data and calculate the correlation a large number of times so that we can determine the true sampling distribution:
cor.samp.dist <- rep(NA, 10000)
for(i in 1:10000){
sigma<-rbind(c(1,.2), c(.2,1))
samp <- as.data.frame(mvrnorm(n=500, mu=c(0,0), Sigma=sigma))
names(samp) <- c("x","y")
cor.samp.dist[i] <- cor(samp$x, samp$y)
}
plot(density(cor.samp.dist))
abline(v=0, lty=2)
sd(cor.samp.dist)[1] 0.0428987But pretend that all we have is our original dataset, d. Our goal is to recreate the actual sampling distribution using information wholly contained within our one sample.
Here is the key insite to the bootstrap: because our sample of data is assumed to be an iid sample of the population, a random sample of our sample is a new sample.
To put it in more reasonable terms: we have seen already that any given sample is a shadow, or is suggested by, the population. That’s true for the distribution of one variable, but is also true for the relationships between variables (or anything else). What we want to do is to use the fact that our sample is representative of the population to generate new samples to explore the range of possibilities that might occur if we truly re-sampled a large number of times.
The way that we are going to actually do this is to repeatedly generate new datasets of the same length as our original by sampling with replacement from our original dataset. Here is one such new sample:
bs.samp <- d[sample(1:nrow(d), replace=T),]
cor(d$x, d$y)[1] 0.2569832cor(bs.samp$x, bs.samp$y)[1] 0.2310803#Not the same!Why is it important that we sample with replacement. Well if we sample without replacement we will just select the rows of the old dataset one at a time until we have selected all of them. In the end we will be left with an exact re-creation of our old dataset in a new order:
bs.samp <- d[sample(1:nrow(d), replace=F),]
cor(d$x, d$y)[1] 0.2569832cor(bs.samp$x, bs.samp$y)[1] 0.2569832#The same!Having repeated entries is key into creating variance that generates different samples everytime you generate a new bootstrap sample.
OK, that’s all fine, but does it work? Let’s repeatedly generate new bootstrap samples and calculate the correlation coefficient.
bs.samp.dist <- rep(NA, 10000)
for(i in 1:10000){
bs.samp <- d[sample(1:nrow(d), replace=T),]
bs.samp.dist[i] <- cor(bs.samp$x, bs.samp$y)
}And let’s compare to the true sampling distribution:
plot(density(cor.samp.dist), col="darkblue")
points(density(bs.samp.dist), col="firebrick", type="l")
abline(v=0, lty=2)
sd(cor.samp.dist)[1] 0.0428987sd(bs.samp.dist)[1] 0.04045926Now, these two sampling distributions are not centered in the same place. But does that really matter? The true sampling distribution is centered around the true population correlation, which we are never going to know. The bootstrap sampling distribution is centered around the correlation in the sample.
But think about when we calculate a standard error in a sample using an equation. In that case all we would get is a standard error, we don’t get to recover the true location of a mean, or difference in mean, or correlation.
What we do get is an extremely reasonable approximation of the standard deviation of the sampling distribution. And what is that: the standard error!
What can we do with our bootstrap estimates?
Usually what people do is to use the bootstrap estimates to form a 95% confidence interval. We can determine this empirically by using the quantile function
quantile(bs.samp.dist, .025) 2.5%
0.1764792 quantile(bs.samp.dist, .975) 97.5%
0.3348584 If we are willing to assume that the sampling distribution is normally distributed, we can also use the standard error to perform a hypothesis test under a null hypothesis:
z.score <- cor(d$x, d$y)/sd(bs.samp.dist)
p <- pnorm(z.score, lower.tail=F)*2
p[1] 2.130112e-10Let’s prove to ourselves that the bootstrap really is doing what we claim by running it on a case where we already know the true sampling distribution.
Let’s generate one sample of data from a known normal distribution:
set.seed(19104)
prime.sample <- rnorm(1000, mean=0, sd=5)
mean(prime.sample)[1] 0.1236452sd(prime.sample)[1] 5.035127From this one sample we can estimate a sampling distribution of this sample mean. To do so we need to calculate a standard error:
se.calc <- sd(prime.sample)/sqrt(1000)Now this sampling distribution is techinically \(t\) distributed, but we have 1000 observations here so the t distribution will have converged on the normal distribution.
So here is this calculated sampling distribution, centered on 0 (though it doesn’t really matter where we center it):
eval <- seq(-1,1, .0001)
plot(eval, dnorm(eval, mean=0, sd=se.calc), col="firebrick", type="l", main="Sampling Distributions")
Now, we know the population in this case, so we can also use that information to calculate the true standard error and sampling distribution:
se.real <- 5/sqrt(1000)
plot(eval, dnorm(eval, mean=0, sd=se.calc), col="firebrick", type="l", main="Sampling Distributions")
points(eval, dnorm(eval, mean=0, sd=se.real), col="darkblue", type="l")
legend("topleft", c("Calculated SampDist", "Real SampDist"), lty=c(1,1), col=c("firebrick","darkblue"))
As we have seen they are approximately the exact same shape.
Now, again, because we know the true population we could also repeatedly sample from this population, calculating a mean in each sample, and that should also give us exactly the same sampling distribution.
samp.dist <- rep(NA, 10000)
for(i in 1:10000){
samp.dist[i] <- mean(rnorm(1000, mean=0, sd=5))
}
plot(eval, dnorm(eval, mean=0, sd=se.calc), col="firebrick", type="l", main="Sampling Distributions")
points(eval, dnorm(eval, mean=0, sd=se.real), col="darkblue", type="l")
points(density(samp.dist), col="forestgreen", type="l")
legend("topleft", c("Calculated SampDist", "Real SampDist", "Empirical SampDist"), lty=c(1,1,1), col=c("firebrick","darkblue", "forestgreen"))
This sampling distribution is also the same.
Finally, what I showed last time is that there is an alternative method to determining the sampling distribution: the bootstrap. To produce the green sampling distribution we repeatedly sampled from our population slightly different sample of 1000, calculating the mean in each, and the resulting distribution is the sampling distribution. What the magic of the bootstrap is that we can get similarly “new samples” by sampling our original data with replacement.
For example here are 9 samples of data compared against the population:
par(mfrow=c(3,3))
eval.pop <- seq(-20,20,.01)
for(i in 1:9){
samp <- rnorm(1000, mean=0, sd=5)
plot(eval.pop, dnorm(eval.pop, mean=0, sd=5), lwd=2, lty=1, type="l",
main=paste("Mean = ", round(mean(samp),2)))
points(density(samp), col="forestgreen", type="l")
}
They are all approximately the same, but with a little bit of variation. That “little bit of variation” causes the variation in the mean that generates a sampling distribution.
Again, the claim of the bootstrap is that we can do exactly the same thing by re-sampling from our existing data with replacement:
par(mfrow=c(3,3))
eval.pop <- seq(-20,20,.01)
for(i in 1:9){
samp <- prime.sample[sample(1:length(prime.sample), length(prime.sample), replace=T)]
plot(eval.pop, dnorm(eval.pop, mean=0, sd=5), lwd=2, lty=1, type="l",
main=paste("Mean = ", round(mean(samp),2)))
points(density(samp), col="magenta", type="l")
}
VERY SIMILIARLY, when we re-sample with replacement from our prime sample we get new samples.
If we do this a large number of times (10,000) we will get a new sampling distribution:
bs.samp.dist <- rep(NA,10000)
for(i in 1:10000){
bs.samp.dist[i] <- mean(prime.sample[sample(1:length(prime.sample), length(prime.sample), replace=T)])
}To compare to the other sampling distributions I’m also going to center this distribution on 0, by subtracting the mean from all values (this is only necessary for visualization, this isn’t a usual step in the bootstrap process):
bs.samp.dist <- bs.samp.dist - mean(bs.samp.dist)plot(eval, dnorm(eval, mean=0, sd=se.calc), col="firebrick", type="l", main="Sampling Distributions")
points(eval, dnorm(eval, mean=0, sd=se.real), col="darkblue", type="l")
points(density(samp.dist), col="forestgreen", type="l")
points(density(bs.samp.dist), col="magenta", type="l")
legend("topleft", c("Calculated SampDist", "Real SampDist", "Empirical SampDist", "BS SampDist"), lty=c(1,1,1,1), col=c("firebrick", "darkblue","forestgreen","magenta"))
Again, it’s the same!
So 4 ways to calculate a sampling distribution:
So in the real world the options for calculating a standard error/sampling distribution are (3) and (4). If (3) is available to us we always prefer it, but if not, option (4) is available.
Let’s consider a real world example.
Here is some data from the American National Election Study, which is the premier academic study of elections.
library(rio)
anes <- import("https://github.com/marctrussler/IIS-Data/raw/main/ANESFinalProjectData.csv")I’m going to look at how people voted in the democratic primary, classifying candidates into centrist and progressives:
#5 and 6 are Sanders and Warren
anes$vote.centrist.dem[anes$V201021 %in% c(1,2,3,4)] <- 1
anes$vote.centrist.dem[anes$V201021 %in% c(5,6)] <- 0
table(anes$vote.centrist.dem)
0 1
731 1695 And I want to know if ideological self placement predicts that:
#1 is very liberal and 7 is very conservative.
anes$ideology <- anes$V201200
anes$ideology[anes$ideology %in% c(-9,-8,99)] <- NA
table(anes$ideology)
1 2 3 4 5 6 7
369 1210 918 1818 821 1492 428 Here is a cross-table of these two things:
table(anes$vote.centrist.dem, anes$ideology)
1 2 3 4 5 6 7
0 165 282 121 89 13 10 4
1 81 492 334 449 74 37 5Even though these are not wide-ranging variables, a correlation can still tell us the degree to which they co-vary:
cor(anes$ideology, anes$vote.centrist.dem, use="pairwise.complete")[1] 0.2734469There is a small correlation between self reported ideology and voting for a centrist candidate in the democratic primary.
Based on recent research by Hakeem Jefferson at Stanford, it might be that the relationship between ideology and vote choice among Black Americans may be different than White Americans. Specifically, Jefferson has found that the cultural context around “liberal” and “conservative” are extraordinarily different for Black Americans such that they are effectively answering a different question, one that (predominantly white) researchers think of completely different.
So let’s look at the correlation among white and among Black Americans sperately.
anes$white.v.black[anes$V201549x==1] <- 1
anes$white.v.black[anes$V201549x==2] <- 0
table(anes$white.v.black)
0 1
726 5963 Let’s see if the correlation is different:
#Among White Americans
cor(anes$ideology[anes$white.v.black==1], anes$vote.centrist.dem[anes$white.v.black==1], use="pairwise.complete")[1] 0.2818327#Among Black Americans
cor(anes$ideology[anes$white.v.black==0], anes$vote.centrist.dem[anes$white.v.black==0], use="pairwise.complete")[1] 0.1446376#Difference
cor.delta <- cor(anes$ideology[anes$white.v.black==1], anes$vote.centrist.dem[anes$white.v.black==1], use="pairwise.complete") - cor(anes$ideology[anes$white.v.black==0], anes$vote.centrist.dem[anes$white.v.black==0], use="pairwise.complete")
cor.delta[1] 0.1371952Yes! It is substantially lower!
Ok but…. We know that every time we take a new sample we are going to get a slightly different answer. And while the ANES overall is pretty big it only has 700 or so Black Americans, not all of whom voted in the democratic primary. So there theoretically would be a lot of variation in that difference in correlations of .14 if we sampled a large number of times. So we want to perform a hypothesis test of some sort to determine if we can reject the null hypothesis that these two correlations are actually equal to one another in the population.
Looking at our 4 ways of calculating standard errors/sampling distributions above, what can we do to calculate a sampling distribution. Well this is real data so we don’t know the population, which rules out 1 and 2. Do we have a formula to calculate the standard error on the difference between two correlations? No! We don’t even have a formula for one correlation, let alone the difference.
So we have no choice but to use the bootstrap.
What we want to do is to sample rows of our data set with replacement to create “new” datasets. In each we will calculate the difference in the correlations between white and black americans.
To substantially reduce the computing power needed for this, i’m going to reduce the dataset down to just the three variables we need, and only keep observations with non-NA data for all three variables.
boot.data <- anes[c("ideology","vote.centrist.dem","white.v.black")]
boot.data <- boot.data[complete.cases(boot.data),]The bootstrap procedure:
bs.cor.delta <- rep(NA,10000)
for(i in 1:10000){
#Generate the bootstrap dataset
bs.data <- boot.data[sample(1:nrow(boot.data), nrow(boot.data), replace=T),]
#Calculat the same estimate as I calculated in the prime dataset, saving the result each time.
bs.cor.delta[i] <- cor(bs.data$ideology[bs.data$white.v.black==1],bs.data$vote.centrist.dem[bs.data$white.v.black==1],use="pairwise.complete") - cor(bs.data$ideology[bs.data$white.v.black==0],bs.data$vote.centrist.dem[bs.data$white.v.black==0], use="pairwise.complete")
}plot(density(bs.cor.delta), main="Bootstrap Sampling Distribution")
abline(v=0, lty=2)
Ok! That gave us a sampling distribution for the difference between these two correlations. What can we do with this?
We can use this empirical distribution to generate a confidence interval. We can literally look at the two values which define the range that contains 95% of bootstrap estimates.
#Our original estimate
cor.delta[1] 0.1371952#CI
quantile(bs.cor.delta, .025) 2.5%
0.01662887 quantile(bs.cor.delta, .975) 97.5%
0.2654626 Given that the confidence interval does not overlap 0 we can conclude that there is a lower than 5% probability of seeing a difference in correlation this extreme if the truth was that the two correlations were equal.
We can also use the standard deviation of these bootstrap samples as a standard error:
se <- sd(bs.cor.delta)
#For my own sanity, make sure that the sample distribution is normally distributed:
eval <- seq(-.1, .5, .0001)
plot(density(bs.cor.delta), main="Bootstrap Sampling Distribution")
points(eval, dnorm(eval, mean=cor.delta, sd=se), col="firebrick", type="l")
abline(v=0, lty=2)
#For null of 0, how many SEs is our test statistic away?
z.score <- (cor.delta-0)/se
#Evalualte under the standard normal
pnorm(z.score, lower.tail=F)*2[1] 0.03125029There is approximately a 3% chance of seeing something as extreme as a difference of .137 if the true difference in correlations was 0.
Cool! That’s genuinely a new thing I didn’t know that we discovered using these data and the bootstrap!
Power gives us the vocabulary to reason about what a hypothesis test can and cannot detect at a given sample size. The bootstrap gives us a computational way to build a sampling distribution around any statistic we compute in a sample. Both tools generalize the inferential machinery we developed for the sample mean to a much wider set of quantities. Next chapter we introduce covariance and correlation proper — the statistics we’ve been using here as a motivating example — and then in the chapter after that we finally get to regression, the single most useful tool in the applied data scientist’s kit.