Power

Power analysis/Sample Size determination

• https://en.m.wikipedia.org/wiki/Power_(statistics)
• Sample size determination from Wikipedia
• Power and Sample Size Determination http://www.stat.wisc.edu/~st571-1/10-power-2.pdf#page=12
• http://biostat.mc.vanderbilt.edu/wiki/pub/Main/AnesShortCourse/HypothesisTestingPart1.pdf#page=40
• Power analysis and sample size calculation for Agriculture (pwr, lmSupport, simr packages are used)
• Why Within-Subject Designs Require Fewer Participants than Between-Subject Designs

Binomial distribution

• Binomial test. Calculating a p-value for a two-tailed test is slightly more complicated, since a binomial distribution isn't symmetric if [math]\displaystyle{ \pi _{0}\neq 0.5 }[/math].
• How To Get The Power Of Test In Hypothesis Testing With Binomial Distribution
• How to Perform a Binomial Test in R
• ?binom.test

Power analysis for default Bayesian t-tests

http://daniellakens.blogspot.com/2016/01/power-analysis-for-default-bayesian-t.html

Using simulation for power analysis

• An example based on a stepped wedge study design
• Framework for power analysis using simulation

Pearson correlation

set.seed(0)

simulate_p_value <- function(sample_size) {
  # Simulate data
  x <- rnorm(sample_size)
  y <- x + rnorm(sample_size)
  
  # Calculate correlation and p-value
  cor.test(x, y)$p.value
}

simulate_p_value_ttest <- function(sample_size) {
  # Simulate data
  x <- rnorm(sample_size, mean = 0)
  y <- rnorm(sample_size, mean = 0.5)
  
  # Calculate p-value
  t.test(x, y)$p.value
}

sample_sizes <- c(10, 30, 50, 100, 300, 500, 1000)
p_values <- sapply(sample_sizes, simulate_p_value)

plot(sample_sizes, p_values, log = "y", type = "o",
     xlab = "Sample Size", ylab = "p-value",
     main = "Effect of Sample Size on p-value")

T-test case

T-test variance could be large if the sample size is small. Below is a case of 30 vs 5 and effect size .5.

set.seed(13)
hist(replicate(100, t.test(rnorm(30), rnorm(5)+.5)$p.value), 
     main='', xlab='pvalue')
set.seed(13); summary(replicate(100, t.test(rnorm(30), rnorm(5)+.5)$p.value))
#      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
# 0.0000157 0.1730494 0.3178627 0.3970372 0.6592276 0.9999164 

Power analysis and sample size calculation for Agriculture

http://r-video-tutorial.blogspot.com/2017/07/power-analysis-and-sample-size.html

Power calculation for proportions (shiny app)

https://juliasilge.shinyapps.io/power-app/

Derive the formula/manual calculation

• One-sample 1-sided test, One sample 2-sided test
• Two-sample 2-sided T test ([math]\displaystyle{ n }[/math] is the sample size in each group)

[math]\displaystyle{ \begin{align} Power & = P_{\mu_1-\mu_2 = \Delta}(\frac{\bar{X}_1 - \bar{X}_2}{\sqrt{\sigma^2/n + \sigma^2/n}} \gt Z_{\alpha /2}) + P_{\mu_1-\mu_2 = \Delta}(\frac{\bar{X}_1 - \bar{X}_2}{\sqrt{\sigma^2/n + \sigma^2/n}} \lt -Z_{\alpha /2}) \\ & \approx P_{\mu_1-\mu_2 = \Delta}(\frac{\bar{X}_1 - \bar{X}_2}{\sqrt{\sigma^2/n + \sigma^2/n}} \gt Z_{\alpha /2}) \\ & = P_{\mu_1-\mu_2 = \Delta}(\frac{\bar{X}_1 - \bar{X}_2 - \Delta}{\sqrt{2 * \sigma^2/n}} \gt Z_{\alpha /2} - \frac{\Delta}{\sqrt{2 * \sigma^2/n}}) \\ & = \Phi(-(Z_{\alpha /2} - \frac{\Delta}{\sqrt{2 * \sigma^2/n}})) \\ & = 1 - \beta =\Phi(Z_\beta) \end{align} }[/math]

Therefore

[math]\displaystyle{ \begin{align} Z_{\beta} &= - Z_{\alpha/2} + \frac{\Delta}{\sqrt{2 * \sigma^2/n}} \\ Z_{\beta} + Z_{\alpha/2} & = \frac{\Delta}{\sqrt{2 * \sigma^2/n}} \\ 2 * (Z_{\beta} + Z_{\alpha/2})^2 * \sigma^2/\Delta^2 & = n \\ n & = 2 * (Z_{\beta} + Z_{\alpha/2})^2 * \sigma^2/\Delta^2 \end{align} }[/math]

# alpha = .05, delta = 200, n = 79.5, sigma=450
1 - pnorm(1.96 - 200*sqrt(79.5)/(sqrt(2)*450)) + pnorm(-1.96 - 200*sqrt(79.5)/(sqrt(2)*450))
# [1] 0.8
pnorm(-1.96 - 200*sqrt(79.5)/(sqrt(2)*450))
# [1] 9.58e-07
1 - pnorm(1.96 - 200*sqrt(79.5)/(sqrt(2)*450)) 
# [1] 0.8

