gss <- readRDS(here::here("data", "gss2024.rds")) |>
haven::zap_labels()
tv <- gss$tvhours[!is.na(gss$tvhours)]6 Hypothesis testing
We’re going to look at testing null hypotheses. The confidence interval of Chapter 5 is closely linked to the idea of the hypothesis test, which uses the sample data to test if there is enough evidence to assert that the population differs from some specific reference value. That reference value is called the null hypothesis.
6.1 Setting up the null model
Null hypothesis: the average American adult in 2024 watched 3 hours of TV per day. By subtracting 3 from the observed value, we can make it so the null hypothesis is 0. That makes it easy to use lm(y ~ 0) to set up the null model that tests whether the mean of tvhours is different from 0.
In the language of Chapter 4, the null model is model C. Model A estimates the mean.
d <- data.frame(tvdev = tv - 3)
m_c <- lm(tvdev ~ 0, data = d) # assume the mean is 3
m_a <- lm(tvdev ~ 1, data = d) # allow the mean to be different from 3
observed_pre <- (deviance(m_c) - deviance(m_a)) / deviance(m_c)
observed_pre[1] 0.008343581The SSE is improved by estimating \(\beta_0\), but was it improved more than we’d expect by chance?
6.2 Creating a null distribution
Here’s a quick example where we make data where the null hypothesis is actually true. Then we can use it to check whether we could get numbers as large as we see in the real data by chance variation alone. We can use rnorm() to make fake data with a given mean and SD.1
1 We could create a better null distribution here by using a different distribution than the normal for our fake tvhours data. As you saw earlier the data itself is actually NOT normally distributed. We could use a count distribution but we’re not ready for that! We’ll talk about that more in Chapter 17.
Convert this idea to a function so we can do this many times. The basic idea is to simulate data where the null is true, then calculate the PRE (which won’t be exactly zero because of sampling variability). Then we create a skeleton and append 1000 null PRE simulations.
set.seed(522)
n_tv <- length(tv)
sd_tv <- sd(tv)
calc_null_pre <- function() {
null_data <- tibble(tvdev = rnorm(n_tv, 0, sd_tv))
m_c <- lm(tvdev ~ 0, data = null_data)
m_a <- lm(tvdev ~ 1, data = null_data)
pre <- (deviance(m_c) - deviance(m_a)) / deviance(m_c)
return(pre) # this returns the observed PRE from that simulation
}
null_sims <- tibble(sim_number = 1:1000) |>
rowwise() |>
mutate(pre = calc_null_pre())6.3 The test
The idea is to compare the real world to the world implied by the null hypothesis. So we’ll show how the observed PRE compares to the distribution of PREs we got from fake data where the null hypothesis is true.
plt(~ pre, data = null_sims, type = "hist", breaks = 40,
xlab = "PRE when the null is true", ylab = "Simulations")
abline(v = observed_pre, lwd = 2)
ggplot(null_sims, aes(x = pre)) +
geom_histogram(bins = 40) +
geom_vline(xintercept = observed_pre, linewidth = 1) +
labs(x = "PRE when the null is true", y = "Simulations")
The solid line shows that the observed PRE is way higher than we’d expect to get by chance. This gives you the intuition for what the F-test is trying to accomplish in an “analytic” way (i.e., by math alone, not by simulation).
What F would that be equivalent to?
fstat <- (observed_pre / 1) / ((1 - observed_pre) / (n_tv - 1))
fstat[1] 18.09805We can compare this F-statistic we calculated to an F-distribution with numerator df = 1 and denominator df = 2151.
Data prep
df1 <- 1
df2 <- n_tv - 1
f_crit <- qf(.95, df1, df2)
x_max <- ceiling(fstat) + 2
# grid for the curve
x <- seq(0, x_max, length.out = 2000)
y <- df(x, df1 = df1, df2 = df2)Show code
plt(y ~ x,
type = "l",
main = paste0("F(1, ", df2, ")"),
xlab = "F",
ylab = "Density",
xlim = c(0, x_max),
lwd = 2)
abline(v = f_crit, lty = 3, lwd = 2)
abline(v = fstat, col = tableau10[2], lwd = 2)
Show code
ggplot(data.frame(x = x, y = y), aes(x = x, y = y)) +
geom_line(linewidth = 1, color = tableau10[1]) +
geom_vline(xintercept = f_crit, linetype = "dotted", linewidth = 1) +
geom_vline(xintercept = fstat, color = tableau10[2], linewidth = 1) +
coord_cartesian(xlim = c(0, x_max)) +
labs(title = paste0("F(1, ", df2, ")"), x = "F", y = "Density")
The dotted line is the “critical value”—that is, the 95th percentile F-statistic we’d expect to get by chance if the null hypothesis is true. The solid line is the observed value from the data. Again, we can see that it’s way higher than what we’d expect by chance. The formula-based approach and the simulation have the same basic logic and the same conclusion.
Before computers, people looked up F* (the critical value for alpha = .05) in a table in the back of a statistics book. Now we can do it using functions in R (as we did in the code for the graph above).
qf(.95, 1, n_tv - 1)[1] 3.845786Remember that there is nothing sacred about 95% or 99% or any of that.
You can also get the p-value from functions rather than from a table.
1 - pf(fstat, 1, n_tv - 1)[1] 2.188091e-05In any case, we will be rejecting the null hypothesis here!
NoteA limit of simulation
The share of our null simulations with a PRE at least as big as the observed one is 0. But with 1,000 simulations, we can’t see anything much smaller than 1/1000. The p-value from the F distribution is way smaller than that. This is one practical reason it’s nice to have formulas!
6.4 The z-score route
We can approach this issue more generally through z-scores, just like in Chapter 5. How “weird” is our result? How many standard errors is it away from the expected value of the null distribution?
se_tv <- sd(tv) / sqrt(n_tv)
z_tv <- (mean(tv) - 3) / se_tv
c(se = se_tv, z = z_tv, p_value = 2 * (1 - pnorm(abs(z_tv)))) se z p_value
7.110872e-02 4.254180e+00 2.098166e-05 That probability is called the p-value of the test.
R’s t.test() will do the whole thing at once. Since we estimated \(\sigma\) from the data, the right distribution to use is actually the \(t\) distribution rather than the normal (more on that in Chapter 7). But with this many people, the difference is tiny.
t.test(tv, mu = 3)
One Sample t-test
data: tv
t = 4.2542, df = 2151, p-value = 2.188e-05
alternative hypothesis: true mean is not equal to 3
95 percent confidence interval:
3.163060 3.441958
sample estimates:
mean of x
3.302509 6.5 Tails and tests
Why are we using the left and right sides? Why not just use the right side and halve the p-value? After all, it is true that we’d only expect to get a value as large as ours that much of the time.
The use of two-tailed tests rather than one-tailed tests is ubiquitous in sociology. It’s regarded as “conservative” even though there is usually not a good rationale for it.
6.6 Alpha level
How do we connect these ideas to a hypothesis test? To conduct a hypothesis test, we need an alpha level (or \(\alpha\) level). This is the proportion of the time we’re willing to falsely assert that the observed data did not come from the null distribution. This is connected to the idea of type-I error or the idea of a false positive.
We’re now ready for the algorithm of the hypothesis test:
- Choose an alpha level (say, .05)
- Calculate the observed sample statistic
- Calculate the absolute difference between the statistic and the expected value under the null
- Convert this difference into a z-score using the SE of the sampling distribution
- Convert the z-score to a p-value
- If the p-value is less than alpha reject the null hypothesis; if the p-value is greater than alpha fail to reject the null hypothesis
Sometimes people write a hypothesis test out formally. Here’s an example.
\[H_0 : \mu = 3 \qquad H_1 : \mu \neq 3\]
6.6.1 Rejecting (or not) the null
This language can feel weird. We can never accept the null hypothesis (\(H_0\)), in part because the probability of an exact value (e.g., 3) being true is basically zero. So we can only either reject the null hypothesis or fail to reject it.
When we reject the null hypothesis, we call a result statistically significant. That’s literally all that phrase means!
6.6.2 Conventional alpha levels
The conventional alpha level for a test is .05. Heuristically speaking, this means we’re willing to falsely reject the null hypothesis 5% of the time. This value is by no means sacred. In fact, it is fundamentally arbitrary.
Just as we saw with confidence intervals, we can pick any value we like, which is both liberating and scary!
Here’s an example. Among widowed respondents, is a majority afraid to walk alone at night in their own neighborhood?
wid <- gss |>
filter(marital == 2, !is.na(fear)) |>
mutate(afraid = as.numeric(fear == 1))
p_wid <- mean(wid$afraid)
n_wid <- nrow(wid)
z_wid <- (p_wid - 0.5) / sqrt(0.25 / n_wid)
c(n = n_wid, p_hat = p_wid, z = z_wid,
p_value = 2 * (1 - pnorm(abs(z_wid)))) n p_hat z p_value
176.00000000 0.41477273 -2.26133508 0.02373852 At \(\alpha = .05\) we would reject the null, but at \(\alpha = .01\) we wouldn’t. Same data, different conclusion!
6.7 Failing to reject
Now let’s ask whether a majority of never-married Americans favor the death penalty for murder. In the case of yes/no questions like this, the most obvious null hypothesis is .5 or 50%. This is because the majority wins.
nm <- gss |>
filter(marital == 5, !is.na(cappun)) |>
mutate(favor = as.numeric(cappun == 1))
p_nm <- mean(nm$favor)
n_nm <- nrow(nm)
z_nm <- (p_nm - 0.5) / sqrt(0.5 * 0.5 / n_nm)
c(n = n_nm, p_hat = p_nm, p_value = 2 * (1 - pnorm(abs(z_nm)))) n p_hat p_value
650.0000000 0.5215385 0.2720952 We fail to reject the null. We’d see a sample majority this big pretty often even if the population were split 50/50. But notice that this doesn’t mean support is 50%. It just means our data can’t tell the difference between 50% and what we observed.2
2 The standard error here uses the null value (.5) rather than \(\hat{p}\). That’s because we’re asking what would happen in a world where the null is true. A confidence interval uses \(\hat{p}\) instead.
6.8 p-values and confidence intervals
There is a close relationship between p-values and confidence intervals. For example, if a 95% confidence interval includes the null value, the p-value of the hypothesis test will be above .05.
t.test(tv, mu = 3)$conf.int[1] 3.163060 3.441958
attr(,"conf.level")
[1] 0.95Since this interval does not include 3, we could decide to reject the null hypothesis on that basis.
In fact, because p-values and CIs can both be used for testing, it’s usually better to use CIs because they convey the uncertainty of the estimate as well.
6.9 p-value pitfalls
People sometimes say and do stupid things with p-values. Here are some tips:
- Don’t use asterisks as informal indicators of “how big” an effect is (do you have an alpha level or not?)
- Don’t mistake “statistical significance” for importance
- Don’t run a bunch of tests and only report the significant ones (at \(\alpha = .05\), about one in twenty will be significant by chance even if nothing is going on!)
- When assessing the plausibility of a hypothesis, you need to know the prior probability of the hypothesis as well as the p-value
6.10 Recap
- a hypothesis test is a model comparison: \(H_0\) fixes a value (model C), \(H_1\) estimates one (model A)
- the p-value is the probability of data at least as extreme as ours given that the null is true
- simulating the null world and computing the p-value analytically give the same answer, up to the resolution of the simulation
- alpha is fundamentally arbitrary
- statistically significant means only “p below threshold”
- a confidence interval contains exactly the nulls a test would not reject