Please beware of the problem of testing too many hypotheses; the more you torture the data, the more likely they are to confess, but confessions obtained under duress may not be admissible in the court of scientific opinion.

—Stephen Stigler

Week 4: Hypothesis testing

Once we have estimated statistical parameters in a model (week 1-3) and are confident that we have appropriately described the data, we can begin to make a statistical inference. Inferential statistics is concerned with inferring properties from a sample of the population to the entire population. In particular, conducting a hypothesis test on a parameter assesses the evidence (from the known observations) that a parameter takes a certain value (or is in a set of values).

Hypothesis testing is a cornerstone of frequentist statistics¹ and modern inductive science². We focus on likelihood ratio tests (LRT), which compare the evidence supporting the null hypothesis versus the alternative hypothesis. There are two cases we consider.

Likelihood Ratio Tests

Simple Hypotheses

A simple-versus-simple hypothesis test occurs the null, and the alternative hypothesis consists of single values. More specifically \[ H_{0}:\theta = \theta_{0}\mbox{ and } H_{1}:\theta = \theta_{1} \] The likelihood ratio is

\[\Lambda = \frac{L(\theta_{0}~\vert~\boldsymbol{y})}{L(\theta_{1}~\vert~\boldsymbol{y})}.\] For this case, the Neyman–Pearson Lemma states that this likelihood-ratio test is the most powerful for all $\alpha$-level tests. Note: $\Lambda$ is a statistic and depends on the data; as such, before we can treat it as a random variable before observing the data. For a fixed constant $c$, the decision rule for this hypothesis test is:

Do not reject $H_{0}$ if $\Lambda > c$
Reject³ $H_{0}$ if $\Lambda < c$

The constant $c$ is chosen with respect to the hypothesis test’s properties (see ``Properties of hypothesis tests’’). As such, the method for determining $c$ is

Assume the null hypothesis, $H_{0}$, is true.
Simplify the rejection region so that it involves a quantity for which the distribution is known.
Choose $c$ to decide the Type I and/or Type II error based on part 2.

The hypothesis test assumes $H_{0}$ is true so that evidence against the null occurs when we make unlikely observations (with respect to the probability distribution).

Question 1 explores this further.

More on this later in the semester. But for a jump start, see here.↩︎
Despite being the dominating paradigm for assessing scientific evidence, like any method, it is not without flaws and can be misused.↩︎
The rejection region is defined as $\Lambda < c$ but can usually be simplified.↩︎

Composite Hypotheses

We cab handle other hypotheses using the generalised likelihood test (GLRT). This hypothesis test partitions the parameter space, $\Theta$, into two sets, $\Theta_{0}$ and $\Theta_{1}$:

The null hypothesis: $H_{0} = \{\theta:\theta \in \Theta_{0}\}$, and
The alternative hypothesis⁴: $H_{1} = \{\theta:\theta \in \Theta_{1}\}$

The ratio for this test is: \[ \Lambda = \frac{\sup_{\theta \in \Theta_{0}} L(\theta~\vert~\boldsymbol{y})}{\sup_{\theta \in \Theta} L(\theta ~\vert~\boldsymbol{y})}. \] We compare the likelihoods for the most likely parameter under the null hypothesis to the most likely parameter in the entire parameter space.

Note that simple likelihood ratio tests are a specific type of GLRT.

Note that $\Theta_{1}$ is the complement of $\Theta_{0}$, $\Theta_{1}^{\complement} = \Theta_{0}$. Hence the two sets cover the entire parameter space.↩︎

Type of Errors

Considering the properties of hypothesis tests (just like the properties of estimators) is also important. Important properties include:

Type I error (Significance): \[ \alpha = P(\text{reject } H_{0}~\vert~H_{0} \text{ is true}) \]
Type II error: \[ \beta = P(\text{do not reject } H_{0}~\vert~H_{1} \text{ is true}) \]
Test power: \[ \text{power} = 1 - \beta = P(\text{reject } H_{0}~\vert~H_{1} \text{ is true}). \] It is important to consider these error rates and power when designing test and experiments. Because Type I error is considered the more serious error, we usually fix that rate and then adjust other parameters (typically difference between the true value and hypothesised value and sample size) to obtain the desired power.

