Confidence Intervals for a Proportion in R

Introduction

[See separate page for general introductory information on confidence intervals for proportions.]

[Note: information about cicalc package will be added to this page soon.]

Data used

The adcibc data stored here was used in this example, creating a binary treatment variable trt taking the values of ACT or PBO and a binary response variable resp taking the values of Yes or No. For this example, a response is defined as a score greater than 4.

The below shows that for the Actual Treatment, there are 36 responders out of 154 subjects = 0.234 (23.4% responders).

# A tibble: 4 × 3
# Groups:   trt [2]
  trt   resp      n
  <chr> <chr> <int>
1 ACT   No      118
2 ACT   Yes      36
3 PBO   No       65
4 PBO   Yes      12

Packages

The {cardx} package is an extension of the {cards} package, providing additional functions to create Analysis Results Data Objects (ARDs)¹. It was developed as part of {NEST} and pharmaverse. This package requires the binary endpoint to be a logical (TRUE/FALSE) vector or a numeric/integer coded as (0, 1) with 1 (TRUE) being the success you want to calculate the confidence interval for.

See here for full description of the {cardx} proportions equations.

If calculating the CI for a difference in proportions, the package requires both the response and the treatment variable to be numeric/integer coded as (0, 1) (or logical vector).

Instead of the code presented below, you can use ard_categorical_ci(data, variables=resp, method ='wilson') for example. This invokes the code below but returns an analysis results dataset (ARD) format as the output. Methods included are waldcc, wald, clopper-pearson, wilson, wilsoncc, strat_wilson, strat_wilsoncc, agresti-coull and jeffreys for one-sample proportions and methods for 2 independent samples, however currently does not have a method for 2 matched proportions.

The {ratesci} package is … [TBC] - note, current development version at https://github.com/petelaud/ratesci has new features in the rateci() function (including more of the CI methods described below) compared to the CRAN release.

The {PropCIs} package produces CIs for methods such as Blaker’s exact method and Midp which aren’t available in {cardx} but are available in SAS. We found results agreed with SAS to the 5th decimal place. The package also calculates CIs for Clopper-Pearson, Wald, Wilson, Agresti-Coull and these align to results obtained in cardx to at least the 7th decimal place. The {PropsCIs} package requires just the number of events (numerator number of successes) & total number of subjects (denominator) as an input dataset. Given Blaker and Midp are rarely used in practice, and {PropsCIs} isn’t a package commonly downloaded from CRAN, further details are not provided here.

The {Hmisc} package produces CIs using the Clopper-Pearson method. In this example (x=36 and n=154), the results match the cardx package. Documentation reports that the method uses F distribution to compute exact intervals based on the binomial cdf. However, if the percentage of responders is 100% then the upper limit is set to 1. Similarly if the percentage of responders is 0%, then the lower limit is set to 0. Hence, in extreme cases there may be differences between this package and the standard implementation of Clopper-Pearson method.

The {RBesT} package (Prior to Version 1.8-0) produces CIs using the Clopper-Pearson method. In this example (x=36 and n=154), the results match the cardx package. However, as described below, there are 2 cases where the results using RBesT package do not match cardx or Hmisc.

x = 0 (0% responders), in which case the lower limit does not match.
x = n (100% responders), in which case the upper limit does not match.

Because of the relationship between the binomial distribution and the beta distribution. This package uses quantiles of the beta distribution to derive exact confidence intervals.

\[ B(\alpha/2;x, n-x+1) < p < B(1-\alpha/2; x+1, n-x)\]

RBesT equations are:
pLow <- qbeta(Low, r + (r == 0), n - r + 1)
pHigh <- qbeta(High, r + 1, n - r + ((n - r) == 0))

In Version 1.8-0 onwards the equations were updated as follows, which then match the Hmisc intervals:
pLow <- qbeta(Low, r, n - r + 1)
pHigh <- qbeta(High, r + 1, n - r)

