Power calculation for two-sample Welch's t test

Question

The R function power.t.test does power calculations (outputs power, sample size, effect size, or whichever parameter you leave out) for t-tests, but only has a single parameter for sample size. The pwr package has a function pwr.t2n.test that performes calculations for a two-sample t-test with different sample sizes (n1,n2). Finally, this suite of stats functions includes a function for Welch's t test (used for samples with different variances), but only includes one parameter for the sample sizes.

I have not been able to find a formula for calculating the power for Welch's t.

Could someone help me out with the formula or an R function for this?

Answer 1

Comment: First, I would suggest you consider carefully whether you have a really good reason to use different sample sizes. Especially if the smaller sample size is used for the group with the larger population, this is not an efficient design.

Second, you can use simulation to get the power for various scenarios. For example, if you use $n_1 = 20,\, \sigma_1 = 15,\,$ $n_2 = 50, \sigma_2 = 10,$ then you have about 75% power for detecting a difference $\delta = 10$ in population means with a Welch test at level 5%.

n1=20; n2=50;  sg1=15; sg2=10;  dlt=10
set.seed(619)
pv = replicate(10^6,
       t.test(rnorm(n1,0,sg1),rnorm(n2,dlt,sg2))$p.val)
mean(pv <= .05)
[1] 0.753043

Because the P-value is taken directly from the procedure t.test in R, results should be accurate to 2 or 3 places, but this style of simulation runs slowly (maybe 2 or 3 min.) with a million iterations.

You might want to use 10,000 iterations if you are doing repeated runs for various sample sizes, and then use a larger number of iterations to verify the power of the final design.

Changing to $n_2 = 20$ gives power 67%, so the extra 30 observations in Group 2 are not 'buying' you as much as you might hope. By contrast, a balanced design with $n_1 = n_2 = 35$ gives about 90% power (with everything else the same).

Answer 2

The article "Optimal sample sizes for Welch’s test under various allocation and cost considerations" from Show-Li Jan & Gwowen Shieh published in Behavior Research Methods in December 2011 has the following code in supplementary material A, slightly modified here for my own ergonomy.

ssize.welch = function(alpha=0.05, power=0.90, mu1, mu2, sigma1, sigma2, n2n1r, use_exact=FALSE)
 {
  mud<-mu1-mu2
  sigsq1<-sigma1^2
  sigsq2<-sigma2^2

  numint<-1000
  dd<-0.00001
  coevec<-c(1,rep(c(4,2),numint/2-1),4,1)
  intl<-(1-2*dd)/numint
  bvec<-intl*seq(0,numint)+dd

  #Z method
  za<-qnorm(1-alpha/2)
  zb<-qnorm(power)
  n1<-ceiling(((sigsq1+sigsq2/n2n1r)*(za+zb)^2)/(mud^2))
  n2<-ceiling(n1*n2n1r)

  if(use_exact) #Exact method
   {
    n1 = n1-1
    n2 = n2-1
    powere<-0

    while(powere<power)
     {
      n1<-n1+1
      n2<-ceiling(n1*n2n1r)
      sigsqt<-sigsq1/n1+sigsq2/n2
      hsigsqt<-sqrt(sigsqt)
      wpdf<-(intl/3)*coevec*dbeta(bvec,(n1-1)/2,(n2-1)/2)
      dft<-n1+n2-2 
      p1<-(n1-1)/dft
      p2<-1-p1
      s1<-sigsq1/n1
      s2<-sigsq2/n2
      b12<-(s1/p1)*bvec+(s2/p2)*(1-bvec)
      r1<-(s1/p1)*bvec/b12
      r2<-1-r1
      dfevec<-1/((r1^2)/(n1-1)+(r2^2)/(n2-1))
      tdfea<-qt(1-alpha/2,dfevec)
      powere<-sum(wpdf*pt(-tdfea*sqrt(b12/sigsqt),dft,mud/hsigsqt))+1-sum(wpdf*pt(tdfea*sqrt(b12/sigsqt),dft,mud/hsigsqt))
      }
    }

  c(n1=n1,n2=n2)
  }

ssize.welch(0.05,0.9,85,105,10,20,3)
ssize.welch(0.05,0.9,85,105,10,20,3,TRUE)

The Z method is also the one used in https://clincalc.com/stats/samplesize.aspx which cites Rosner B. Fundamentals of Biostatistics. 7th ed. Boston, MA: Brooks/Cole; 2011. Weirdly, it forces you to input only one variance but the formula it gives can use two (and is the same as the paper above). It is in the spirit of the usual way of computing sample sizes, but I'm not too sure about when the underlying approximations start to be false.

After some fiddling around, both methods give similar results unless you have very unbalanced groups (but in that case, at some point you might want to just approximate the super large group as a known population).

Hope this helps.

Edit : just realized this does not exactly answer your question about power rather than sample size, but you can easily flip the Z method formula to compute power (exact method seems more hairy ; worst case, numeric trial and error should work since the relationship is monotonous).

Answer 3

It seems to me the power calculation is pretty straight forward. Let's consider a 2 tailed alternative.

The null is rejected if $t<t_{\alpha/2}$ or $t>t_{1-\alpha/2}$. So Power = $P(T'<t_{t_\alpha/2})+P(T'>t_{1-\alpha/2})$

The true effect is $\eta=\frac{\mu_1-\mu_2}{\sqrt{s_1^2/n_1+s_2^2/n_2}}$

The true sampling distribution of $T'$ is a noncentral T distribution centered at $\eta$. Let $T = T'-\eta$ follow a t distribution centered at zero. So

$$\text{power} = P(T<t_{t_\alpha/2}-\eta)+P(T>t_{1-\alpha/2}-\eta)$$

I just figured this out on my own, so I may have made a mistake, but It seems to make sense to me. This formula cannot be reversed to give you the necessary sample size though, it would have to be solved by numerical means to get a pre-specified power given effect size (or given $\mu_1, \mu_2, \sigma_1, \sigma_2$. It would be easiest to solve it for $n_1=n_2=n$, but in theory you could specify (for example) $n_1$ and solve for $n_2$.