Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

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:

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)

share|improve this question
add comment

2 Answers

up vote 6 down vote accepted

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)

share|improve this answer
    
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. –  B_Miner Oct 15 '11 at 21:51
    
I found my error and adjusted above in my question. Thanks again! –  B_Miner Oct 15 '11 at 22:16
1  
@B_Miner I've updated my answer to address your question. –  caracal Oct 15 '11 at 22:47
add comment

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

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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