SSiz MulMeans

Content Disclaimer
Copyright @2020.
All Rights Reserved.

StatsToDo: Sample Size for Comparing Means from multiple Groups (One Way Analysis of Variance)

Links : Home Index (Subjects) Contact StatsToDo

Explanations and References Sample Size Tables Javascript Programs

Data

Sample Size is a table of 5 columns
  - Each row contains data from a separate study
  - Col 1 = Probability of Type I Error (α)
  - Col 2 = Power (1-β)
  - Col 3 = The Number of Groups in the study
  - Col 4 = Maximum Difference Between Any Two Groups to be Detected
  - Col 5 = Expected within group Standard Deviation

Power Estimation is a table of 5 columns
  - Each row contains data from a separate study
  - Col 1 = Probability of Type I Error α
  - Col 2 = Averaged numbers per group used
  - Col 3 = The Number of Groups in the study
  - Col 4 = Maximum Difference Between Any Two Groups observed
  - Col 5 = Observed within group SD (sqrt(observed esidual variance))

Pilot is for a single plan, a single column with 5 rows
  - Row 1 : Probability of Type I Error usually 0.05
  - Row 2 : The Number of Groups in the study
  - Row 3 : Largest Difference Between Groups to be detected
  - Row 4 : Expected within group or population Standard Deviation
  - Row 5 : Sample size interval for tabulation of changing values of Power
  - Row 6 : Maximum sample size for tabulation of changing values of Power

R Codes

This panel provides R codes for sample size related to comparing means in 2 or more grops, the One Way Analysis of Variance (OWAV)

Ref: Cohen J (1988) Statistical power analysis for the behavioral sciences. Second edition. Lawrence Erlbaum Associates, Pubishers. London. ISBN 0-8058-0283-5 p. 276, 550

Subroutine functions used by all 3 programs

# Subroutine calculating the 3 Effect sizes 
# f1: Model 1. Minimal difference between means
# f2: Model 2. Differences between means evenly distributed
# f3: Model 3. All differences between means close to declared maximum
CalEffectSize <- function(k, diff, sd)   # nroup, average diff, sd
{
  d = diff / sd                                # eq 8.2.5. p. 276
  f1 = d * sqrt(1 / (2 * k))                   # eq 8.2.8. p. 277
  f2 = d / 2 * sqrt((k + 1) / (3 * (k - 1)))   # 8.2.10 p. 278
  f3 = d / 2                                   # 8.2.12 p. 279
  if(k %% 2 != 0)
  {
    f3 = d * ((sqrt(k**2 - 1) / (2 * k)))      # 8.2.13 p. 279
  }
  return (c(f1, f2, f3))
} 

# Subroutine estimating beta, probability of Type II Error
CalBeta <- function(alpha, nn, k, f)   # alpha,ssiz,ngrp,f
{
  n = round(nn)
  u = (k - 1)
  v = k * (n + 1)                      # denominator equation 12.8.3
  lambda = f^2 * n * k                 # Lambda equation 12.8.2 p. 550
  Fc = qf(1 - alpha, u, v)             # F for alpha
  # Laubscher's equation
  top1 = sqrt(2.0 * (u + lambda) - (u + 2 * lambda) / (u + lambda))
  top2 = sqrt((2.0 * v - 1) * u * Fc / v)
  bot = sqrt(u * Fc / v + (u + 2 * lambda) / (u + lambda)) 
  zb = (top1 - top2) / bot;
  return (1-pnorm(zb))
}

Program 1: Sample size

The program calculates sample size repeatedly, so the actual calculations are in a subroutine, which is called from the main program

  
# pgm 1: Sample size