Calculating required sample size in R and SAS

pwr package is used. For two-sided test, the formula for sample size is

[math]\displaystyle{ n_{\mbox{each group}} = \frac{2 * (Z_{\alpha/2} + Z_\beta)^2 * \sigma^2}{\Delta^2} = \frac{2 * (Z_{\alpha/2} + Z_\beta)^2}{d^2} }[/math]

where [math]\displaystyle{ Z_\alpha }[/math] is value of the Normal distribution which cuts off an upper tail probability of [math]\displaystyle{ \alpha }[/math], [math]\displaystyle{ \Delta }[/math] is the difference sought, [math]\displaystyle{ \sigma }[/math] is the presumed standard deviation of the outcome, [math]\displaystyle{ \alpha }[/math] is the type 1 error, [math]\displaystyle{ \beta }[/math] is the type II error and (Cohen's) d is the effect size - difference between the means divided by the pooled standard deviation.

# An example from http://www.stat.columbia.edu/~gelman/stuff_for_blog/c13.pdf#page=3
# Method 1.
require(pwr)
pwr.t.test(d=200/450, power=.8, sig.level=.05,
           type="two.sample", alternative="two.sided")
#
#     Two-sample t test power calculation 
#
#              n = 80.4
#              d = 0.444
#      sig.level = 0.05
#          power = 0.8
#    alternative = two.sided
#
# NOTE: n is number in *each* group

# Method 2.
2*(qnorm(.975) + qnorm(.8))^2*450^2/(200^2)
# [1] 79.5
2*(1.96 + .84)^2*450^2 / (200^2)
# [1] 79.4

And stats::power.t.test() function.

power.t.test(n = 79.5, delta = 200, sd = 450, sig.level = .05,
             type ="two.sample", alternative = "two.sided")
#
#     Two-sample t test power calculation 
#
#              n = 79.5
#          delta = 200
#             sd = 450
#      sig.level = 0.05
#          power = 0.795
#    alternative = two.sided
#
# NOTE: n is number in *each* group

R package related to power analysis

CRAN Task View: Design of Experiments

• powerAnalysis w/o vignette
• powerbydesign w/o vignette
• easypower w/ vignette
• pwr w/ vignette, https://www.statmethods.net/stats/power.html. The reference is Cohen's book.
  • Power analysis in Statistics with R
• powerlmm Power Analysis for Longitudinal Multilevel/Linear Mixed-Effects Models.
• ssize.fdr w/o vignette
• samplesize w/o vignette
• ssizeRNA w/ vignette
• power.t.test(), power.anova.test(), power.prop.test() from stats package

RNA-seq

• Feasibility of sample size calculation for RNA-seq studies Poplawski 2018.
  • The ‘Scotty’ (MATLAB) tool performed best
  • ‘Scotty’, ssizeRNA 2016 and PROPER 2015 generated comparable results.
  • Bi et al. showed that ssizeRNA provided a more accurate estimate of power/sample size than RnaSeqSampleSize; ssizeRNA and RnaSeqSampleSize provided results much faster than PROPER. See Power and sample size calculations for high-throughput sequencing-based experiments 2018
• Empirical assessment of the impact of sample number and read depth on RNA-Seq analysis workflow performance
• Number of replicates needed is dependent on the following that varies from gene to gene
  • Within group variance
  • Read coverage
  • Desired detectable effect size
• powsimR: power analysis for bulk and single cell RNA-seq experiments Vieth 2017
• RnaSeqSampleSize
• RNASeqPower. It was used in GREIN: An Interactive Web Platform for Re-analyzing GEO RNA-seq Data 2019

ScRNA-seq

scRNA-seq

Russ Lenth Java applets

https://homepage.divms.uiowa.edu/~rlenth/Power/index.html

Bootstrap method

The upstrap Crainiceanu & Crainiceanu, Biostatistics 2018

Multiple Testing Case

Optimal Sample Size for Multiple Testing The Case of Gene Expression Microarrays

Unbalanced randomization

Can unbalanced randomization improve power?

Yes, unbalanced randomization can improve power, in some situations

Prediction

Minimum sample size for external validation of a clinical prediction model with a binary outcome Riley 2021

High dimensional data

• Simulation To Understand Needed Sample Sizes from the online book Biostatistics for Biomedical Research
• MKmisc package
  • Sample size planning for developing classifiers using high-dimensional DNA microarray data Biostatistics 2007 & online calculator.
    • Goal: Determine the sample size required to ensure that the expected correct classification probability for a predictor developed from training data is within γ of the optimal expected correct classification probability for a linear classification problem (p6).
    • Input: standard fold-change, number of genes, prevalence in the larger group, tolerance (cannot control). The standard fold-change is calculated by [math]\displaystyle{ 0.8 * \max|{t_g}| * \sqrt{1/n1 + 1/n2} }[/math] where [math]\displaystyle{ t_g }[/math] is the T-statistic for gene [math]\displaystyle{ g }[/math].
  • ssize.pcc(). Input: standard fold-change, number of genes, prevalence in the larger group, tolerance, number of significant genes.
  • Vignette ssize.pcc, which is based on the probability of correct classification (PCC)

Survival data

Sample size and predictive performance of machine learning methods with survival data: A simulation study 2023. R code is in the supplementary.