-
Notifications
You must be signed in to change notification settings - Fork 1
/
cpval.R
96 lines (63 loc) · 2.63 KB
/
cpval.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
93
94
95
96
set.seed(1234)
## 400 true null generated by standard normal, 100 false null generated by mean 3.5 norm
Zstat <- c(rnorm(400), rnorm(n=100, mean = 10))
pval <- 1-pnorm(Zstat)
## Fisher's combining function
# fixed selection cutoff t
t <- qchisq(1-0.05/500,200)
b <- c()
for (i in seq_along(pval)){
b[i] <- min(c(exp(-t/2)/prod(pval[-i]),1))
# b[i] <- exp(-t/2)/prod(pval[-i])
}
cpval <- pval/b
## Stouffer's
## minP combining function
t = 0.5
b <- c()
for (i in seq_along(pval)){
# b[i] <- min(c(exp(-t/2)/prod(pval[-i]),1))
b[i] <- ifelse(min(pval[-i])>t, t, 1)
}
cpval <- pval/b
hist(pval)
## Fisher Combining Function
tf <- seq(0,10)
p2 <- seq(0,1,0.01)
f <- function(p2,tf){ pmin(exp(-tf/2)/p2,1)}
b1 <- outer(p2,tf,f)
plot(p2, b1[,1],ylim = c(0,1), xlab= expression(p[2]), ylab = expression(b[1]), type = "l", main="inflation factor versus the other p-value")
for (i in 2:length(t)){
lines(p2,b1[,i])
}
text(c(0.95,0.77,0.595,0.44, 0.32,0.22,0.14, 0.04),c(0.97,0.75,0.58,0.47,0.38,0.315, 0.2, 0.08),
labels = c("t=0", "t=1","t=2", "t=3", "t=4", "t=5", ". . . . . .", "t=10"))
plot(t, b1[1,],ylim = c(0,1), xlab= "t", ylab = expression(b[1]), type = "l", main="inflation factor versus selection cutoff")
for (i in c(0.2,0.4, 0.6, 1:10)){
lines(t,b1[i*10+1,])
}
text(c(9.7,7.6,6.5,6.2,5.7,5,4.7, 4,3.5),c(0.97,0.88,0.74,0.6,0.475,0.335, 0.245, 0.19,0.13),
labels = c("p2=0", "p2=0.02","p2=0.04", "p2=0.06", "p2=0.1", "p2=0.2", "p2=0.3", ". . . . . .", "p2=1"))
persp(p2, t, b1, theta = 60, phi = 30, expand = 0.5, col = "lightblue",
ltheta = 120, shade = 0.75, ticktype = "detailed",
xlab = "p2", ylab = "t", zlab = "b1",main = "3D Surface plot for inflation factor b1"
)
## MinP Combining Function
t <- seq(0,1,0.1)
p2 <- seq(0,1,length=100)
f <- function(p2,t){ ifelse(p2>t, t, 1) }
b1 <- outer(p2,t,f)
plot(p2, b1[,1],ylim = c(0,1), xlab= expression(p[2]), ylab = expression(b[1]), type = "l", main="inflation factor versus the other p-value")
for (i in 2:length(t)){
lines(p2,b1[,i])
}
#text(c(0.95,0.77,0.595,0.44, 0.32,0.22,0.14, 0.04),c(0.97,0.75,0.58,0.47,0.38,0.315, 0.2, 0.08),
# labels = c("t=0", "t=1","t=2", "t=3", "t=4", "t=5", ". . . . . .", "t=10"))
persp(p2, t, b1, theta = -45, phi = 30, expand = 0.5, col = "lightblue",
ltheta = 120, shade = 0.75, ticktype = "detailed",
xlab = "p2", ylab = "t", zlab = "b1",main = "3D Surface plot for inflation factor b1"
)
persp(p2, t, b1, phi = 45, theta = 45,
xlab = "p2", ylab = "t", zlab="b1",
main = "3D Surface plot for inflation factor b1"
)