Bioequivalence and Bioavailability Forum

Dear Ben,

❝ ... at the moment I don't see a difference to CVfromCI.

Correct . The bk apply here also if you use the same number of subjects in the (sequence) groups.

❝ For the parallel and the 2x2 crossover design the results are the same as those from nQuery (sometimes it differs, but only +-1, probably due to rounding/higher precision).

Congratulation for the agreement with nQuery. Also I don't have nQuery and I haven't found the formulas they use in the nQuery handbooks available on the net I strongly argue that the differences are more probably due to the different number of degrees of freedom, namely df=n-1 in the so-called 'paired' design. You (having access to nQuery) can test it with my code suggestions below.
BTW: Jerry Hinze's PASS also has a module "Confidence Intervals for Paired Means with Tolerance Probability" with df=n-1.

My checks with your code against the tables in the Joulious book show total agreement for the parallel group design as well as for the 2x2 crossover.

Here some minor code suggestions:
ss_ci <- function(w, sigma, design="2x2", alpha=0.1, gamma=0.05, debug=FALSE)
{# w is the desired precision (CI halfwidth)
 # 1-alpha is the confidence level
 # 1-gamma is the coverage probability (PASS calls this tolerance probability)  
  if (gamma<=0) stop("gamma must lie between >0 and <0.5 !", call. = FALSE)
  if (gamma>0.5) stop("gamma > 0.5 does'nt make sense!", call. = FALSE)
  if (alpha>=1 | alpha<=0) stop("alpha must be 0 < alpha <1 !", call. = FALSE)
  if (sigma<=0 ) stop("sigma must be > 0 !", call. = FALSE)
 
  if (design != "paired") { 
  # get degrees of freedom and design constant
  # (copied from CVfromCI)  
  d.no <- .design.no(design)
  if (is.na(d.no)) stop("Unknown design!", call. = FALSE)
    desi <- .design.props(d.no)
    dfe  <- parse(text = desi$df[1], srcfile = NULL)
    bk   <- desi$bk
    step <- desi$steps
  } else {
    # this should give the results for PASS, nQuery
    # in PASS the sigma of the paired diffs has to be given =sqrt(2)*sigma
    bk <- 2; step <- 1
    dfe  <- parse(text = "n-1", srcfile = NULL)
  }
  n  <- 1
  # DL: while loop replaced
  n <- -eval(dfe) + 2
  # DL addition:
  # normal approximation as starting value
  # gives for all examples checked sampsiz lower than final result
  # but avoids brute force counter from n=3 on
  # this may result in thousands if iterations: check w=0.05, sigma=1 with debug=TRUE
  z  <- qnorm(1-alpha/2)
  n1 <- trunc(bk * z^2*(sigma/w)^2 + z^2/4) #Julious (8.20)
  n  <- max(n,n1)
  # DL: make an multiple of step size (see below)
  n  <- step*trunc(n/step)
  df <- eval(dfe)
  x  <- (w/sigma)^2*eval(dfe)*n/bk/qf(1-alpha, 1, df)
  prob <- pchisq(x,df)
  if (debug) cat("n=",n,"prob=",prob,"\n") # for testing purposes
  while ( prob < 1-gamma ) {
    # DL: since the 'design constants' are derived for equal numbers in the
    # (sequence) groups it would be a fine idea to increment by step size
    n <- n+step # this gives the Julious numbers rounded up to even
    df <- eval(dfe)
    x <- (w/sigma)^2*df*n/bk/qf(1-alpha, 1, df)
    prob <- pchisq(x, df)
    if (debug) cat("n=",n,"prob=",prob,"\n") # for testing purposes
  }
  return(n)
}