Confidence Intervals for a Single Proportion

Loading the packages:

# Loading the tidyverse package
library(tidyverse)

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

Summary Statistics

Let’s start by getting the needed statistics from the sample: \(y\), \(n\), and \(\hat{\pi}\)

# total touchdowns
td_total <- sum(drives$drive_end == "Touchdown")

# sample size
n <- nrow(drives)

# Sample proportion
td_prop <- mean(drives$drive_end == "Touchdown")

c('y' = td_total, 'n' = n, 'pi_hat' = td_prop)

##       y       n  pi_hat 
##  300.00 1250.00    0.24

Confidence Intervals for a proportion

If we reject our null hypothesis, we often followup with a confidence interval to narrow down what \(\pi\) could be.

We’ll look at 3 different intervals:

1. Wald interval - the most common interval

2. Score interval - Equivalent to the score test

3. Agresti-Coull interval - Tends to work better than the previous two for smaller samples when the sample proportion is close to 0 or 1

We’ll be using a 95% confidence level for our examples, but these work for any chosen level!

Wald Interval

\[\hat{\pi} \pm z_{\alpha/2} \sqrt{\frac{\hat{\pi}(1-\hat{\pi})}{n}}\]

# Start by finding the critical value of a 95% CI:
z95 <- qnorm(1 - 0.05/2)

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

# Next the confidence interval!
wald_ci_td <- 
  c(
    lower = td_prop - z95*wald_se,
    upper = td_prop + z95*wald_se
  )

wald_ci_td |> signif(digits = 3)

## lower upper 
## 0.216 0.264

Using a Wald confidence interval, we can be 95% confident that between 21.6% to 26.4% of all NFL drives end in touchdowns.

Score CI

How do we calculate a Score confidence interval?

We find the values of \(\pi_0\) that would cause us to not reject the null hypothesis. The easiest way to find that is find the values of \(\pi_0\) that give a test stat less than or equal to our critical value we found earlier!

score_df <- 
  tibble(
    pi0 = seq(0, 1, by = 0.001),    # Create a range of pi_0 values to test
    se0 = sqrt(pi0*(1-pi0)/n), # Calculating the standard errors for the different pi0
    z = (td_prop - pi0)/se0    # Finding the Score test stat for the various pi0
  )        

score_df

## # A tibble: 1,001 × 3
##      pi0      se0     z
##    <dbl>    <dbl> <dbl>
##  1 0     0        Inf  
##  2 0.001 0.000894 267. 
##  3 0.002 0.00126  188. 
##  4 0.003 0.00155  153. 
##  5 0.004 0.00179  132. 
##  6 0.005 0.00199  118. 
##  7 0.006 0.00218  107. 
##  8 0.007 0.00236   98.8
##  9 0.008 0.00252   92.1
## 10 0.009 0.00267   86.5
## # ℹ 991 more rows

# Next we want to find the values of pi0 that have a test stat below z95
score_df |> 
  filter(abs(z) <= z95) |> 
  summarize(
    lower = min(pi0),
    upper = max(pi0)
  )

## # A tibble: 1 × 2
##   lower upper
##   <dbl> <dbl>
## 1 0.218 0.264

When rounding to 3 decimal places, we get a very similar confidence interval to the Wald interval earlier!

There is also a formula we can use in place of conducting a grid search:

\[\frac{y +\left(z^*\right)^2/2}{n + \left(z^*\right)^2} \pm \frac{z^*}{n + \left(z^*\right)^2} \sqrt{\frac{y(n - y)}{n} + \frac{\left(z^*\right)^2}{4}}\]

where \(y\) is the count, not proportion

# First, find the "statistic"
score_stat <- (td_total + z95^2/2)/(n + z95^2)


# Next find the standard error
score_se <- sqrt(td_total *(n - td_total)/n + z95^2/4)


# Then margin of error
me_formula <- z95/(n + z95^2) * score_se


# Then interval:
tribble(
      ~'result',                ~'value',
            'y',                td_total,
            'n',                       n,
   'score stat',              score_stat,
     'Score SE',                score_se,
     'Score ME',              me_formula,
  'lower bound', score_stat - me_formula,
  'upper bound', score_stat + me_formula
) |> 
  mutate(
    value = round(value, 3)
  )

## # A tibble: 7 × 2
##   result         value
##   <chr>          <dbl>
## 1 y            300    
## 2 n           1250    
## 3 score stat     0.241
## 4 Score SE      15.1  
## 5 Score ME       0.024
## 6 lower bound    0.217
## 7 upper bound    0.264

Our two answers are slight different due to rounding and only searching to the nearest 0.001. If we used enough digits in our “search”, the answer will be identical!

Agresti-Coull

The Agresti-Coull (AC) interval ‘mixes’ the score and Wald intervals together.

How? It keeps the same general form of the Wald interval:

\[\hat{\pi} \pm z_{\alpha/2} \sqrt{\frac{\hat{\pi}(1-\hat{\pi})}{n}}\]

but it updates \(y\) and \(n\):

\[\tilde{y} = y + (z^*)^2/2 \\ \tilde{n} = n + (z^*)^2 \\ \tilde{\pi} = \tilde{y} / \tilde{n}\]

# Calculating the 'new' point estimate
y_tilde <- td_total + (z95^2)/2
n_tilde <-  n + z95^2
pi_tilde <- y_tilde / n_tilde

# The standard error
AG_se <- sqrt(pi_tilde * (1-pi_tilde) / n)

# Lower and upper bounds
AG_ci <- 
  c('lower' = pi_tilde - z95 * AG_se,
    'upper' = pi_tilde + z95 * AG_se
    )

AG_ci

##     lower     upper 
## 0.2170939 0.2644992

Just like the score and Wald interval, it’s very, very similar!

Built-in Functions

If you want one function to rule them all, we can use the binom.confint() function in the binom package.

It works similarly to prop.test(), in that we need to specify x and n. However, we can also specify methods as a way to tell it which method to build the confidence interval:

• ‘asymptotic’ = Wald
• ‘wilson’ = score
• ‘ac’ = Agresti-Coull

library(binom)
binom.confint(
  x = td_total,
  n = n,
  conf.level = 0.95,
  methods = 'asymptotic'
)

##       method   x    n mean     lower     upper
## 1 asymptotic 300 1250 0.24 0.2163242 0.2636758

binom.confint(
  x = td_total,
  n = n,
  conf.level = 0.95,
  methods = 'wilson'
)

##   method   x    n mean     lower     upper
## 1 wilson 300 1250 0.24 0.2171436 0.2644495

binom.confint(
  x = td_total,
  n = n,
  conf.level = 0.95,
  methods = 'ac'
)

##          method   x    n mean     lower     upper
## 1 agresti-coull 300 1250 0.24 0.2171302 0.2644629

You can also give it all three methods in one function by giving methods a vector:

binom.confint(
  x = td_total,
  n = n,
  conf.level = 0.95,
  methods = c('asymptotic', 'wilson', 'ac')
)

##          method   x    n mean     lower     upper
## 1 agresti-coull 300 1250 0.24 0.2171302 0.2644629
## 2    asymptotic 300 1250 0.24 0.2163242 0.2636758
## 3        wilson 300 1250 0.24 0.2171436 0.2644495