[Rd] Updated power.t.test

[Claus Ekstroem]

Greetings,

I've tried to update the power.t.test function to allow for different sample sizes and sample variances in the case of a two-sample t test. I'd gladly update 
the corresponding help page if the code is to replace the current power.t.test function. The modified power.t.test code is included below. 

The changes are as follows:

 - Added three new arguments to the function
   * ratio     : the ratio between group sizes (defaults to 1)
   * sd.ratio  : the ratio between group sd's (defaults to 1)
   * df.method : the method used to calculate the df (defaults to Welch's/Satterthwaite as in the t.test function)

The three new arguments are unly used for the unpaired two-sample case (i.e., when type="two.sample"), and they have default values such that function calls in 
"old code" should work. However, I've changed the output for the group sizes and sd's so they now return a 2-vector with the corresponding group sizes and 
sd's. For example 

# Old code 
> power.t.test(power=.8, delta=1, sd=1)

     Two-sample t test power calculation 

              n = 16.71477
          delta = 1
             sd = 1
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

 NOTE: n is number in *each* group 


# New code
> power.t.test(power=.8, delta=1, sd=1)

     Two-sample t test power calculation 

              n = 16.71477, 16.71477
          delta = 1
             sd = 1, 1
      sig.level = 0.05
          power = 0.8
    alternative = two.sided
      df.method = welch

 NOTE: n is vector of number in each group 


This change of output may break some existing code and is not necessarily a "good idea". Another (possibly better) solution to this would be to return 
both ratios and then let print.power.htest present the result in the correct way. Also, the df.method in the result list could be removed for the one-sample 
and the paired case.

Please comment,

Claus Ekstrøm



power.t.test <-
  function (n = NULL, delta = NULL, sd = 1, sig.level = 0.05, power = NULL,
            ratio = 1, sd.ratio = 1, 
            type = c("two.sample", "one.sample", "paired"),
            alternative = c("two.sided", "one.sided"),
            df.method = c("welch", "classical"),
            strict = FALSE) 
{ 
  type <- match.arg(type)

  if (type == "two.sample") {
    
    if (sum(sapply(list(n, delta, sd, power, sig.level, ratio, sd.ratio), is.null)) !=   1) 
      stop("exactly one of n, delta, sd, power, sig.level, ratio and sd.ratio must be NULL")
    
    if (!is.null(ratio) && ratio <= 0)
      stop("ratio between group sizes must be positive")

    if (!is.null(sd.ratio) && sd.ratio <= 0)
      stop("sd.ratio between group sd's must be positive")
  }
  else {

    if (sum(sapply(list(n, delta, sd, power, sig.level), is.null)) !=   1) 
      stop("exactly one of n, delta, sd, power, and sig.level must be NULL")
  }
  
  alternative <- match.arg(alternative)
  df.method <- match.arg(df.method)
  tsample <- switch(type, one.sample = 1, two.sample = 2, paired = 1)
  tside <- switch(alternative, one.sided = 1, two.sided = 2)
  if (tside == 2 && !is.null(delta)) 
    delta <- abs(delta)
  p.body <- quote({
    nu <- switch(tsample, n-1, switch(df.method, welch=(sd^2/n + (sd*sd.ratio)^2/(n*ratio))^2/((sd^2/n)^2/(n-1) + ((sd*sd.ratio)^2/(ratio*n))^2/(n*ratio-1)), 
classical=(1+ratio)*n-2)) 
    pt(qt(sig.level/tside, nu, lower = FALSE), nu, ncp = switch(tsample, sqrt(n/tsample), sqrt(n/(1+sd.ratio/ratio))) * delta/sd, lower = FALSE)
  })
  if (strict & tside == 2) 
    p.body <- quote({
      nu <- switch(tsample, n-1, switch(df.method, welch=(sd^2/n + (sd*sd.ratio)^2/(n*ratio))^2/((sd^2/n)^2/(n-1) + ((sd*sd.ratio)^2/(ratio*n))^2/(n*ratio-1)), 
classical=(1+ratio)*n-2)) 
      qu <- qt(sig.level/tside, nu, lower = FALSE)
      ncp <- switch(tsample, sqrt(n/tsample), sqrt(n/(1+sd.ratio/ratio))) * delta/sd
      pt(qu, nu, ncp = ncp, lower = FALSE) + 
        pt(-qu, nu, ncp = ncp, lower = TRUE)
    })

  if (is.null(power)) 
    power <- eval(p.body)
  else if (is.null(n)) 
    n <- uniroot(function(n) eval(p.body) - power, c(2, 1e+07))$root
  else if (is.null(sd)) 
    sd <- uniroot(function(sd) eval(p.body) - power, delta * c(1e-07, 1e+07))$root
  else if (is.null(delta)) 
    delta <- uniroot(function(delta) eval(p.body) - power, sd * c(1e-07, 1e+07))$root
  else if (is.null(sig.level)) 
    sig.level <- uniroot(function(sig.level) eval(p.body) - power, c(1e-10, 1 - 1e-10))$root
  else if (is.null(ratio))
    ratio <- uniroot(function(ratio) eval(p.body) - power, c(2/n, 1e+07))$root
  else if (is.null(sd.ratio))
    sd.ratio <- uniroot(function(sd.ratio) eval(p.body) - power, c(1e-07, 1e+07))$root
  else stop("internal error")
  NOTE <- switch(type, paired = "n is number of *pairs*, sd is std.dev. of *differences* within pairs", 
                 two.sample = "n is vector of number in each group", NULL)

  n <- switch(type, paired=n, two.sample=c(n, n*ratio), one.sample=n)
  sd <- switch(type, paired=sd, two.sample=c(sd, sd*sd.ratio), one.sample=sd)
  
  METHOD <- paste(switch(type, one.sample = "One-sample", two.sample = "Two-sample", 
                         paired = "Paired"), "t test power calculation")
  structure(list(n = n, delta = delta, sd = sd, sig.level = sig.level, 
                 power = power, alternative = alternative, note = NOTE, 
                 method = METHOD), class = "power.htest")
}


-- 
*****************************************
Claus Thorn Ekstrøm <ekstrom at dina.kvl.dk>
Dept of Mathematics and Physics, KVL
Thorvaldsensvej 40
DK-1871 Frederiksberg C
Denmark
Phone:[+45] 3528 2341
Fax:  [+45] 3528 2350