Wilk’s theorem

For large samples (i.e. $n\rightarrow \infty$), a log-likelihood ratio has an approximate Chi-squared distribution \[ 2 \left( \ell(\hat{\boldsymbol{\theta}}) - \ell(\boldsymbol{\theta}) \right) \sim \chi^{2}_{p} \] where $p$ is the number of parameters estimated. Often it is neccessary to use a Taylor series approximation on $2 \left( \ell(\hat{\boldsymbol{\theta}}) - \ell(\boldsymbol{\theta}) \right)$ because it is hard to invert. This leads to the approximate distribution \[ \left( \hat{\boldsymbol{\theta}} - \boldsymbol{\theta} \right)^{\top}\mathcal{J} \left( \hat{\boldsymbol{\theta}} - \boldsymbol{\theta} \right) \sim \chi^{2}_{p} \] where $\mathcal{J}$ is the observed information matrix.

Theory questions

Question 1

Let $y_{1},y_{2},\ldots,y_{n}$ be a random sample of observations from a population described by a normal distribution with mean equal to zero and variance $\sigma^2$.

What is the probability distribution of $T = \frac{1}{\sigma^2}\sum_{i=1}^n y_i^2$? State what it is, no need to show proof.
Determine the likelihood ratio, $\Lambda$, for the most powerful test of $H_0: \sigma^2=\sigma^2_0$ against $H_1:\sigma^2=\sigma^2_1$.
Simplify the rejection region⁵, $\Lambda < c$, in terms of the statistic $T$. Note that $c$ is a constant⁶ with $0<c<1$.
Determine when the null hypothesis will be rejected using significance level $\alpha$.

Hint: Consider cases $\sigma^2_1 > \sigma^2_0$ and $\sigma^2_1 < \sigma^2_0$ separately.↩︎
Hint: The constant $c$ may be combined with other constants to simplify the derivation.↩︎

Question 2

Let $y_1,y_2,\ldots,y_n$ be a random sample from a population described by the normal $N(\theta,a\theta^2)$ distribution, that is the variance of the distribution is $a\theta^2$ for $a > 0$. Determine the form of the generalised likelihood ratio test of $H_0:a=1$ against $H_1:a\ne 1$. Start by finding the MLEs of the distribution⁷.

Solution:

