-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path18-03-2024_cass_BHRC_PRS_covar_bias_testing.R
92 lines (86 loc) · 2.07 KB
/
18-03-2024_cass_BHRC_PRS_covar_bias_testing.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
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
source("functions_to_source.R")
## Does the PRS have association with...
test_bias <- function(prs_column) {
opt <- NULL
df <-
select(prs, IID, {{prs_column}}) %>%
rename(PRS = {{prs_column}})
## State?
## Mann-Whitney U test/Wilcox test
for_test <- inner_join(df, state, by = "IID")
state_test <- wilcox.test(PRS ~ popID, data = for_test)
opt <-
data.frame(
"Statistics" = state_test$statistic,
"P_value" = state_test$p.value
)
## Sex?
## Mann-Whitney U test
for_test <- inner_join(df, sex, by = "IID")
sex_test <- wilcox.test(PRS ~ sex, data = for_test)
opt <-
rbind(opt,
data.frame(
"Statistics" = sex_test$statistic,
"P_value" = sex_test$p.value
))
## Age?
## Kendall test
for_test <- inner_join(df, ages, by = "IID")
agW0_test <- cor.test(for_test$PRS, for_test$age_W0, method = "kendall")
agW1_test <- cor.test(for_test$PRS, for_test$age_W1, method = "kendall")
agW2_test <- cor.test(for_test$PRS, for_test$age_W2, method = "kendall")
opt <-
rbind(opt,
data.frame(
"Statistics" = agW0_test$statistic,
"P_value" = agW0_test$p.value
))
opt <-
rbind(opt,
data.frame(
"Statistics" = agW1_test$statistic,
"P_value" = agW1_test$p.value
))
opt <-
rbind(opt,
data.frame(
"Statistics" = agW2_test$statistic,
"P_value" = agW2_test$p.value
))
rownames(opt) <-
c(
"Wilcox test state",
"Wilcox test sex",
"Kendall test W0",
"Kendall test W1",
"Kendall test W2"
)
return(opt)
}
tests_output <-
list(
AD = test_bias(AD),
ADHD = test_bias(ADHD),
AN = test_bias(AN),
ANX = test_bias(ANX),
ASD = test_bias(ASD),
BD = test_bias(BD),
CD = test_bias(CD),
EA = test_bias(EA),
HEIGHT = test_bias(HEIGHT),
MDD = test_bias(MDD),
OCD = test_bias(OCD),
PTSD = test_bias(PTSD),
SA = test_bias(SA),
SCZ = test_bias(SCZ),
SD = test_bias(SD),
SWB = test_bias(SWB),
TS = test_bias(TS)
)
str(tests_output)
lapply(tests_output, function(data) {
if(any(data$P_value < 0.05)) {
print("Deu merda")
}
})