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:

enter image description here

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


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

#Under null
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
t<-(diff)/sqrt((sd1^2/n1)+ (sd2^2/n2))

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).


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

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)


2 Answers 2


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.

enter image description here

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)

  • 1
    $\begingroup$ thanks! The one thing is, doesn't this formula for df assume that the population standard deviations are equal? See page 3 of the following (where I got the Satterthwaite df): stata-journal.com/sjpdf.html?articlenum=st0062. Supposedly, SAS uses this approximation in the proc I posted. $\endgroup$
    – B_Miner
    Oct 15, 2011 at 21:51
  • 1
    $\begingroup$ I found my error and adjusted above in my question. Thanks again! $\endgroup$
    – B_Miner
    Oct 15, 2011 at 22:16
  • 2
    $\begingroup$ @B_Miner I've updated my answer to address your question. $\endgroup$
    – caracal
    Oct 15, 2011 at 22:47

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).


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.