Inference for 1 Categorical Variable - Single Proportion

# Loading the tidyverse package
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.1     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.3     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Reading in the NFL drives data from github

# Read in the nfl drive data
drives <- read.csv("https://raw.githubusercontent.com/Shammalamala/STA4504/refs/heads/main/data/ch1/nfl%20drives.csv")

# Looking at the different ways drives can end
unique(drives$drive_end)

## [1] "Field Goal" "Turnover"   "Touchdown"  "Punt"

We will be focusing on the variable drive_end

There are 4 ways in the data that a touchdown can end: Touchdown, Field Goal, Punt, Turnover

Hypothesis Test

Let’s perform a hypothesis test to answer the question “Do 20% of all drives end in touchdowns?”

\[H_0: \pi = 0.20 \\ H_1: \pi \ne 0.20\]

For a single parameter, test statistics follow the same general pattern:

\[z = \frac{\textrm{statistic} - \textrm{null value}}{\textrm{standard error}}\]

If we are interested in learning about a single proportion, our test statistic is:

\[z = \frac{\hat{\pi} - \pi_0}{\sqrt{\frac{\pi(1-\pi)}{n}}}\]

We have \(\pi_0\), but need to calculate our sample proportion:

# Calculating the total number of tds
td_total <- sum(drives$drive_end == "Touchdown")

# Calculating the sample proportion using mean(drives$drive_end = "Touchdown")
td_prop <- mean(drives$drive_end == "Touchdown")

td_prop

## [1] 0.24

Let’s save the sample size, \(n\), and the null hypothesis value, \(\pi_0\):

n <- nrow(drives); pi0 <- 0.20

We have a sample proportion of \(\hat{\pi} = 0.24\) (24%). Our null hypothesis is 0.20 (\(\pi_0 = 0.20\)).

So we have the top of our test statistic fraction, but what about the denominator:

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

We don’t have \(\pi\) (hence the hypothesis test). So what do we replace it with?

There are two choices:

1. Replace the unknown parameter with the sample statistic: \(\pi \rightarrow \hat{\pi}\)

• This is our Wald standard error

2. Replace the unknown parameter with the null hypothesis value: \(\pi \rightarrow \pi_0\)$

• This is our Score, sometimes called Wilson, standard error

For larger sample sizes, the two won’t be that different. But for smaller sample sizes, the standard errors can be very different.

When performing a hypothesis test, a Score standard error is typically used and the test is then called a Score Test (unsurprising!)

Wald Test and Score Tests

Let’s start with the Wald test:

\[z = \frac{\hat{\pi} - \pi_0}{\sqrt{\frac{\hat{\pi}(1-\hat{\pi})}{n}}}\]

wald_se <- sqrt(td_prop * (1 - td_prop)/n)
wald_se

## [1] 0.01207974

Next, we’ll find the test statistic:

wald_z <- (td_prop - pi0) / wald_se
wald_z

## [1] 3.311331

Then we can find the p-value: \(P(|Z| > 3.311)\)

wald_pval <- 2 * pnorm(abs(wald_z), lower.tail = F)
wald_pval

## [1] 0.0009285334

The Wald test rejects the null hypothesis!

Score Test

Up next: Score Test

\[z = \frac{\hat{\pi} - \pi_0}{\sqrt{\frac{\pi_0(1-\pi_0)}{n}}}\]

score_se <- sqrt(pi0 * (1 - pi0)/n)
score_se

## [1] 0.01131371

Next, we’ll find the test statistic:

score_z <- (td_prop - pi0) / score_se
score_z

## [1] 3.535534

Then we can find the p-value: \(P(|Z| > 3.536)\)

score_pval <- 2 * pnorm(abs(score_z), lower.tail = F)
score_pval

## [1] 0.000406952

The test statistic is larger and the p-value smaller for the score test compared to the Wald test!

Likelihood Ratio Test

The likelihood ratio test (LRT) takes a different approach than the typical \(\frac{\hat{\pi} - pi_0}{SE}\) test statistic. Instead, for a binomial random variable, it is a ratio of two different binomial distributions:

