Types of errors

Hypothesis testing is not error-proof. We start from the assumption that \(H_0\) is true unless there is strong evidence to reject \(H_0\). Rejecting \(H_0\) when \(H_0\) is true is a type I error (false positive), while retaining \(H_0\) when \(H_1\) is true is called a type II error (false negative).

Power function

The power function of a test with rejection region \(R\) is defined by

\[\beta(\theta) = \Pr_\theta (X \in R)\]

The size of a test is defined to be

\[\alpha = \text{sup}_{\theta \in \Theta_0} \beta(\theta)\]

• \(\text{sup}\) - supremum, or the largest value that \(\beta(\theta)\) could take on in the given \(\theta \in \Theta_0\)

A test is said to have level \(\alpha\) if its size is less than or equal to \(\alpha\).

Sided tests

A test of the form

\[H_0: \theta = \theta_0 \quad \text{versus} \quad H_1: \theta \neq \theta_0\]

is called a two-sided test, or a simple hypothesis. A test of the form

\[H_0: \theta \leq \theta_0 \quad \text{versus} \quad H_1: \theta > \theta_0\]

or

\[H_0: \theta \geq \theta_0 \quad \text{versus} \quad H_1: \theta < \theta_0\]

is called a one-sided test, or a composite hypothesis.

Example hypothesis test

Let \(X_1, \ldots, X_n \sim N(\mu, \sigma^2)\) where \(\sigma\) is known. We want to test \(H_0: \mu \leq 0\) versus \(H_1: \mu > 0\). Hence, \(\Theta_0 = (-\infty, 0]\) and \(\Theta_1 = (0, \infty]\). Consider the test

\[\text{reject } H_0 \text{ if } T>c\]

where \(T = \bar{X}\). The rejection region is

\[R = \left\{(x_1, \ldots, x_n): T(x_1, \ldots, x_n) > c \right\}\]

Let \(Z\) denote the standard Normal random variable. The power function is

\[ \begin{align} \beta(\mu) &= \Pr_\mu (\bar{X} > c) \\ &= \Pr_\mu \left(\frac{\sqrt{n} (\bar{X} - \mu)}{\sigma} > \frac{\sqrt{n} (c - \mu)}{\sigma} \right) \\ &= \Pr_\mu \left(Z > \frac{\sqrt{n} (c - \mu)}{\sigma} \right) \\ &= 1 - \Phi \left( \frac{\sqrt{n} (c - \mu)}{\sigma} \right) \end{align} \]

This function is increasing in \(\mu\):

data_frame(x = seq(-3, 3, by = 0.01),
           y = pnorm(x),
           h = x > -0) %>%
  ggplot(aes(x, y, color = h)) +
  geom_line(size = 1) +
  geom_vline(xintercept = 0) +
  scale_color_brewer(type = "qual", labels = expression(H[0], H[1])) +
  annotate(geom = "text", x = -0.25, y = .8, label = expression(beta(mu))) +
  labs(title = "Power function for test statistic",
       x = expression(mu),
       y = expression(alpha),
       color = NULL)

## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'

Hence

\[\alpha = \text{sup}_{\mu \leq 0} \beta(\mu) = \beta(0) = 1 - \Phi \left( \frac{\sqrt{n} (c)}{\sigma} \right)\]

For a size \(\alpha\) test, we set this equal to \(\alpha\) and solve for \(c\) to get

\[c = \frac{\sigma \Phi^{-1} (1 - \alpha)}{\sqrt{n}}\]

We reject \(H_0\) when

\[\bar{X} > \frac{\sigma \Phi^{-1} (1 - \alpha)}{\sqrt{n}}\]

Equivalently, we reject when

\[\frac{\sqrt{n}(\bar{X} - 0)}{\sigma} > z_\alpha\]

where \(z_\alpha = \Phi^{-1} (1 - \alpha)\).

Ideally we find the test with the highest power under \(H_1\) among all size \(\alpha\) tests. In practice, we use many of the same commonly used tests.

Power of the Wald test

Suppose the true value of \(\theta\) is \(\theta_* \neq \theta_0\). The power \(\beta(\theta_*)\) – the probability of correctly rejecting the null hypothesis – is given (approximately) by

\[1 - \Phi \left( \frac{\hat{\theta} - \theta_0}{\widehat{\se}} + z_{\alpha/2} \right) + \Phi \left( \frac{\hat{\theta} - \theta_0}{\widehat{\se}} - z_{\alpha/2} \right)\]