The MLEs are $\hat{\theta}=\bar{y}$ and $\hat{\sigma}^2 = \frac{1}{n}\sum_{i=1}^n (y_i-\bar{y})^2$. Thus if $\sigma^2=a\theta^2$, then $\hat{a}=\hat{\sigma}^2/\hat{\theta}^2=\frac{\sum_{i=1}^n (y_i-\bar{y})^2}{\sum_{i=1}^n y_i}$. The likelihood for the null hypothesis is $L_{0}$, which is found by writing the likelihood of the model with $a = 1$ and $\theta = \hat{\theta}.$: \[ L_{0} = \prod_{i=1}^{n}(2\pi \hat{\theta}^2)^{-1/2}\exp\left\{-\frac{(y_{i}-\hat{\theta})^{2}}{2\hat{\theta}^2}\right\} \] \[ =(2\pi \bar{y}^2)^{-n/2}\exp\left\{-\sum_{i=1}^{n}\frac{(y_{i}-\bar{y})^{2}}{2\bar{y}^{2}}\right\} \] \[ =(2\pi \bar{y}^2)^{-n/2}\exp\left\{-\frac{1}{2}\left(\frac{\sum_{i=1}^{n}y_{i}^{2}}{\bar{y}^{2}} - n\right)\right\} \] The likelihood for the alternative hypothesis is $L_{1}$: \[ L_{1} = \prod_{i=1}^{n}(2\pi \hat{a}\hat{\theta}^{2})^{-1/2}\exp\left\{-\frac{(y_{i}-\hat{\theta})^{2}}{2\hat{a}\hat{\theta}}\right\} \] \[ =(2\pi\hat{a} \bar{y}^2)^{-n/2}\exp\left\{- \frac{n}{2}\right\}\\ \] The GLRT can be simplified as \[ \Lambda =\frac{L_0}{L_1} = \hat{a}^{n/2}\exp \left(-\frac{1}{2} \frac{\sum_{i=1}^{n}y_{i}^{2}}{\bar{y}^{2}} \right). \] Note that \[ \hat{a} = \frac{\sum_{i=1}^{n}(y_{i} - \bar{y})^{2}}{n\bar{y}^{2}} = \frac{\sum_{i=1}^{n}y_{i}^{2}}{n\bar{y}^{2}} - 1 \] To use the Wilk's theorem for the asymptotic test, we need \[ w = -2 \log(\Lambda) = n\log\hat{a} - \frac{\sum_{i=1}^{n}y_{i}^{2}}{\bar{y}^{2}} \] and we would reject $H_{0}$ if $w > c$ where $c$ is chosen such that $P(X > c) = \alpha$ for an $\alpha$-level test and $X \sim \chi^{2}_{1}$. For an exact (non-asymptotic) test, the critical region is \[ \hat{a}^{n/2}\exp \left(-\frac{1}{2} \frac{\sum_{i=1}^{n}y_{i}^{2}}{\bar{y}^{2}} \right) < c_{1}\] or equivalently \[ n\log\left(\frac{\sum_{i=1}^{n}y_{i}^{2}}{n\bar{y}^{2}} - 1\right) - \frac{\sum_{i=1}^{n}y_{i}^{2}}{\bar{y}^{2}} < c_{2} \] Which we doesn't simplify any further. We can use this exact test if the distribution of the LHS can be found, which is difficult; but, it has been formulated in terms of a sufficient statistic - namely, $t = \frac{\sum_{i=1}^{n}y_{i}^{2}}{\bar{y}^{2}}$ - for which approximate methods could be used.

Hint: the maximum likelihood estimates from the normal distribution $N(\mu,\sigma^{2})$ are $\hat{\mu} = \bar{y}$ and $\hat{\sigma}^{2} = n^{-1}\sum_{i=1}^{n} (y_{i} - \bar{y})^{2}$. The MLEs have an invariance property, so that the MLEs of $\theta$ and $a$ can be determined from those of $\theta$ and $a\theta$ by transformation.↩︎

Practical questions

Note on random variables in R

For this section, you will need to use random variables in R. In base R; these have a particular naming convention:

d<rv>(): probability density (or mass) function. The univariate normal distribution has dnorm(x, mean = 0, sd = 1) for example.
p<rv>(): distribution function, $P(X < q)$. The chi-squared has pchisq(q, df, lower.tail = TRUE) for example.
q<rv>(): quantile function, $q$ such that $P(X < q) = p$. The binomial distribution has qbinom(p, size, prob, lower.tail = TRUE) for example.
r<rv>(): random variable generator. The exponential distribution has rexp(n, rate = 1) for example.

Question 3

Suppose you observe the data, $\boldsymbol{y}$ given below, and you wish to perform the hypothesis test in Question 1 (theory) with \[ H_{0}: \sigma^{2} = 100$ \mbox{ versus } $H_{1}: \sigma^{2} = 20. \]

Do you reject the null hypothesis at $\alpha = 0.05$? What is the $p$-value⁸ of the test?
Plot the distributions relating to the null and alternative hypothesis, along with where the observations occurred. Discuss how this graph relates to the hypothesis test.

The data are: $\mathbf{y}=(9.91, 4.37, -1.06, 7.63, 0.78, 2.73, -1.29, 2.78, -4.98, 2.45)$

