-
Notifications
You must be signed in to change notification settings - Fork 0
/
diff_ini.R
71 lines (53 loc) · 1.78 KB
/
diff_ini.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
##### Simulations to investigate the sensitivity of LSMM to initial parameter specification #####
# Vary alpha=0.2, 0.4, 0.6,
# pi1=0.01, 0.05, 0.1, 0.15, 0.2 to get
# Figures S24 in Supplementary Document
library(LSMM)
library(MASS)
# function to generate data
generate_data_TGM <- function(M, alpha, pi1){
index <- sample(M, M*pi1)
gamma <- numeric(M)
gamma[index] <- 1
Pvalue <- runif(M)
Pvalue[index] <- rbeta(length(index), alpha, 1)
return(list(Pvalue = Pvalue, gamma = gamma))
}
# evaluate FDR
perf_FDR <- function(true, est, fdr){
fdr <- as.numeric(fdr)
t <- table(true, est)
if (sum(est)==0){
FDR.fit <- 0
}
else if (sum(est)==length(est)){
FDR.fit <- t[1]/(t[1]+t[2])
}
else{
FDR.fit <- t[1,2]/(t[1,2]+t[2,2])
}
return(FDR.fit)
}
M <- 100000 # No. of SNPs
alpha <- 0.2 # parameter in the Beta distribution
pi1 <- 0.1 # proportion of risk SNPs
rep <- 50 # repeat times
result <- matrix(0, rep, 4)
for (i in 1:rep){
cat(i, "out of", rep, "\n")
data <- generate_data_TGM(M, alpha, pi1)
# default
fit1 <- LSMM(data$Pvalue)
assoc.SNP1 <- assoc.SNP(fit1, FDRset = 0.1, fdrControl = "global")
result[i, 1] <- fit1$pi1_
result[i, 2] <- perf_FDR(data$gamma, assoc.SNP1$gamma, 1-fit1$pi1)
# random initial values
alpha_ini <- runif(1, 0.1, 0.6)
pi1_ini <- runif(1, 0, 0.3)
fit2 <- LSMM(data$Pvalue, alpha = alpha_ini, pi1_ = pi1_ini)
assoc.SNP2 <- assoc.SNP(fit2, FDRset = 0.1, fdrControl = "global")
result[i, 3] <- fit2$pi1_
result[i, 4] <- perf_FDR(data$gamma, assoc.SNP2$gamma, 1-fit2$pi1)
}
result <- as.data.frame(result)
names(result) <- c("est.default", "FDR.default", "est.random", "FDR.random")