-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathExercise_PowerCalc_5.3.2021.Rmd
More file actions
40 lines (37 loc) · 1.08 KB
/
Exercise_PowerCalc_5.3.2021.Rmd
File metadata and controls
40 lines (37 loc) · 1.08 KB
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
# 5.3.2020
# This script calculates power for observing effects of exercise!
# thuyvy's power function:
```{r}
gotta_power = function(data, trait, power_per) {
#'@param trait - ECG trait
#'@param power_per - power, as a decimal
ptab <- cbind(NULL, NULL)
for (i in seq(from = 0.02, to = 1, by = 0.02)) {
pwrt <-
power.t.test(
delta = ((mean(data[,trait])*i)+mean(data[,trait]))-mean(data[,trait]),
sd = sd(data[,trait]),
power = power_per,
type = "two.sample",
alternative = "two.sided"
)
ptab <- rbind(ptab, cbind(pwrt$delta, pwrt$n, i*100))
}
power_table = as.data.frame(ptab)
names(power_table) = c("delta", "n", "percent")
power_table<<-power_table
}
```
# read data
```{r}
df = readRDS('/Volumes/JHPCE 2/dcs01/active/projects/SHAPE/R_objects/final_shape35_dataset.rds')
noNA = df[-which(is.na(df$mtDNA_change)),]
gotta_power(noNA, 'mtDNA_change', 0.8)
power.t.test(
n = nrow(noNA),
sd = sd(noNA$mtDNA_change),
power = 0.8,
type = "two.sample",
alternative = "two.sided"
)
```