Unable to reproduce paper results (sample size)

Question

I am interested in using the results from this paper (Modified exact sample size for a binomial proportion with special emphasis on diagnostic rest parameter estimation by Geoffrey Fosgate) to calculate sample sizes. There is an implementation of the algorithm used to calculate the sample size in R's binomSamSize package. I have checked the code used to calculate this and it seems to follow the paper exactly, but the results shown in the paper's results section do not seem to be reproducible. The sample sizes turn out to be too small.

The paper requires a subscription to be viewed in it's entirety, but here is a snippet to show how the sample size is calculated:

Here are the results I want to reproduce:

Here is the code I am using (lifted from the binomSamSize package):

ciss.midp <- function(p0, d, alpha, nMax = 1e+06){
  pi.L <- p0 - d
  pi.U <- p0 + d
  if (pi.L < 0) 
    stop("p0 - d is below zero!")
  if (pi.U > 1) 
    stop("p0 + d is above one!")
  n <- floor(max(1/p0, 1/(1 - p0)))
  done <- FALSE
  while (!done & (n < nMax)) {
    n <- n + 1
    x <- round(p0 * n)
    lhs2 <- 1/2 * dbinom(x, size = n, prob = pi.L) + 1/2 * 
      dbinom(x, size = n, prob = pi.U) + 
      pbinom(x, size = n, prob = pi.L, lower.tail = FALSE) + 
      pbinom(x - 1, size = n, prob = pi.U)
    if (!is.na(lhs2)) {
      done <- (lhs2 < alpha)
    }
  }
  return(n)
}

What I end up getting to produce the circled column is this:

> sapply(seq(0.5, 0.9, 0.05), function(i) ciss.midp(p0=i, d=0.1, alpha=0.1))
[1] 68 67 65 61 57 50 42 34 23

PS: If there is anything else I can provide to make this easier to answer, please let me know in the comments.

Answer 1

Maybe the author's implementation - despite the description - starts with a large n and works back and stops at the first time we are above alpha, thus returning the last n where the sum was below alpha? At least an implementation of such an algo would be consistent with the results of the paper (test code below).

There appear to be more than one n, where a switch from the sum being above alpha for n and being below alpha for n+1... However, I do think it's better to go from small n and up.

##Go the other way...not my 1st choice though...
ciss.midp2 <- function (p0, d, alpha, nMax = 1e+05) {
    pi.L <- p0 - d
    pi.U <- p0 + d
    if (pi.L < 0)
        stop("p0 - d is below zero!")
    if (pi.U > 1)
        stop("p0 + d is above one!")
    n <- nMax
    done <- FALSE
    while (!done & (n > 0)) {
        n <- n - 1
        x <- round(p0 * n)
        lhs2 <- 1/2 * dbinom(x, size = n, prob = pi.L) +
          pbinom(x, size = n, prob = pi.L, lower.tail = FALSE) +
          (1 - pbinom(x-1, size = n, prob = pi.U, lower.tail=FALSE)) +
          1/2*dbinom(x, size = n, prob = pi.U)

        if (!is.na(lhs2)) {
            done <- (lhs2 > alpha)
        }
    }
    return(n+1)
}

library("binomSamSize")
p.grid <- seq(0.5, 0.9, 0.05)
sapply(p.grid, function(i) ciss.midp(p0=i, d=0.1, alpha=0.1))
[1] 68 67 65 61 57 50 42 34 23
sapply(p.grid, function(i) ciss.midp2(p0=i, d=0.1, alpha=0.1)
[1] 68 67 65 61 57 52 45 39 29