-
Notifications
You must be signed in to change notification settings - Fork 0
/
README.Rmd
executable file
·117 lines (94 loc) · 3.88 KB
/
README.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
---
output:
github_document:
pandoc_args: [
]
---
<!-- README.md is generated from README.Rmd. Please edit that file -->
```{r, echo = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.path = "man/figures/README-"
# fig.path = "docs/images/README-"
# fig.path = "images/README-"
# base.dir = "./"
)
```
# mthapower <img src="./man/figures/mthapower.png" align="right" width="100px" />
[![CRAN version](https://www.r-pkg.org/badges/version/mthapower)](https://cran.r-project.org/package=mthapower){target="blank"}
[![Downloads from Rstudio mirror](https://cranlogs.r-pkg.org/badges/grand-total/mthapower)](https://www.r-pkg.org/pkg/mthapower){target="blank"}
[![DOI](https://zenodo.org/badge/95907604.svg)](https://zenodo.org/badge/latestdoi/95907604){target="blank"}
Calculate sample size and power for association studies involving mitochondrial DNA haplogroups - Based on Samuels et al. AJHG, 2006. 78(4):713-720. [DOI:10.1086/502682](https://www.ncbi.nlm.nih.gov/pmc/PMC1424681){target="blank"}
## Installation
- From CRAN:
```{r cran-installation, eval = FALSE}
install.packages("mthapower")
```
- From GitHub:
```{r gh-installation, eval = FALSE}
# install.packages("devtools")
devtools::install_github("aurora-mareviv/mthapower")
```
- Package DOI: [10.5281/zenodo.3252674](https://doi.org/10.5281/zenodo.3252674){target="blank"}
- Cite as: Baluja A. (2019, June 22). Mthapower - R package version 0.1.1 (Version 0.1.1). Zenodo. DOI: 10.5281/zenodo.3252674
## Shiny app
- Run in Shinyapps.io: [mtDNA_power_calc](https://aurora.shinyapps.io/mtDNA_power_calc/){target="blank"}
- Run locally from Gist:
```{r gh-rungist, eval = FALSE}
# install.packages("shiny")
shiny::runGist('5895082')
```
## Examples
### Sample size estimation
- Determine the minimum number of cases (`Ncmin`), required to detect: either a change from `p0` (haplogroup frequency in controls) to `p1` (haplogroup frequency in cases), or a given OR, with a predefined confidence interval, in a study with `Nh` haplogroups.
```{r example, message=FALSE, warning=FALSE}
library(mthapower)
library(dplyr)
mydata <- mthacases(p0=0.445, Nh=11,
OR.cas.ctrl=c(2), power=80,
sig.level=0.05) # Baudouin study
mydata <- mthacases(p0=0.445, Nh=11,
OR.cas.ctrl=c(1.25,1.5,1.75,2,2.25,2.5,2.75,3),
power=80, sig.level=0.05)
mydata <- mydata[c(2,6)]
mydata %>%
knitr::kable()
plot(mydata)
```
### Power estimation
- For a given study size, determine the minimum effect size that can be detected with the desired power and significance level, in a study with `Nh` haplogroups.
```{r example2a, message=FALSE, warning=FALSE}
# Example 2a:
# library(mthapower)
pow <- mthapower(n.cases=203, p0=0.443, Nh=13, OR.cas.ctrl=2.33, sig.level=0.05)
pow %>%
knitr::kable()
```
```{r example2b, message=FALSE, warning=FALSE}
# Example 2b:
# Create data frames
pow.H150 <- mthapower(n.cases=seq(50,1000,by=50), p0=0.433, Nh=11,
OR.cas.ctrl=1.5, sig.level=0.05)
pow.H175 <- mthapower(n.cases=seq(50,1000,by=50), p0=0.433, Nh=11,
OR.cas.ctrl=1.75, sig.level=0.05)
pow.H200 <- mthapower(n.cases=seq(50,1000,by=50), p0=0.433, Nh=11,
OR.cas.ctrl=2, sig.level=0.05)
pow.H250 <- mthapower(n.cases=seq(50,1000,by=50), p0=0.433, Nh=11,
OR.cas.ctrl=2.5, sig.level=0.05)
# Bind the three data frames:
bindata <- rbind(pow.H150,pow.H175,pow.H200,pow.H250)
# Adds column OR to binded data frame:
bindata$OR <- rep(factor(c(1.50,1.75,2,2.5)),
times = c(nrow(pow.H150),
nrow(pow.H175),
nrow(pow.H200),
nrow(pow.H250)))
# Create plot:
# install.packages("car")
library(car)
scatterplot(power~ncases | OR, regLine=FALSE,
smooth=FALSE,
boxplots=FALSE, by.groups=TRUE,
data=bindata)
```