The {ExactCIdiff} package produces exact CIs for two dependent proportions (matched pairs).

The {DescTools} package has a function BinomDiffCI which produces CIs for two independent proportions (unmatched pairs) including methods for Agresti/Caffo, Wald, Wald with Continuity correction, Newcombe Score, Newcombe score with continuity correction, and more computationally intensive methods such as Miettinen and Nurminen, Mee, Brown Li’s Jeffreys, Hauck-Anderson and Haldane. See here for more detail.

The {presize} package has a function prec_prop() which also calculates CIs for 2 independent samples using the Wilson, Agresti-Coull, Exact or Wald approaches. The package is not described in further detail here since in most cases {DescTools} will be able to compute what is needed. However, it’s mentioned due to other functionality it has available such as sample size and precision calculations for AUC, correlations, cronbach’s alpha, intraclass correlation, Cohen’s kappa, likelihood ratios, means, mean differences, odds ratios, rates, rate ratios, risk differences and risk ratios.

Methods for Calculating Confidence Intervals for a single proportion using cardx

For more technical derivation and reasons for use of each of the methods listed below, see the corresponding SAS page.

Let’s start by calculating a Confidence interval for the proportion of successes observed in the Active Treatment group (a single sample).

Clopper-Pearson (Exact or binomial CI) Method

Clopper-Pearson Exact CI is one of the most popular methods, it is often good for small sample sizes when the proportion is not close to the tails (0,1), but it can be too conservative (too wide an interval compared to the interval containing the true population proportion 95% of the time).

The cardx package calculates the Clopper-Pearson score by calling stats::binom.test() function.

cardx::proportion_ci_clopper_pearson(act2, conf.level = 0.95) |>
  as_tibble()

# A tibble: 1 × 11
      N     n conf.level estimate statistic  p.value parameter conf.low
  <int> <dbl>      <dbl>    <dbl>     <dbl>    <dbl>     <dbl>    <dbl>
1   154    36       0.95    0.234        36 2.21e-11       154    0.169
# ℹ 3 more variables: conf.high <dbl>, method <chr>, alternative <chr>

Normal Approximation Method (Also known as the Wald Method)

In large random samples from independent trials, the sampling distribution of proportions approximately follows the normal distribution. The expectation of a sample proportion is the corresponding population proportion. Therefore, based on a sample of size \(n\), a \((1-\alpha)\%\) confidence interval for population proportion can be calculated using normal approximation as follows:

\(p\approx \hat p \pm z_\alpha \sqrt{\hat p(1-\hat p)}/{n}\), where \(\hat p\) is the sample proportion, \(z_\alpha\) is the \(1-\alpha/2\) quantile of a standard normal distribution corresponding to level \(\alpha\), and \(\sqrt{\hat p(1-\hat p)}/{n}\) is the standard error.

For more technical information see the corresponding SAS page.

Example code

The following code calculates a confidence interval for a binomial proportion using normal approximation equation manually. This is replicated exactly using the cardx::proportion_ci_wald function which also allows the continuity correction to be applied.

# sample proportion by trt
summary <- adcibc |>
  filter(trt == "ACT") |>
  group_by(resp) |>
  tally() |>
  ungroup() |>
  mutate(
    total = sum(n),
    p = n / total
  )

# Calculate standard error and 95% wald confidence intervals for population proportion
waldci <- summary |>
  filter(resp == "Yes") |>
  mutate(
    se = sqrt(p * (1 - p) / total),
    lower_ci = (p - qnorm(1 - 0.05 / 2) * se),
    upper_ci = (p + qnorm(1 - 0.05 / 2) * se)
  )
waldci

# A tibble: 1 × 7
  resp      n total     p     se lower_ci upper_ci
  <chr> <int> <int> <dbl>  <dbl>    <dbl>    <dbl>
1 Yes      36   154 0.234 0.0341    0.167    0.301