Remember this is a two-tailed test. Essentially we are collecting the probability mass in the center of the standard normal distribution and subtracting that from 1, to get the area in the tails of the distribution. Hence, two-tailed test.

Recall that \(\widehat{\se}\) tends to 0 as the sample size increases. So we can note that:

• The power is large if \(\theta_*\) is far from \(\theta_0\)
• The power is large if the sample size is large

Example: comparing two means

Let \(X_1, \ldots, X_m\) and \(Y_1, \ldots, Y_n\) be two independent samples from populations with means \(\mu_1, \mu_2\) respectively. Let’s test the null hypothesis that \(\mu_1 = \mu_2\). Write this as

\[H_0: \delta = 0 \quad \text{versus} \quad H_1: \delta \neq 0\]

where \(\delta = \mu_1 - \mu_2\). The estimate of \(\delta\) is \(\hat{\delta} = \bar{X} - \bar{Y}\) with estimated standard error

\[\widehat{\se} = \sqrt{\frac{s_1^2}{m} + \frac{s_2^2}{n}}\]

where \(s_1^2\) and \(s_2^2\) are the sample variances. The size \(\alpha\) Wald test rejects \(H_0\) when \(|W| > z_{\alpha / 2}\) where

\[W = \frac{\hat{\delta} - 0}{\widehat{\se}} = \frac{\bar{X} - \bar{Y}}{\sqrt{\frac{s_1^2}{m} + \frac{s_2^2}{n}}}\]

\(t\)-test

expand.grid(x = seq(-4, 4, by = 0.01),
            df = c(1, 3, 5, 10, 30, Inf)) %>%
  as_tibble() %>%
  mutate(y = dt(x, df),
         df = factor(df, labels = c(1, 3, 5, 10, 30, "Normal"))) %>%
  ggplot(aes(x, y, color = df)) +
  geom_line() +
  scale_color_brewer(type = "qual") +
  labs(title = "t-Distribution",
       x = expression(X),
       y = "PDF",
       color = "Degrees of freedom") +
  theme(legend.position = "bottom")

Relationship to confidence intervals

There is a relationship between the Wald test and the \(1 - \alpha\) asymptotic confidence interval \(\hat{\theta} \pm \widehat{\se} z_{\alpha/2}\). The size \(\alpha\) Wald test rejects \(H_0: \theta = \theta_0 \quad \text{versus} \quad \theta \neq \theta_0\) if and only if \(\theta_0 \notin C\) where

\[C = (\hat{\theta} - \widehat{\se}z_{\alpha / 2}, \hat{\theta} + \widehat{\se}z_{\alpha / 2})\]

Thus, testing the hypothesis is equivalent to checking whether the null value is in the confidence interval.

Statistical vs. scientific significance

Rejecting \(H_0\) indicates the result is statistically significant. That is, we have strong evidence to reject \(H_0\). The result or effect size can still be small if our test is powerful. In that situation, we have statistical significance but not necessarily scientific/substantive/practical significance. You should always be concerned with both of these types of significance. Statistical significance alone is not necessarily a useful or informative finding.

Interpreting \(p\)-values

• These values are informal standards. There is no rhyme or reason they have to be so
• A large \(p\)-value is not strong evidence in favor of \(H_0\)
  • \(H_0\) could be true
  • \(H_0\) is false but the test has low power
• \(p\)-value is not \(\Pr (H_0 | \text{Data})\). The \(p\)-value is not the probability that the null hypothesis is true

Calculating \(p\)-values

Suppose that the size \(\alpha\) test is of the form

\[\text{reject } H_0 \text{ if and only if } T(X_n) \geq c_\alpha\]

Then,

\[\text{p-value} = \text{sup}_{\theta \in \Theta_0} \Pr_\theta (T(X^n) \geq T(x^n))\]

where \(x^n\) is the observed value of \(X^n\). If \(\Theta_0 = \{ \theta_0 \}\) then

\[\text{p-value} = \Pr_{\theta_0} (T(X^n) \geq T(x^n))\]

Informally, the \(p\)-value is the probability under \(H_0\) of observing a value of the test statistic the same as or more extreme then what was actually observed.

data_frame(x = seq(-3, 3, by = 0.01),
           y = dnorm(x),
           tails = x < qnorm(p = 0.025) | x > qnorm(p = 0.975)) %>%
  ggplot(aes(x, y)) +
  geom_linerange(aes(ymax = y, ymin = 0, color = tails)) +
  scale_color_brewer(type = "qual", labels = expression(NULL, alpha/2)) +
  geom_vline(xintercept = qnorm(p = c(0.025, 0.975)), linetype = 2) +
  annotate(geom = "text", x = qnorm(p = c(0.025, 0.975)) - .25, y = .35, 
           label = expression(paste("|", -w, "|"), paste("|", w, "|"))) +
  labs(title = "P-value for the Wald test statistic",
       x = NULL,
       y = NULL,
       color = NULL) +
  theme_void()

## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'

Example: cholesterol data

Consider a set of 371 individuals in a health study examining cholesterol levels (in mg/dl). 320 individuals have narrowing of the arteries, while 51 patients have no evidence of heart disease. Is the mean cholesterol different in the two groups?

Let the estimated mean cholesterol levels for the first group be \(\bar{X} = 216.2\) and for the second group \(\bar{Y} = 195.3\). Let the estimated standard error for each group be \(\widehat{\se}(\hat{\mu}_1) = 5.0\) and \(\widehat{\se}(\hat{\mu}_2) = 2.4\). The Wald test statistic is

\[W = \frac{\hat{\delta} - 0}{\widehat{\se}} = \frac{\bar{X} - \bar{Y}}{\sqrt{\widehat{\se}_1^2 + \widehat{\se}_2^2}} = \frac{216.2 - 195.3}{\sqrt{5^2 + 2.4^2}} = 3.78\]

To compute the \(p\)-value, let \(Z \sim N(0,1)\) denote a standard Normal random variable. Then

\[\text{p-value} = \Pr (|Z| > 3.78) = 2 \Pr(Z < -3.78) = 0.0002\]

which is very strong evidence against the null hypothesis.

Example: attitudes towards abortion

• \(H_A\) - In a comparison of individuals, liberals are more likely to favor allowing a woman to obtain an abortion for any reason than conservatives
• \(H_0\) - There is no difference in support between liberals and conservatives for allowing a woman to obtain an abortion for any reason. Any difference is the result of random sampling error.
• Unlike in previous examples, we want to account for moderates

• Say the null hypothesis is correct - there are no differences between ideological groups and attitudes towards abortion. What would the table look like?

  Right to Abortion	Liberal	Moderate	Conservative	Total
  Yes	40.8%	40.8%	40.8%	40.8%
  	(206.45)	(289.68)	(271.32)	(768)
  No	59.2%	59.2%	59.2%	59.2%
  	(299.55)	(420.32)	(393.68)	(1113)
  Total	26.9%	37.7%	35.4%	100%
  	(506)	(710)	(665)	(1881)

• In truth, what does our table actually look like?

  Right to Abortion	Liberal	Moderate	Conservative	Total
  Yes	62.6%	36.6%	28.7%	40.8%
  	(317)	(260)	(191)	(768)
  No	37.4%	63.4%	71.28%	59.2%
  	(189)	(450)	(474)	(1113)
  Total	26.9%	37.7%	35.4%	100%
  	(506)	(710)	(665)	(1881)

• How can we test if these differences are statistically significant? That is, how do we test to see if we can reject the null hypothesis?

• Calculating test statistic
  • \(\chi^2=\sum{\frac{(X_j - E_j)^2}{E_j}}=145.27\)
  • \(\text{Degrees of freedom} = (\text{number of rows}-1)(\text{number of columns-1})=2\)
• Calculating \(p\)-value
  • \(\text{p-value} = \Pr (\chi_2^2 > 145.27) = 0\)
  • The probability that the null hypothesis is true and the observed frequencies are the result of random sampling error is less than 1 in a quintillion. Extremely extremely unlikely the null hypothesis is true.

Example: Mendel’s peas

Mendel bred peas with round yellow seeds and wrinkled green seeds. There are four types of progeny: round yellow, wrinkled yellow, round green, and wrinkled green. The number of each type is multinomial with probability \(p = (p_1, p_2, p_3, p_4)\). His history of inheritance predicts that \(p\) is equal to

\[p_0 \equiv \left(\frac{9}{16}, \frac{3}{16}, \frac{3}{16}, \frac{1}{16} \right)\]

In \(n = 556\) trials he observed \(X = (315,101,108,32)\). The likelihood ratio test for \(H_0: p = p_0 \quad \text{versus} \quad H_1: p \neq p_0\) is

