0

I was wondering how I could check via simulation in R that the 95% Confidence Interval obtained from a binomial test with 5 successes in 15 trials when TRUE p = .5 has a 95% "Coverage Probability" in the long-run?

Here is the 95% CI for such a test using R (how can show that the following CI has a 95% coverage in the long-run if TRUE p = .5):

as.numeric(binom.test(x = 5, n = 15, p = .5)[[4]])
# > [1] 0.1182411 0.6161963 (in the long-run 95% of the time, ".5" is contained within these
#                            two numbers, how to show this in R?)
rnorouzian
  • 7,397
  • 5
  • 27
  • 72
  • Yes, I understand that, but there are several formulas for the binomial proportion's CI. Which one do you want to use? – Rui Barradas Aug 08 '17 at 15:57
  • @RuiBarradas, does that now sound clearer what I want to achieve? – rnorouzian Aug 08 '17 at 16:16
  • @RuiBarradas, I guess your comments seem irrelavant to OP's question! OP just wants to know the coverage probability of the current method for CI in R in the long-run. Please delete your comments! –  Aug 08 '17 at 16:22
  • FWIW, binom.test uses Clopper-Pearson intervals, which overcover: see [this figure](http://imgur.com/a/N1Ff5) from Agresti's *Categorical Analysis* book – Ben Bolker Aug 08 '17 at 23:41

1 Answers1

1

Something like this?

fun <- function(n = 15, p = 0.5){
    x <- rbinom(1, size = n, prob = p)
    res <- binom.test(x, n, p)[[4]]
    c(Lower = res[1], Upper = res[2])
}

set.seed(3183)
R <- 10000
sim <- t(replicate(R, fun()))

Note that binom.test when called with 5 successes, 15 trials and p = 0.5 will always return the same value, hence the call to rbinom. The number of successes will vary. We can compute the proportion of cases when p is between Lower and Upper.

cov <- mean(sim[,1] <= .5 & .5 <= sim[,2])
rnorouzian
  • 7,397
  • 5
  • 27
  • 72
Rui Barradas
  • 70,273
  • 8
  • 34
  • 66