# cardx package Wald method without continuity correction
cardx::proportion_ci_wald(act2, conf.level = 0.95, correct = FALSE) |>
  as_tibble()

# A tibble: 1 × 7
      N     n estimate conf.low conf.high conf.level method                     
  <int> <dbl>    <dbl>    <dbl>     <dbl>      <dbl> <glue>                     
1   154    36    0.234    0.167     0.301       0.95 Wald Confidence Interval w…

# cardx package Wald method with continuity correction
cardx::proportion_ci_wald(act2, conf.level = 0.95, correct = TRUE) |>
  as_tibble()

# A tibble: 1 × 7
      N     n estimate conf.low conf.high conf.level method                     
  <int> <dbl>    <dbl>    <dbl>     <dbl>      <dbl> <glue>                     
1   154    36    0.234    0.164     0.304       0.95 Wald Confidence Interval w…

Wilson Method (Also known as the Score method or the Altman, Newcombe method³ )

The cardx package calculates the Wilson (score) method by calling stats::prop.test() function. This method is often used as a compromise between the Clopper-Pearson and the Wald given it was found to be accurate for most parameter values (even those close to 0 and 1), and it does not suffer from being over-conservative. For more technical information see the corresponding SAS page.

The package also contains a function for proportion_ci_strat_wilson() which calculates the stratified Wilson CIs for unequal proportions as described on page 47 here.

# cardx package Wilson method without continuity correction
cardx::proportion_ci_wilson(act2, conf.level = 0.95, correct = FALSE) |>
  as_tibble()

# A tibble: 1 × 11
      N     n conf.level estimate statistic  p.value parameter conf.low
  <int> <dbl>      <dbl>    <dbl>     <dbl>    <dbl>     <int>    <dbl>
1   154    36       0.95    0.234      43.7 3.90e-11         1    0.174
# ℹ 3 more variables: conf.high <dbl>, method <glue>, alternative <chr>

# cardx package Wilson method with continuity correction
cardx::proportion_ci_wilson(act2, conf.level = 0.95, correct = TRUE) |>
  as_tibble()

# A tibble: 1 × 11
      N     n conf.level estimate statistic  p.value parameter conf.low
  <int> <dbl>      <dbl>    <dbl>     <dbl>    <dbl>     <int>    <dbl>
1   154    36       0.95    0.234      42.6 6.70e-11         1    0.171
# ℹ 3 more variables: conf.high <dbl>, method <glue>, alternative <chr>

Agresti-Coull Method

The cardx package calculates the Agresti-Coull method using the equation from the published method by Alan Agresti & Brent Coull based on adding 2 successes and 2 failures before computing the wald CI. The CI is truncated, when it overshoots the boundary (<0 or >1).

# cardx package agresti_coull method
cardx::proportion_ci_agresti_coull(act2, conf.level = 0.95) |>
  as_tibble()

# A tibble: 1 × 7
      N     n estimate conf.low conf.high conf.level method                     
  <int> <dbl>    <dbl>    <dbl>     <dbl>      <dbl> <chr>                      
1   154    36    0.234    0.174     0.307       0.95 Agresti-Coull Confidence I…

Jeffreys Method

Jeffreys method is a particular type of Bayesian Highest Probability Density (HPD) Method. For proportions, the beta distribution is generally used for the prior, which consists of two parameters alpha and beta. Setting alpha=beta=0.5 is called Jeffrey’s prior. NOTE: if you want to use any other priors, you can use binom.bayes which estimates a credible interval for proportions.

# cardx package jeffreys method
cardx::proportion_ci_jeffreys(act2, conf.level = 0.95) |>
  as_tibble()

# A tibble: 1 × 7
      N     n estimate conf.low conf.high conf.level method           
  <int> <dbl>    <dbl>    <dbl>     <dbl>      <dbl> <glue>           
1   154    36    0.234    0.172     0.305       0.95 Jeffreys Interval

Continuity Adjusted Methods

[TBC]

Consistency with hypothesis tests