\[ \begin{align} \lambda &= 2 \log \left( \frac{\Lagr (\hat{p})}{\Lagr(p_0)} \right) \\ &= 2 \sum_{i=1}^4 X_j \log \left( \frac{\hat{p}_j}{p_{0j}} \right) \\ &= 2 \left[ 315 \log \left( \frac{\dfrac{315}{556}}{\dfrac{9}{16}} \right) + 101 \log \left( \frac{\dfrac{101}{556}}{\dfrac{3}{16}} \right) + 108 \log \left( \frac{\dfrac{108}{556}}{\dfrac{3}{16}} \right) + 32 \log \left( \frac{\dfrac{32}{556}}{\dfrac{1}{16}} \right)\right] \\ &= 0.48 \end{align} \]

Under \(H_1\) there are four parameters. However, three parameters must sum to one so the dimension of the parameter space is three. Under \(H_0\) there are no free parameters so the dimension of the restricted parameter space is zero. The difference of these dimension is three. Therefore the limiting distribution of \(\lambda\) under \(H_0\) is \(\chi_3^2\) and the p-value is

\[\text{p-value} = \Pr (\chi_3^2 > 0.48) = 0.923\]

A similar result is achieved using the \(\chi^2\) test. When both tests are applicable, they will lead to similar results as long as the sample size is large.

Correcting \(p\)-values

Consider \(m\) hypothesis tests:

\[H_{0i} \quad \text{versus} \quad H_{1i}, i = 1, \ldots, m\]

and let \(P_1, \ldots, P_m\) denote the \(m\) \(p\)-values for these tests. The Bonferroni method approaches the problem as follows: given \(p\)-values \(P_1, \ldots, P_m\), reject null hypothesis \(H_{0i}\) if

\[P_i \leq \frac{\alpha}{m}\]

For example, if a study tested \(m=20\) hypotheses with a desired \(\alpha = 0.05\) (the standard threshold), then the Bonferroni correction would test each individual hypothesis at \(\alpha = \frac{0.05}{20} = 0.0025\).

For example, draw 100 observations from a normal distribution \(N(0,1)\) and test the null hypothesis \(H_0: \mu = 0\) using a t-test.¹ If we simulate this process multiple times, we should reject \(H_0\) approximately 5% of the time.

sim_norm_null <- 1000 %>%
  rerun(rnorm(n_obs)) %>%
  map(~ t.test(x = .x, mu = 0)) %>%
  map_dbl(~ .x$p.value) %>%
  as_tibble %>%
  mutate(sig = value < .05)

mean(sim_norm_null$value)

## [1] 0.498

ggplot(sim_norm_null, aes(value, fill = sig)) +
  geom_histogram(binwidth = .025, boundary = 0) +
  labs(title = "Distribution of p-values for single test",
       x = expression(P),
       y = NULL) +
  theme(legend.position = "none")

Now let’s simulate 5 random variables and test the null hypothesis that all means are simultaneously 0, then the probability of at least one significant result is \(1 - (1 - 0.05)^5 = 0.226\):

sim_norm_mult <- 1000 %>%
  rerun(5 %>%
          rerun(rnorm(n_obs)) %>%
          map(~ t.test(x = .x, mu = 0)) %>%
          map_dbl(~ .x$p.value) %>%
          as_tibble %>%
          mutate(sig = value < .05)) %>%
  bind_rows(.id = "sim") %>%
  group_by(sim) %>%
  rename(raw = value) %>%
  mutate(correct = raw < (.05 / n()))

sim_norm_mult %>%
  summarize(sig = any(raw < .05)) %>%
  ungroup %>%
  summarize(mean(sig))

## # A tibble: 1 x 1
##   `mean(sig)`
##         <dbl>
## 1       0.239

sim_norm_mult %>%
  filter(raw == min(raw)) %>%
  ggplot(aes(raw, fill = sig)) +
  geom_histogram(binwidth = .01, boundary = 0) +
  labs(title = "Distribution of p-values for multiple tests",
       x = expression(P),
       y = NULL) +
  theme(legend.position = "none")

But if we use the Bonferroni correction:

sim_norm_mult %>%
  summarize(sig = any(correct)) %>%
  ungroup %>%
  summarize(mean(sig))

## # A tibble: 1 x 1
##   `mean(sig)`
##         <dbl>
## 1       0.052

sim_norm_mult %>%
  filter(raw == min(raw)) %>%
  ggplot(aes(raw, fill = correct)) +
  geom_histogram(binwidth = .01, boundary = 0) +
  labs(title = "Distribution of p-values for multiple tests",
       subtitle = "With Bonferroni correction",
       x = expression(P),
       y = NULL) +
  theme(legend.position = "none")

The Bonferroni correction is a conservative adjustment, and errs on the side of caution. You minimize the risk of a false positive, but therefore increase the risk of a false negative. Alternative correction methods such as the Benjamini–Hochberg procedure are less conservative.