-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path22-04-2024_cass_BHRC_PRS_association_tests.R
73 lines (62 loc) · 1.88 KB
/
22-04-2024_cass_BHRC_PRS_association_tests.R
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
source("functions_to_source.R")
vADHD <-
prs_v2 %>%
select(IID, ADHD) %>%
rename(PRS = ADHD)
aADHD <-
ages %>%
tidyr::pivot_longer(
cols = c(2, 3, 4),
names_to = "wave",
values_to = "age")
aADHD$wave <- gsub("age_", "", aADHD$wave)
pADHD <-
pheno %>%
select(IID, wave, dcanyhk) %>%
rename(diagnosis = 3)
data <-
plyr::join_all(list(sex, vADHD, state, aADHD), by = "IID", type = "inner") %>%
inner_join(., pADHD, by = c("IID", "wave")) %>%
filter(IID %in% ages$IID)
non_redudant_data <-
data %>%
select(IID, sex, popID, PRS) %>%
unique()
### Teste de normalidade
y <- rnorm(1363)
ks.test(unique(data$PRS), y)
### Teste de associação
## sexo
wilcox.test(PRS ~ sex, non_redudant_data)
library(rstatix)
library(coin)
wilcox_effsize(non_redudant_data, PRS ~ sex)
## estado
wilcox.test(PRS ~ popID, non_redudant_data)
wilcox_effsize(non_redudant_data, PRS ~ popID)
## diagnóstico (falta a força da associação)
# All waves (do i put :age to see only the PRS and diagnosis association?)
wilcox.test(PRS ~ diagnosis, data)
wilcox_effsize(data, PRS ~ diagnosis)
# W0
wilcox.test(PRS ~ diagnosis, filter(data, wave == "W0"))
wilcox_effsize(filter(data, wave == "W0"), PRS ~ diagnosis)
# W1
wilcox.test(PRS ~ diagnosis, filter(data, wave == "W1"))
wilcox_effsize(filter(data, wave == "W1"), PRS ~ diagnosis)
# W2
wilcox.test(PRS ~ diagnosis, filter(data, wave == "W2"))
wilcox_effsize(filter(data, wave == "W2"), PRS ~ diagnosis)
## Parent diagnosis association
data_pt <-
filter(parents, IID %in% data$IID) %>%
inner_join(., vADHD, by = "IID")
wilcox.test(PRS ~ ADHD, data_pt)
wilcox_effsize(data_pt, PRS ~ ADHD)
### Age association
# W0
summary(lm(PRS ~ age, filter(data, wave == "W0")))$coefficients
# W1
summary(lm(PRS ~ age, filter(data, wave == "W1")))$coefficients
# W2
summary(lm(PRS ~ age, filter(data, wave == "W2")))$coefficients