# subroutine to calculate sample size.  It needs CalEffectSize and CalBeta
# the algorithm iterates with sample size until the required beta is reached
CalSSiz <- function(alpha, beta, k, f) # k= number of groups, f is as estimated by CalEffectSize
{ 
  nL = 20                                
  nR = 1e5                              
  pL = CalBeta(alpha, nL, k, f) - beta
  pR = CalBeta(alpha, nR, k, f) - beta
  oldnM = nR
  nM = (nL + nR) / 2.0;  
  pM = CalBeta(alpha, nM, k, f) - beta;
  i = 0;
  while(abs(oldnM - nM)>1 & i<20 & abs(pM)>0.0001)
  {
    oldnM = nM
    if(pL * pM > 0) 
    { 
      nL = nM
    }
    else 
    { 
      nR = nM 
    } 
    nM = (nL + nR) / 2.0 
    pL = CalBeta(alpha, nL, k, f) - beta
    pM = CalBeta(alpha, nM, k, f) - beta
    #print(c(oldnM, nM, pM))
    i = i + 1
  }
  return (ceiling(nM))
}

# Main program: needs subroutines CalEffectSize, CalSSiz
# data entry
dat = ("
Alpha  Power NGrp      Diff      SD
0.05	0.8	2	0.5	1.0
0.05	0.8	3	0.5	1.0
0.05	0.8	4	0.5	1.0
0.01	0.8	2	0.5	1.0
0.01	0.8	3	0.5	1.0
0.01	0.8	4	0.5	1.0
       ")
df <- read.table(textConnection(dat),header=TRUE)  # conversion to data frame    

# vector for results
F1 <- vector()
SSiz1 <- vector()
F2 <- vector()
SSiz2 <- vector()
F3 <- vector()
SSiz3 <- vector()

# Calculate
for(i in 1 : nrow(df))
{
  alpha = df$Alpha[i]
  beta = 1 - df$Power[i]
  ngrp = df$NGrp[i]
  diff = df$Diff[i]
  sd = df$SD[i]
  arES = CalEffectSize(ngrp, diff, sd)          # calculate effect sizes f1, f2, f3
  f1 = arES[1]                                  # ssiz by f1
  ssiz1 = CalSSiz(alpha, beta, ngrp, f1)
  F1 <- append(F1, sprintf(f1, fmt="%#.4f"))
  SSiz1 <- append(SSiz1, ssiz1)
  f2 = arES[2]                                  # ssiz by f2
  ssiz2 = CalSSiz(alpha, beta, ngrp, f2)
  F2 <- append(F2, sprintf(f2, fmt="%#.4f"))
  SSiz2 <- append(SSiz2, ssiz2)                 # ssiz by f3  
  f3 = arES[3]
  ssiz3 = CalSSiz(alpha, beta, ngrp, f3)
  F3 <- append(F3, sprintf(f3, fmt="%#.4f"))
  SSiz3 <- append(SSiz3, ssiz3)
}
# combine results in data frame and show
df$F1 <- F1
df$SSiz1 <- SSiz1
df$F2 <- F2
df$SSiz2 <- SSiz2
df$F3 <- F3
df$SSiz3 <- SSiz3
df # data frame with input data and results

The results are as follows

F1, F2, and F3 are effect sizes of the 3 models
SSiz1, SSiz2, and SSiz3 are the sample size (per group) according to the 3 models

> df # data frame with input data and results
  Alpha Power NGrp Diff SD     F1 SSiz1     F2 SSiz2     F3 SSiz3
1  0.05   0.8    2  0.5  1 0.2500    62 0.2500    62 0.2500    62
2  0.05   0.8    3  0.5  1 0.2041    78 0.2041    78 0.2357    58
3  0.05   0.8    4  0.5  1 0.1768    88 0.1863    79 0.2500    46
4  0.01   0.8    2  0.5  1 0.2500    94 0.2500    94 0.2500    94
5  0.01   0.8    3  0.5  1 0.2041   111 0.2041   111 0.2357    85
6  0.01   0.8    4  0.5  1 0.1768   125 0.1863   111 0.2500    62

Program 2: Power

# Pgm 2: Power: needs subroutines CalEffectSize, CalBeta
# data entry
dat = ("
Alpha SSiz   NGrp      Diff      SD
0.05	62	2	0.5	1.0
0.05	58	3	0.5	1.0
0.05	44	4	0.5	1.0
0.01	94	2	0.5	1.0
0.01	86	3	0.5	1.0
0.01	63	4	0.5	1.0
       ")
df <- read.table(textConnection(dat),header=TRUE)  # conversion to data frame    

# result vectors
F1 <- vector()
Power1 <- vector()
F2 <- vector()
Power2 <- vector()
F3 <- vector()
Power3 <- vector()
# calculations
for(i in 1 : nrow(df))
{
  alpha = df$Alpha[i]
  ssiz = df$SSiz[i]
  ng = df$NGrp[i]
  diff = df$Diff[i]
  sd = df$SD[i]
  arES = CalEffectSize(ng, diff, sd)                     # establish the 3 effect sizes f1, f2, f3
  f1 = arES[1]                                           # power according to f1
  power1 = 1 - CalBeta(alpha, ssiz, ng, f1)
  F1 <- append(F1, sprintf(f1, fmt="%#.4f"))
  Power1 <- append(Power1, sprintf(power1, fmt="%#.4f"))
  f2 = arES[2]                                           # power according to f2
  power2 = 1 - CalBeta(alpha, ssiz, ng, f2)
  F2 <- append(F2, sprintf(f2, fmt="%#.4f"))
  Power2 <- append(Power2, sprintf(power2, fmt="%#.4f"))
  f3 = arES[3]                                           # power according to f3 
  power3 = 1 - CalBeta(alpha, ssiz, ng, f3)
  F3 <- append(F3, sprintf(f3, fmt="%#.4f"))
  Power3 <- append(Power3, sprintf(power3, fmt="%#.4f"))
}  
# append results to data frame and display
df$F1 <- F1
df$Power1 <- Power1
df$F2 <- F2
df$Power2 <- Power2
df$F3 <- F3
df$Power3 <- Power3
df # show input data and results

The results are as follow.

F1, F2, and F3 are the effect sizes of the 3 models
Power1, Power2, and Power3 are the corresponding power (1 - β)

> df # show input data and results
  Alpha SSiz NGrp Diff SD     F1 Power1     F2 Power2     F3 Power3
1  0.05   62    2  0.5  1 0.2500 0.7986 0.2500 0.7986 0.2500 0.7986
2  0.05   58    3  0.5  1 0.2041 0.6725 0.2041 0.6725 0.2357 0.8005
3  0.05   44    4  0.5  1 0.1768 0.4794 0.1863 0.5252 0.2500 0.7992
4  0.01   94    2  0.5  1 0.2500 0.8007 0.2500 0.8007 0.2500 0.8007
5  0.01   86    3  0.5  1 0.2041 0.6526 0.2041 0.6526 0.2357 0.8117
6  0.01   63    4  0.5  1 0.1768 0.4014 0.1863 0.4575 0.2500 0.8018

Program 3: Pilot Study

The program calculates power with increasing sample size

# Pgm 3. Pilot study: needs subroutines CalEffectSize, CalBeta

# Parameters
alpha = 0.05         # Probability of Type I Error
ngrp = 3             # number of group
diff = 0.5           # largest difference between any 2 groups
sd = 1.0             # within group (population) Standard Deviation
intv = 5             # interval of sample size
maxN = 100           # maximum sample size

# estimating the effect sizes f1, f2, f3 for the 3 models
arES <- CalEffectSize(ngrp, diff, sd)   # nroup, average diff, sd
f1 = arES[1]  
f2 = arES[2]  
f3 = arES[3]
  
# vectors for results
SSiz <- vector()     # sample size
Power1 <- vector()   # Power for Effect Size 1
Diff1 <- vector()    # difference in power1 from previous row
IncCs1 <- vector()   # increase in power1 per case increase
PICs1 <- vector()    # % increase in power1 per case increase
Power2 <- vector()   # Power for Effect Size 2
Diff2 <- vector()    # difference in power2 from previous row
IncCs2 <- vector()   # increase in power2 per case increase
PICs2 <- vector()    # % increase in power2 per case increase
Power3 <- vector()   # Power for Effect Size 3
Diff3 <- vector()    # difference in power3 from previous row
IncCs3 <- vector()   # increase in power3 per case increase
PICs3 <- vector()    # % increase in power3 per case increase

# Calculations
# first row
n = intv                                                     # sample size
SSiz <- append(SSiz,n)
power1 = 1 - CalBeta(alpha, n, ngrp, f1)                     # power 1   
Power1 <- append(Power1, sprintf(power1, fmt="%#.4f"))  
Diff1 <- append(Diff1,0)                                     # difference in power1 from previous row 
IncCs1 <- append(IncCs1,0)                                   # increase in power1 per case increase
PICs1 <- append(PICs1,0)                                     # % increase in power1 per case increase
power2 = 1 - CalBeta(alpha, n, ngrp, f2)                     # Power 2
Power2 <- append(Power2, sprintf(power2, fmt="%#.4f"))  
Diff2 <- append(Diff2,0)                                     # difference in power2 from previous row
IncCs2 <- append(IncCs2,0)                                   # increase in power2 per case increase
PICs2 <- append(PICs2,0)                                     # % increase in power2 per case increase
power3 = 1 - CalBeta(alpha, n, ngrp, f3)                     # power 3   
Power3 <- append(Power3, sprintf(power3, fmt="%#.4f"))  
Diff3 <- append(Diff3,0)                                     # difference in power3 from previous row
IncCs3 <- append(IncCs3,0)                                   # increase in power3 per case increase
PICs3 <- append(PICs3,0)                                     # % increase in power3 per case increase
# subsequent rows
while(n < maxN)
{
  n = n + intv                                                              # sample size
  SSiz <- append(SSiz,n)
  oldpower1 = power1
  power1 = 1 - CalBeta(alpha, n, ngrp, f1)                                  # power 1   
  Power1 <- append(Power1, sprintf(power1, fmt="%#.4f")) 
  diff1 = power1 - oldpower1
  Diff1 <- append(Diff1,sprintf(diff1, fmt="%#.4f"))
  incCs1 = diff1 / intv
  IncCs1 <- append(IncCs1,sprintf(incCs1, fmt="%#.4f"))
  PICs1 <- append(PICs1,sprintf(incCs1 / oldpower1 * 100, fmt="%#.1f"))
  oldpower2 = power2 
  power2 = 1 - CalBeta(alpha, n, ngrp, f2)                                  # power 2
  Power2 <- append(Power2, sprintf(power2, fmt="%#.4f")) 
  diff2 = power2 - oldpower2
  Diff2 <- append(Diff2,sprintf(diff2, fmt="%#.4f"))
  incCs2 = diff2 / intv
  IncCs2 <- append(IncCs2,sprintf(incCs2, fmt="%#.4f"))
  PICs2 <- append(PICs2,sprintf(incCs2 / oldpower2 * 100, fmt="%#.1f"))
  oldpower3 = power3
  power3 = 1 - CalBeta(alpha, n, ngrp, f3)                                  # power 3   
  Power3 <- append(Power3, sprintf(power3, fmt="%#.4f")) 
  diff3 = power3 - oldpower3
  Diff3 <- append(Diff3,sprintf(diff3, fmt="%#.4f"))
  incCs3 = diff3 / intv
  IncCs3 <- append(IncCs3,sprintf(incCs3, fmt="%#.4f"))
  PICs3 <- append(PICs3,sprintf(incCs3 / oldpower3 * 100, fmt="%#.1f"))
}
# Combine all vectors to data frame for display
df <- data.frame(SSiz,Power1,Diff1,IncCs1,PICs1,
                      Power2,Diff2,IncCs2,PICs2,
                      Power3,Diff3,IncCs3,PICs3)
df # display results in data frame

The results are as follows

SSiz = sample size per group
Power (1,2,3) power of the 3 models at that sample size
Diff (1,2,3) difference (increase) in power from the previous row
IncCs (1,2,3) increase in power per case increase in sample size
PICs (1,2,3) percent(%) increase (from the previous row) in power per case increase in sample size

df # display results in data frame
   SSiz Power1  Diff1 IncCs1 PICs1 Power2  Diff2 IncCs2 PICs2 Power3  Diff3 IncCs3 PICs3
1     5 0.0898      0      0     0 0.0898      0      0     0 0.1053      0      0     0
2    10 0.1443 0.0545 0.0109  12.1 0.1443 0.0545 0.0109  12.1 0.1809 0.0756 0.0151  14.4
3    15 0.2043 0.0600 0.0120   8.3 0.2043 0.0600 0.0120   8.3 0.2626 0.0817 0.0163   9.0
4    20 0.2663 0.0619 0.0124   6.1 0.2663 0.0619 0.0124   6.1 0.3448 0.0821 0.0164   6.3
5    25 0.3282 0.0619 0.0124   4.6 0.3282 0.0619 0.0124   4.6 0.4242 0.0794 0.0159   4.6
6    30 0.3887 0.0606 0.0121   3.7 0.3887 0.0606 0.0121   3.7 0.4989 0.0747 0.0149   3.5
7    35 0.4470 0.0583 0.0117   3.0 0.4470 0.0583 0.0117   3.0 0.5678 0.0688 0.0138   2.8
8    40 0.5024 0.0554 0.0111   2.5 0.5024 0.0554 0.0111   2.5 0.6301 0.0624 0.0125   2.2
9    45 0.5544 0.0520 0.0104   2.1 0.5544 0.0520 0.0104   2.1 0.6858 0.0557 0.0111   1.8
10   50 0.6028 0.0484 0.0097   1.7 0.6028 0.0484 0.0097   1.7 0.7349 0.0491 0.0098   1.4
11   55 0.6475 0.0446 0.0089   1.5 0.6475 0.0446 0.0089   1.5 0.7776 0.0428 0.0086   1.2
12   60 0.6884 0.0409 0.0082   1.3 0.6884 0.0409 0.0082   1.3 0.8146 0.0370 0.0074   1.0
13   65 0.7255 0.0372 0.0074   1.1 0.7255 0.0372 0.0074   1.1 0.8462 0.0316 0.0063   0.8
14   70 0.7591 0.0336 0.0067   0.9 0.7591 0.0336 0.0067   0.9 0.8731 0.0269 0.0054   0.6
15   75 0.7893 0.0302 0.0060   0.8 0.7893 0.0302 0.0060   0.8 0.8958 0.0227 0.0045   0.5
16   80 0.8163 0.0270 0.0054   0.7 0.8163 0.0270 0.0054   0.7 0.9148 0.0190 0.0038   0.4
17   85 0.8403 0.0240 0.0048   0.6 0.8403 0.0240 0.0048   0.6 0.9306 0.0158 0.0032   0.3
18   90 0.8616 0.0213 0.0043   0.5 0.8616 0.0213 0.0043   0.5 0.9437 0.0131 0.0026   0.3
19   95 0.8803 0.0188 0.0038   0.4 0.8803 0.0188 0.0038   0.4 0.9545 0.0108 0.0022   0.2
20  100 0.8968 0.0165 0.0033   0.4 0.8968 0.0165 0.0033   0.4 0.9634 0.0089 0.0018   0.2

StatsToDo: Sample Size for Comparing Means from multiple Groups (One Way Analysis of Variance)

Example

Pilot study

Sample size

Power

References

Subroutine functions used by all 3 programs

Program 1: Sample size

Program 2: Power

Program 3: Pilot Study