Power for two sample t test

Question

I am trying to understand power calculation for the case of the two independent sample t-test (not assuming equal variances so I used Satterthwaite).

Here is a diagram that I found to help understand the process:

So I assumed that given the following about the two populations and given the sample sizes:

mu1<-5
mu2<-6
sd1<-3
sd2<-2
n1<-20
n2<-20

I could compute the critical value under the null relating to having 0.05 upper tail probability:

df<-(((sd1^2/n1)+(sd2^2/n2)^2)^2) / ( ((sd1^2/n1)^2)/(n1-1) + ((sd2^2/n2)^2)/(n2-1)  )
CV<- qt(0.95,df) #equals 1.730018

and then calculate the alternative hypothesis (which for this case I learned is a "non central t distribution"). I calculated beta in the diagram above using the non central distribution and the critical value found above. Here is the full script in R:

#under alternative
mu1<-5
mu2<-6
sd1<-3
sd2<-2
n1<-20
n2<-20


#Under null
Sp<-sqrt(((n1-1)*sd1^2+(n2-1)*sd2^2)/(n1+n2-2))
df<-(((sd1^2/n1)+(sd2^2/n2)^2)^2) / ( ((sd1^2/n1)^2)/(n1-1) + ((sd2^2/n2)^2)/(n2-1)  )
CV<- qt(0.95,df)


#under alternative
diff<-mu1-mu2
t<-(diff)/sqrt((sd1^2/n1)+ (sd2^2/n2))
ncp<-(diff/sqrt((sd1^2/n1)+(sd2^2/n2)))


#power
1-pt(t, df, ncp)

This gives a power value of 0.4935132.

Is this the correct approach? I find that if I use other power calculation software (like SAS, which I think I have set up equivalently to my problem below) I get another answer (from SAS it is 0.33).

SAS CODE:

proc power;
      twosamplemeans test=diff_satt
         meandiff = 1
         groupstddevs = 3 | 2
         groupweights = (1 1)
         ntotal = 40
         power = .
        sides=1;
   run;

Ultimately, I would like to get an understanding that would allow me to look at simulations for more complicated procedures.

EDIT: I found my error. should have been

1-pt(CV, df, ncp) NOT 1-pt(t, df, ncp)

Answer 1

You're close, some small changes are required though:

• The true difference in means is typically taken as $\mu_{2} - \mu_{1}$, not the other way around.
• G*Power uses $n_{1}+n_{2} - 2$ as degrees of freedom for the $t$-distribution in this case (different variances, same group sizes), following a suggestion from Cohen as explained here
• SAS might use Welch's formula or Satterthwaite's formula for the df given unequal variances (found in this pdf you cited) - with only 2 significant digits in the result one cannot tell (see below)

With n1, n2, mu1, mu2, sd1, sd2 as defined in your question:

> alpha   <- 0.05
> dfGP    <- n1+n2 - 2                     # degrees of freedom (used by G*Power)
> cvGP    <- qt(1-alpha, dfGP)             # crit. value for one-sided test (under the null)
> muDiff  <- mu2-mu1                       # true difference in means
> sigDiff <- sqrt((sd1^2/n1) + (sd2^2/n2)) # true SD for difference in empirical means
> ncp     <- muDiff / sigDiff              # noncentrality parameter (under alternative)
> 1-pt(cvGP, dfGP, ncp)                    # power
[1] 0.3348385

This matches the result from G*Power which is a great program for these questions. It displays the df, critical value, ncp as well, so you can check all these calculations separately.

Edit: Using Satterthwaite's formula or Welch's formula doesn't change much (still 0.33*):

# Satterthwaite's formula
> var1  <- sd1^2
> var2  <- sd2^2
> num   <- (var1/n1 + var2/n2)^2
> denST <- var1^2/((n1-1)*n1^2) + var2^2/((n2-1)*n2^2)
> (dfST <- num/denST)
[1] 33.10309

> cvST <- qt(1-alpha, dfST)
> 1-pt(cvST, dfST, ncp)
[1] 0.3336495

# Welch's formula
> denW <- var1^2/((n1+1)*n1^2) + var2^2/((n2+1)*n2^2)
> (dfW <- (num/denW) - 2)
[1] 34.58763

> cvW   <- qt(1-alpha, dfW)
> 1-pt(cvW, dfW, ncp)
[1] 0.3340453

(note that I slightly changed some variable names as t, df, and diff are also names of built-in functions, also note that the numerator of your code for df is wrong, it has a misplaced ^2, and one ^2 too many, it should be ((sd1^2/n1) + (sd2^2/n2))^2)

Answer 2

If you are mainly interested in computing the power (rather than learning through doing it by hand) and you are already using R then look at the pwr package and either the pwr.t.test or pwr.t2n.test functions. (these can be good to verify your results even if you do it by hand to learn).