Solution


    
    y <- c(9.91, 4.37, -1.06, 7.63, 0.78, 
           2.73, -1.29, 2.78, -4.98, 2.45)
    
    n <- length(y)
    
    sigma2_null <- 100
    sigma2_alt <- 20
    
    alpha_level <- 0.05
    
    
    test_stat <- sum(y^2) / sigma2_null
    
    if( sigma2_null > sigma2_alt){
      # P(X < c) = alpha
      c_critical_val <- qchisq(p = alpha_level, 
                               df = n, lower.tail = TRUE)
      reject_decision <- test_stat < c_critical_val
      
      pval <- pchisq(q = test_stat, 
                     df = n, lower.tail = TRUE)
      
    } else {
      # P(X > c) = alpha
      c_critical_val <- qchisq(p = alpha_level, 
                               df = n, lower.tail = FALSE)
      reject_decision <- test_stat > c_critical_val
      
      pval <- pchisq(q = test_stat, 
                     df = n, lower.tail = FALSE)
      
    }
    
    paste("Decision:", 
          ifelse(reject_decision, "Reject", "Do not reject"), 
          "null hypothesis")
    
    # p-value
    # you can use the p-value to decide on rejection too by
    # comparing it to alpha
    round(pval, 3)

The $p$-value of a test is the probability of observing a test statistic equal to or more extreme than what was actually observed. In this case, the critical region is in the left-tail ($\sigma^{2}_{0} < \sigma^{2}_{1}$), so the p-value is $P(T < t~\vert~H_{0} \text{ is true})$ where $t$ is the observed statistic.↩︎

Question 4

Suppose you wish to conduct an experiment to gauge whether male and female Leonberger dogs have different rates of hip dysplasia⁹. Given you have access to a random sample of Leonbergers where 4 out of 20 males and 2 out of 25 females have hip dysplasia, answer the following questions

Conduct a generalised likelihood ratio test¹⁰. to consider the evidence that the probability of disease occurrence is the same for male and female Leonbergers. Use Wilk’s theorem to approximate the distribution of the ratio.
What is the $p$-value of the test?
Assuming the MLEs are the true parameter values, conduct a simulation study to approximate the power of this (asymptotic) hypothesis test.
Imagine you wish to conduct another sample of Leonberger dogs to test for the occurrence of dysplasia. By adapting the simulation in c. can you determine the sample size $n$ (number of dogs) required for the power of the test to be 0.9? Assume that you will recruit an equal number of male and female dogs.

Wee jock and Loulou with deflated ball.

Solution

