Estimates the average power and the number of replicates required by a future experiment, starting from the effect sizes and the variances observed in a differential analysis. The multiplicity is handled by controlling the False Discovery Rate following the approach of Jung (2005) and its closed-form simplification by Liu & Hwang (2007). The calculation is self-contained: it only uses the fold changes and the per-group standard deviations stored in the results table, and requires no external power/sample-size package.
estimate.power(
DEprot.analyses.object,
contrast = 1,
sample.size.range = c(2, 30),
target.power = 0.8,
fdr = NULL,
pi0 = NULL,
pi0.lambda = 0.5,
effect.size = "empirical",
desired.FC = NULL,
sd.quantile = 0.5,
min.effect.size = NULL,
hedges.correction = TRUE,
approximation = "t",
max.iterations = 100,
tolerance = 1e-06,
line.color = "black",
power.threshold.line.color = "firebrick"
)An object of class DEprot.analyses.
Number indicating the position of the contrast to use. Default: 1.
Numeric vector of length 2 indicating the minimum and maximum number of samples per group to test. Default: c(2,30).
Number between 0 and 1 indicating the average power to reach. Default: 0.8.
Number between 0 and 1 indicating the False Discovery Rate to control. Default: NULL, the padj.th stored in the object is used.
Number between 0 and 1 indicating the proportion of non-differential proteins. Default: NULL, estimated from the p-value distribution.
Number between 0 and 1 indicating the tuning parameter used to estimate pi0 (Storey's estimator: mean(p > lambda)/(1 - lambda)). Default: 0.5.
Either the string 'empirical', to average the power over the effect sizes observed for the responsive proteins, or a numeric value indicating a single standardized effect size (Cohen's d) to plan for. Default: "empirical".
Numeric value indicating a desired fold change (linear scale) to plan for. When provided, the effect size is computed as log2(desired.FC)/sd, where sd is the quantile of the observed pooled standard deviations defined by sd.quantile. It overrides effect.size. Default: NULL.
Number between 0 and 1 indicating the quantile of the observed pooled standard deviations to use together with desired.FC. Default: 0.5 (median).
Numeric value indicating the minimum absolute effect size to keep among the responsive proteins. Default: NULL, no filtering.
Logical value indicating whether the small-sample bias of the standardized effect sizes should be corrected (Hedges' g). Default: TRUE.
String indicating the distribution used for the power calculation. One among: 't' (non-central t, exact) and 'normal' (normal approximation, as in the original Jung's formulation). Default: "t".
Numeric value indicating the maximum number of iterations allowed to solve the average power at each sample size. Default: 100.
Numeric value indicating the tolerance used to solve the average power at each sample size. Default: 1e-6.
String indicating any R-supported color to use for the curves. Default: "black".
String indicating any R-supported color to use for the target power line. Default: "firebrick".
An object of class DEprot.power.
For a given number of replicates n per group, the average power over the m1 responsive
proteins and the per-test significance level alpha are linked by the FDR:
alpha = FDR * m1 * (1-beta) / ((1 - FDR) * m0), with m1 = m * (1 - pi0) and m0 = m - m1.
Since the two quantities depend on each other, at each sample size the couple is obtained by root finding
on power(alpha(x)) - x, taking the largest root because x = 0 is always a trivial solution.
The realized FDR reported in the results is therefore equal to the level requested by construction.
The power of the single protein is obtained from a two-sided two-sample t-test with non-centrality parameter
d * sqrt(n/2) (d * sqrt(n) and df = n-1 for paired designs).
The standardized effect sizes are computed as log2(FoldChange) divided by the pooled standard deviation
reconstructed from the sd.<group> columns, and the per-protein group sizes are recovered from the
sd/sem pair, so that proteins with missing values are treated correctly. For paired analyses
the effect size is a Cohen's dz, obtained from lfcSE; when lfcSE is not available
(stat.test = "wilcoxon") the pooled standard deviation is used instead, which ignores the
within-pair correlation and makes the estimate conservative.
The function works on the output of all the differential functions of DEprot (diff.analyses,
diff.analyses.limma, diff.analyses.prolfqua and diff.analyses.proDA), since it only needs
the fold changes and the per-group dispersions. For diff.analyses.proDA the standard deviations are
computed on the observed values only, so for the proteins affected by dropout they are estimated on the upper
part of the intensity distribution and are under-estimated: restricting the estimation to the proteins
quantified in all the samples (min.effect.size is not enough for this) gives a more conservative answer.
Two further limitations are worth keeping in mind. Effect sizes estimated from a pilot experiment are inflated for the
proteins that were selected because they came out significant (winner's curse), so the empirical mode is
optimistic when pi0 is high; planning on desired.FC avoids the problem. Moreover, when the
differential analyses have been run on imputed counts the variance is compressed and every estimate below
becomes optimistic: a warning is raised in that case.
Jung SH (2005), Sample size for FDR-control in microarray data analysis. Bioinformatics 21(14):3097-3104.
Liu P, Hwang JTG (2007), Quick calculation for sample size while controlling false discovery rate with application to microarray analysis. Bioinformatics 23(6):739-746.
pwr <- estimate.power(DEprot::test.toolbox$diff.exp.limma, contrast = 1)
#> Warning: The differential analyses have been performed on imputed counts ('imputed'): the imputation compresses the variance and all the estimates below are optimistic.
#> Consider re-running the differential analyses on the normalized counts, restricted to the proteins quantified in all the samples of the contrast.
pwr
#> Contrast | FBS.vs.6h.DMSO (paired)
#> Counts used | imputed (limma)
#> Proteins tested | 50 (m1 = 32, pi0 = 0.36)
#> Effect size | empirical, median |d| = 0.87
#> FDR | 0.05
#> Current n/group | 4
#> Required n/group | 13 (average power 0.8)
#>
pwr@power.table
#> n.per.group alpha average.power expected.TP expected.FP
#> 1 2 9.356725e-10 0.00000001 0.00000032 1.684211e-08
#> 2 3 9.356725e-10 0.00000001 0.00000032 1.684211e-08
#> 3 4 4.076387e-03 0.04356638 1.39412425 7.337496e-02
#> 4 5 2.805862e-02 0.29987653 9.59604892 5.050552e-01
#> 5 6 4.194378e-02 0.44827415 14.34477286 7.549880e-01
#> 6 7 5.138867e-02 0.54921646 17.57492673 9.249961e-01
#> 7 8 5.839984e-02 0.62414825 19.97274415 1.051197e+00
#> 8 9 6.384600e-02 0.68235416 21.83533322 1.149228e+00
#> 9 10 6.819880e-02 0.72887469 23.32398995 1.227578e+00
#> 10 11 7.174826e-02 0.76680949 24.53790375 1.291469e+00
#> 11 12 7.468648e-02 0.79821179 25.54277734 1.344357e+00
#> 12 13 7.714737e-02 0.82451250 26.38439997 1.388653e+00
#> 13 14 7.922790e-02 0.84674814 27.09594045 1.426102e+00
#> 14 15 8.100045e-02 0.86569232 27.70215412 1.458008e+00
#> 15 16 8.252043e-02 0.88193706 28.22198584 1.485368e+00
#> 16 17 8.383101e-02 0.89594388 28.67020410 1.508958e+00
#> 17 18 8.496642e-02 0.90807862 29.05851599 1.529396e+00
#> 18 19 8.595417e-02 0.91863514 29.39632444 1.547175e+00
#> 19 20 8.681657e-02 0.92785213 29.69126830 1.562698e+00
#> 20 21 8.757197e-02 0.93592543 29.94961387 1.576295e+00
#> 21 22 8.823552e-02 0.94301707 30.17654616 1.588239e+00
#> 22 23 8.881987e-02 0.94926237 30.37639588 1.598758e+00
#> 23 24 8.933558e-02 0.95477406 30.55276979 1.608041e+00
#> 24 25 8.979168e-02 0.95964863 30.70875609 1.616250e+00
#> 25 26 9.019577e-02 0.96396724 30.84695170 1.623524e+00
#> 26 27 9.055433e-02 0.96779940 30.96958066 1.629978e+00
#> 27 28 9.087296e-02 0.97120472 31.07855092 1.635713e+00
#> 28 29 9.115644e-02 0.97423444 31.17550197 1.640816e+00
#> 29 30 9.140898e-02 0.97693343 31.26186986 1.645362e+00
#> expected.discoveries expected.FDR
#> 1 3.368421e-07 0.05
#> 2 3.368421e-07 0.05
#> 3 1.467499e+00 0.05
#> 4 1.010110e+01 0.05
#> 5 1.509976e+01 0.05
#> 6 1.849992e+01 0.05
#> 7 2.102394e+01 0.05
#> 8 2.298456e+01 0.05
#> 9 2.455157e+01 0.05
#> 10 2.582937e+01 0.05
#> 11 2.688713e+01 0.05
#> 12 2.777305e+01 0.05
#> 13 2.852204e+01 0.05
#> 14 2.916016e+01 0.05
#> 15 2.970735e+01 0.05
#> 16 3.017916e+01 0.05
#> 17 3.058791e+01 0.05
#> 18 3.094350e+01 0.05
#> 19 3.125397e+01 0.05
#> 20 3.152591e+01 0.05
#> 21 3.176479e+01 0.05
#> 22 3.197515e+01 0.05
#> 23 3.216081e+01 0.05
#> 24 3.232501e+01 0.05
#> 25 3.247048e+01 0.05
#> 26 3.259956e+01 0.05
#> 27 3.271426e+01 0.05
#> 28 3.281632e+01 0.05
#> 29 3.290723e+01 0.05