\[\frac{\ell_1}{\ell_0} = \frac{{n \choose y} \hat{\pi}^y(1-\hat{\pi})^{n-y}}{{n \choose y} (\pi_0)^y(1-\pi_0)^{n-y}}\]

We can find the numerator and denominator using dbinom(...) with prob = the corresponding probabilities

# unrestricted likelihood
ell1 <- dbinom(td_total, size = n, prob = td_prop)

# restricted likelihood
ell0 <- dbinom(td_total, size = n, prob = pi0)

# LRT test stat 
lrt_test_stat <- ell1 / ell0
lrt_test_stat

## [1] 390.6599

In order to find a p-value, we need to know that distribution the test statistic follows.

While \(\ell_1 / \ell_0\) itself doesn’t have a defined distribution, we can use Wilk’s Theorem to find a test statistic and distribution:

\[2\log\left(\frac{\ell_1}{\ell_0}\right) \sim \chi^2_v\]

For a single binomial random, there is one ‘unrestrained’ parameter in \(\ell_1\), so \(v = 1\)

The LRT p-value is:

\[P(\chi^2_1 > 11.936)\]

log_lrt_test_stat <- 2 * log(lrt_test_stat)
lrt_pval <- pchisq(q = log_lrt_test_stat, df = 1, lower.tail = F)
lrt_pval

## [1] 0.0005506918

The p-value for the LRT test is similar to the score test statistic.

Comparing the three tests:

If we compare the three exams:

tibble(
  test = c('Wald', 'Score', 'LRT'),
  test_stat = round(c(wald_z, score_z, log_lrt_test_stat), 5),
  p_value = round(c(wald_pval, score_pval, lrt_pval), 5)
)

## # A tibble: 3 × 3
##   test  test_stat p_value
##   <chr>     <dbl>   <dbl>
## 1 Wald       3.31 0.00093
## 2 Score      3.54 0.00041
## 3 LRT       11.9  0.00055

Built-in functions for tests:

Note: Regardless of which test we use, the functions will return the \(\chi^2\) version of the test statistic, \(z^2\).

Wald test:

There isn’t one :(

That’s because using the Wald test for a single proportion is a bad idea. We did it (and you’ll do it on the homework) to show that it exists, but you shouldn’t use it!

Why?

What happens if \(y = n\)?

Score test:

For the score test, you can use prop.test() in base R, or prop_test() in rstatix. prop_test() is very similar to prop.test(), but it gives the results in a data frame that is easier to work with than the object created by prop.test().

Both functions have the following arguments:

• x = the number of successes (we denote as \(y\))
• n = the sample size
• p = the null hypothesis value, \(\pi_0\)
• correct = F to have it do a score test without a binomial correction

prop.test(
  x = td_total,
  n = n, 
  p = pi0,
  alternative = 'two.sided',
  correct = F
)

## 
##  1-sample proportions test without continuity correction
## 
## data:  td_total out of n, null probability pi0
## X-squared = 12.5, df = 1, p-value = 0.000407
## alternative hypothesis: true p is not equal to 0.2
## 95 percent confidence interval:
##  0.2171436 0.2644495
## sample estimates:
##    p 
## 0.24

rstatix::prop_test(
  x = td_total,
  n = n, 
  p = pi0,
  correct = F
) 

## # A tibble: 1 × 5
##       n statistic    df        p p.signif
## * <int>     <dbl> <int>    <dbl> <chr>   
## 1  1250      12.5     1 0.000407 ***

LRT Built-in:

To conduct a likelihood ratio test using a built-in function, we need to install the DescTools package and use the GTest(...)

The GTest(...) function doesn’t do a one-proportion test, but it will do a goodness-of-fit test. What good is that to us?

If the number of categories in a GoF test is 2, then it is equivalent to a one prop test!

That does mean instead of just giving it the number of successes and the hypothesized probability, we need to give it both \(y\) and \(n-y\) along with \(\pi_0\) and \(1-\pi_0\) in two vectors for x and p, as seen below:

DescTools::GTest(
  x = c(td_total, n - td_total),
  p = c(pi0, 1-pi0)
)

## 
##  Log likelihood ratio (G-test) goodness of fit test
## 
## data:  c(td_total, n - td_total)
## G = 11.936, X-squared df = 1, p-value = 0.0005507