a. and b. From previous worksheets, the Bernoulli likelihood is \[ L_{0}(p~\vert~s) = p^{s}(1-p)^{n-s} \] where $n$ is the number of observations and $s = \sum_{i=1}^{n}y_{i}$ is the number of successes. The MLE is $\hat{p} = s/n$. The likelihood in this scenario has two independent populations (males and females) so the likelihood is \[ L_{1}(p_{1},p_{2}~\vert~s) = p_{1}^{s_{1}}(1-p_{1})^{n_{1}-s_{1}}p_{2}^{s_{2}}(1-p_{2})^{n_{2}-s_{2}} \] where subscript 1 denotes the parameters/values for the females dogs, and subscript 2 is for the male dogs. Since the parameters $p_{1}$ and $p_{2}$ don't depend on each other, the MLEs will be $\hat{p}_{1} = s_{1}/n_{1}$ and $\hat{p}_{2} = s_{2}/n_{2}$. The hypothesis test should be $H_{0}: p_{1}=p_{2}$ and $H_{1}: p_{1}\neq p_{2}$. The GLRT is therefore: \[ \Lambda = \frac{L_{0}(\hat{p}~\vert~s) }{L_{1}(\hat{p}_{1},\hat{p}_{2}~\vert~s)} \] and the (asymptotic) test statistic is \[ w = -2\log(\Lambda) = -2 \left[ \log L_{0}(\hat{p}~\vert~s) - \log L_{1}(\hat{p}_{1},\hat{p}_{2}~\vert~s) \right] \] we would reject $H_{0}$ if $w > c$ where $c$ is chosen such that $P(X > c) = \alpha$ for an $\alpha$-level test and $X \sim \chi^{2}_{1}$. The degrees of freedom is $2 - 1 = 1$ because the number of parameters in $H_{1}$ is 2 and in $H_{0}$ is 1.


    (a) and (b) solution
    # likelihood for Bernoulli observations
    bern_likelihood <- function(p,s,n){
      ( p^s ) * ( (1-p)^(s-n) )
    }
    
    # critical value level
    alpha <- 0.05
    
    # sufficient statistics (number of successes)
    s_1 <- 4
    s_2 <- 2
    s <- s_1 + s_2
    
    # number of observations
    n_1 <- 20
    n_2 <- 25
    n <- n_1 + n_2
    
    # MLEs under the two hypotheses
    p_hat_null <- s/n
    p_1_hat_alt <- s_1/n_1
    p_2_hat_alt <- s_2/n_2
    
    # likelihoods under null and alt hypothesis
    L_0 <- bern_likelihood(p = p_hat_null, s = s, n = n)
    L_1 <- bern_likelihood(p = p_1_hat_alt, s = s_1, n = n_1) *
      bern_likelihood(p = p_2_hat_alt, s = s_2, n = n_2)
    
    # - 2 log(LR)
    w_log_ratio <- - 2 * (log(L_0) - log(L_1))
    
    # cutoff based on chi-squared distribution
    test_stat_chi <- qchisq(p = alpha, df = 1, lower.tail = FALSE)
    
    reject_decision <- w_log_ratio > test_stat_chi
    
    paste("Decision:", ifelse(reject_decision, "Reject", "Do not reject"), 
          "null hypothesis (Caution: based on large n approximation)") 
    
    # part (b): (approximate) p-value
    pval <- pchisq(q = w_log_ratio, df = 1, lower.tail = FALSE)
    # you can use the p-value to decide on rejection too by
    # comparing it to alpha
    round(pval,3)


    
    (c) solution
    sim_bern_hypothesis_test <- function(true_p_1, true_p_2, n_1, n_2){
  
    s_1 <- rbinom(n = 1, size = n_1, prob = true_p_1)
    s_2 <- rbinom(n = 1, size = n_2, prob = true_p_2)
    s <- s_1 + s_2
    
    n <- n_1 + n_2
    
    # MLEs under the two hypotheses
    p_hat_null <- s/n
    p_1_hat_alt <- s_1/n_1
    p_2_hat_alt <- s_2/n_2
    
    # likelihoods under null and alt hypothesis
    L_0 <- bern_likelihood(p = p_hat_null, s = s, n = n)
    L_1 <- bern_likelihood(p = p_1_hat_alt, s = s_1, n = n_1) *
      bern_likelihood(p = p_2_hat_alt, s = s_2, n = n_2)
    
    # - 2 log(LR)
    w_log_ratio <- - 2 * (log(L_0) - log(L_1))
    
    # cutoff based on chi-squared distribution
    test_stat_chi <- qchisq(p = alpha, df = 1, lower.tail = FALSE)
    
    reject_decision <- w_log_ratio > test_stat_chi
    
    return(reject_decision)
    
  }
  
  sim_results <- 
      rerun(.n = 1000, 
           sim_bern_hypothesis_test(
              true_p_1 = p_1_hat_alt, 
              true_p_2 = p_2_hat_alt,
              n_1 = n_1,
              n_2 = n_2
              ) 
            ) %>%
      simplify()
    
    # Very low power under these sample sizes
    mean(sim_results)

d. One way to answer is to repeat the simulation above with $n_{1} = n_{2} = b$, with a large enough $b$ to get the desired power.

Hip dysplasia is a serious condition in many giant dog breeds, see here, but the occurrence has been reduced in several types of dogs by careful breeding.↩︎
Hint: Assume a Bernoulli probability distribution for occurrence of the disease.↩︎