EpipwR
Introduction
EpipwR is a quasi-simulation based approach to power analysis for
epigenome-wide association studies (EWAS) with a continuous or binary
outcome variable. It is fast, efficient and allows the exact level of
precision to be set by the user. EpipwR is the first open source tool
able to calculate power for EWAS with a continuous outcome. For binary
outcomes, it re-implements the pwrEWAS methodology (Graw et al., 2019 -
currently available on GitHub) with a few changes and significantly
improves the computation time. EpipwR is simple in that it has three
main functions: get_power_cont() and
get_power_cc calculate power for various settings with
continuous and binary outcomes (respectively) and
EpipwR_plot() plots the output. The following vignette
describes in detail the EpipwR workflow for an analysis of a continuous
outcome.
Methodology (Continuous Outcome)
Suppose an upcoming EWAS plans to analyze the methylation levels at K CpG sites and some continuous outcome Y. Of the K CpG sites, it is assumed that only a subest of these (Km ) have non-zero correlation to Y. The goal of EpipwR is to determine the level of power achieved by this study at a given sample size n, where power is defined as the proportion of the Km CpG sites are found to have a statistically significant association with Y.
To begin, EpipwR generates data for the Km CpG methylation levels and for Y. Users specify μρ and σρ, the average and standard deviation of the correlation between each of the Km methylation levels and Y. If users prefer that all Km methylation levels have equal correlation with Y, then σρ can be set to 0. Like many power analysis tools, EpipwR assumes no correlation across CpG sites (i.e., they are conditionally independent given Y). Users also select which empirical EWAS data to use to generate initial methylation levels. Full details of the data generation process can be found in the corresponding EpipwR article (Barth and Reynolds, 2024).
Once a dataset has been generated, p-values are calculated by performing the specified hypothesis test (Pearson, Spearman or Kendall correlation) for each Km site. Please note that in absence of any “vertical” correlation or additional covariates, performing a standard analysis in limma is equivalent to using the Pearson correlation test. Users can specify either a false discovery rate or family-wise type I error rate correction for determining statistical significance. If the false discovery rate is selected, p-values for the remaining null tests are generated via a Uniform(0, 1) distribution.
Since the process described above is based heavily on random sampling, it must be repeated on several data sets to produce a reliable power estimate. Whereas traditional tools require the user to specify the number of repetitions, EpipwR sets 20 to the minimum number of data sets, allows users to set a maximum number of data sets and a level of precision in the form of a 95% confidence interval margin of error. Anytime after the first 20 datasets, EpipwR terminates once one of these two conditions are met. The default margin of error is 3% (it is recommended that this be between 1 and 5%). For example, if the target margin of error is 3% and the maximum number of data sets is 1000, then the algorithm stops when N = 1000 or once the following conditions are met:
$$1.96\cdot\frac{s_{(N)}}{\sqrt{N}}<.03,\qquad N\ge20 $$
where s(N) is the standard deviation of average power over N data sets.
Function Inputs
The below table describes the details of the
get_power_cont inputs.
| Parameter | Description |
dm |
The total number of non-null tests (i.e., the number of CpG sites with methylation having a non-zero correlation with the phenotype). |
Total |
The total number of tests (i.e., total number of CpG sites in the study). |
n |
The sample size of the study - accepts a vector of values. |
fdr_fwer |
Either the false discovery rate (FDR) or the family-wise type I
error rate (FWER), depending on use_fdr. |
rho_mu |
The average correlation of the non-null tests - accepts a vector of values. |
rho_sd=0 |
The standard deviation of the correlation for non-null tests. |
Tissue="Saliva" |
Tissue corresponding to empirical EWAS data (see
?get_power_cont for valid options). |
Nmax=1000 |
The maximum number of data sets that EpipwR will use to determine power. Cannot be less than 20. |
MOE=.03 |
The target margin of error of a 95% confidence interval for average power. |
test="Pearson" |
Type of correlation/test utilized. Acceptable inputs are
"Pearson", "Kendall" and
"Spearman" . |
use_fdr=T |
Indicates that the false discovery rate should be used as the correction. Otherwise, uses FWER. |
Suppress_Updates=F |
If T , removes intermediate status updates. |
Methodology (Binary outcome)
For the case-control setting, EpipwR re-implements the pwrEWAS
methodology (Graw et al., 2019) with a few key changes to improve speed
and interpretation, namely (1) keeping the precision equal between
groups rather than variance, (2) requiring users to specify mean and
variance for an effect size distribution (rather than a target maximal
difference) and (3) generating p-values directly for null tests. Like in
get_power_cont, get_power_cc calculates the
optimal amount of data sets needed to determine power at a
user-specified level of precision.
EpipwR Workflow
Suppose the user is planning an EWAS of 100,000 CpG sites with an estimated 100 of those having a non-zero correlation with continuous phenotype. They would like to test sample sizes of 70, 80, 90, 100 and fixed correlations of 0.3, 0.4 and 0.5. They plan to use an FDR of 5% on Pearson correlation tests, the Saliva empirical EWAS and 1,000 as the maximum number of data sets. Since they hope to have very accurate results, the margin of error is set to be 1%. They would then execute the following code:
output <- get_power_cont(
dm=100
,Total=100000
,n=c(70,80,90,100)
,fdr_fwer=.05
,rho_mu = c(0.3, 0.4, 0.5)
,rho_sd = 0
,Tissue = "Saliva"
,Nmax=1000
,MOE=.01
,test="pearson"
,use_fdr=T
,Suppress_updates=F
)
outputwhich produces the following data frame.
Note that executing this code block will produce very similar but not
identical results, since EpipwR is based on random sampling. The first
two columns track the setting for which power was calculated (each row
is a unique combination of n and rho_mu ). The
next four columns show the average power for this setting, the standard
deviation of power, the number of data sets used for this setting, and
finally the standard error of the average power (column 4 divided by the
square root of column 5). Note that when power is very close to 0, often
the minimum number of data sets is sufficient to achieve the desired
level of precision - the same is true when power is close to 1.
Finally, users can also summarize the output with error bar plots
with EpipwR_plot(output). For our previous example, this
produces the following plot:
Please note that EpipwR_plot can only run with the data
frame output of get_power_cont or get_power_cc
- it will not work with data in any other format.
# Load sessioninfo package
library(sessioninfo)
# Display session information
session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#> setting value
#> version R version 4.4.1 (2024-06-14)
#> os Ubuntu 24.04.1 LTS
#> system x86_64, linux-gnu
#> ui X11
#> language (EN)
#> collate C
#> ctype en_US.UTF-8
#> tz Etc/UTC
#> date 2024-10-30
#> pandoc 3.2.1 @ /usr/local/bin/ (via rmarkdown)
#>
#> ─ Packages ───────────────────────────────────────────────────────────────────
#> package * version date (UTC) lib source
#> AnnotationDbi 1.69.0 2024-10-30 [2] https://bioc.r-universe.dev (R 4.4.1)
#> AnnotationHub 3.15.0 2024-10-30 [2] https://bioc.r-universe.dev (R 4.4.1)
#> Biobase 2.67.0 2024-10-30 [2] https://bioc.r-universe.dev (R 4.4.1)
#> BiocFileCache 2.15.0 2024-10-30 [2] https://bioc.r-universe.dev (R 4.4.1)
#> BiocGenerics 0.53.0 2024-10-30 [2] https://bioc.r-universe.dev (R 4.4.1)
#> BiocManager 1.30.25 2024-08-28 [2] RSPM (R 4.4.0)
#> BiocVersion 3.21.1 2024-10-30 [2] https://bioc.r-universe.dev (R 4.4.1)
#> Biostrings 2.75.0 2024-10-30 [2] https://bioc.r-universe.dev (R 4.4.1)
#> bit 4.5.0 2024-09-20 [2] RSPM (R 4.4.0)
#> bit64 4.5.2 2024-09-22 [2] RSPM (R 4.4.0)
#> blob 1.2.4 2023-03-17 [2] RSPM (R 4.4.0)
#> bslib 0.8.0 2024-07-29 [2] RSPM (R 4.4.0)
#> buildtools 1.0.0 2024-10-28 [3] local (/pkg)
#> cachem 1.1.0 2024-05-16 [2] RSPM (R 4.4.0)
#> cli 3.6.3 2024-06-21 [2] RSPM (R 4.4.0)
#> crayon 1.5.3 2024-06-20 [2] RSPM (R 4.4.0)
#> curl 5.2.3 2024-09-20 [2] RSPM (R 4.4.0)
#> DBI 1.2.3 2024-06-02 [2] RSPM (R 4.4.0)
#> dbplyr 2.5.0 2024-03-19 [2] RSPM (R 4.4.0)
#> digest 0.6.37 2024-08-19 [2] RSPM (R 4.4.0)
#> dplyr 1.1.4 2023-11-17 [2] RSPM (R 4.4.0)
#> EpipwR * 1.1.0 2024-10-30 [1] https://bioc.r-universe.dev (R 4.4.1)
#> EpipwR.data 0.99.0 2024-10-22 [2] Bioconductor 3.20 (R 4.4.1)
#> evaluate 1.0.1 2024-10-10 [2] RSPM (R 4.4.0)
#> ExperimentHub 2.13.1 2024-10-28 [2] https://bioc.r-universe.dev (R 4.4.1)
#> fansi 1.0.6 2023-12-08 [2] RSPM (R 4.4.0)
#> fastmap 1.2.0 2024-05-15 [2] RSPM (R 4.4.0)
#> filelock 1.0.3 2023-12-11 [2] RSPM (R 4.4.0)
#> generics 0.1.3 2022-07-05 [2] RSPM (R 4.4.0)
#> GenomeInfoDb 1.41.2 2024-10-02 [2] https://bioc.r-universe.dev (R 4.4.1)
#> GenomeInfoDbData 1.2.13 2024-10-30 [2] Bioconductor
#> glue 1.8.0 2024-09-30 [2] RSPM (R 4.4.0)
#> htmltools 0.5.8.1 2024-04-04 [2] RSPM (R 4.4.0)
#> httr 1.4.7 2023-08-15 [2] RSPM (R 4.4.0)
#> IRanges 2.39.2 2024-10-25 [2] https://bioc.r-universe.dev (R 4.4.1)
#> jquerylib 0.1.4 2021-04-26 [2] RSPM (R 4.4.0)
#> jsonlite 1.8.9 2024-09-20 [2] RSPM (R 4.4.0)
#> KEGGREST 1.45.1 2024-10-16 [2] https://bioc.r-universe.dev (R 4.4.1)
#> knitr 1.48 2024-07-07 [2] RSPM (R 4.4.0)
#> lifecycle 1.0.4 2023-11-07 [2] RSPM (R 4.4.0)
#> magrittr 2.0.3 2022-03-30 [2] RSPM (R 4.4.0)
#> maketools 1.3.1 2024-10-28 [3] Github (jeroen/maketools@d46f92c)
#> memoise 2.0.1 2021-11-26 [2] RSPM (R 4.4.0)
#> pillar 1.9.0 2023-03-22 [2] RSPM (R 4.4.0)
#> pkgconfig 2.0.3 2019-09-22 [2] RSPM (R 4.4.0)
#> png 0.1-8 2022-11-29 [2] RSPM (R 4.4.0)
#> R6 2.5.1 2021-08-19 [2] RSPM (R 4.4.0)
#> rappdirs 0.3.3 2021-01-31 [2] RSPM (R 4.4.0)
#> rlang 1.1.4 2024-06-04 [2] RSPM (R 4.4.0)
#> rmarkdown * 2.28 2024-08-17 [2] RSPM (R 4.4.0)
#> RSQLite 2.3.7 2024-05-27 [2] RSPM (R 4.4.0)
#> S4Vectors 0.43.2 2024-10-17 [2] https://bioc.r-universe.dev (R 4.4.1)
#> sass 0.4.9 2024-03-15 [2] RSPM (R 4.4.0)
#> sessioninfo * 1.2.2 2021-12-06 [2] RSPM (R 4.4.0)
#> sys 3.4.3 2024-10-04 [2] RSPM (R 4.4.0)
#> tibble 3.2.1 2023-03-20 [2] RSPM (R 4.4.0)
#> tidyselect 1.2.1 2024-03-11 [2] RSPM (R 4.4.0)
#> UCSC.utils 1.1.0 2024-10-29 [2] https://bioc.r-universe.dev (R 4.4.1)
#> utf8 1.2.4 2023-10-22 [2] RSPM (R 4.4.0)
#> vctrs 0.6.5 2023-12-01 [2] RSPM (R 4.4.0)
#> xfun 0.48 2024-10-03 [2] RSPM (R 4.4.0)
#> XVector 0.45.0 2024-10-02 [2] https://bioc.r-universe.dev (R 4.4.1)
#> yaml 2.3.10 2024-07-26 [2] RSPM (R 4.4.0)
#> zlibbioc 1.51.2 2024-10-21 [2] Bioconductor 3.20 (R 4.4.1)
#>
#> [1] /tmp/RtmpdAcBbP/Rinst14dc142bef8b
#> [2] /github/workspace/pkglib
#> [3] /usr/local/lib/R/site-library
#> [4] /usr/lib/R/site-library
#> [5] /usr/lib/R/library
#>
#> ──────────────────────────────────────────────────────────────────────────────References
Barth, J., Reynolds, A.W. (2024). EpipwR: Efficient Power Analysis for EWAS with Continuous Outcomes [preprint]. bioRxiv, 2024.09.06.611713; doi: https://doi.org/10.1101/2024.09.06.611713
Graw, S., Henn, R., Thompson, J. A., and Koestler, D. C. (2019). pwrEWAS: A user-friendly tool for comprehensive power estimation for epigenome wide association studies (EWAS). BMC Bioinformatics, 20(1):218.