-
Notifications
You must be signed in to change notification settings - Fork 1
/
ART.r
68 lines (64 loc) · 1.83 KB
/
ART.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
# Sample code for ART, ART-A, and RTP
# https://doi.org/10.1101/667238
#
library(MASS)
library(mvtnorm)
# RTP via Equ. 2
RTP <- function(Z, k, L) { # Z <- sum(-log(P[1:k]))
integrate(function(x,y,m,n) 1-pgamma(log(qbeta(x,m+1,n-m))*m+y,m),0,1,Z,k,L)$va
}
ART.A <- function(P, k, L) { # "ART-A"
wgt <- rep(1,k)
z <- P
z[1] <- ( 1 - P[1] )^L
for(j in 2:L) z[j] <- ((1-P[j]) / (1-P[j-1]))^((L-(j-1)))
p <- (1-z)[1:k]
k = length(p)
sumZ <- rep(0, k)
y <- qnorm(p)
z <- y
gz <- z[1] * wgt[1]
sumZ[1] <- gz
for(i in 2:k) {
gz <- p[ 1 : i ]
for(j in 1:i) gz[j] <- z[j] * wgt[j]
sumZ[i] <- sum(gz)
}
Lo = diag(k); Lo[lower.tri(Lo)] <- 1
pSg <- Lo %*% diag(wgt[1:k]^2) %*% t(Lo)
pCr <- cov2cor(pSg)
sZ <- sumZ
for(i in 1:k) {
sZ[i] <- sumZ[i] / sqrt(diag(pSg))[i]
}
ppZ <- pmvnorm(lower = rep(-Inf,k), upper = rep(max(sZ), k), sigma = pCr)[1]
c(ppZ, which.min(sZ))
}
ART <- function(lW, Pk, k, L) {
d = (k-1)*(digamma(L+1) - digamma(k))
ak = (k-1)*log(Pk) - lW + qgamma(1-pbeta(Pk, k, L-k+1), shape=d)
1 - pgamma(ak, shape=k+d-1)
}
## read in summary statistics and LD matrix
d <- read.csv("Z_score.csv", header=TRUE, as.is=TRUE)
S <- read.csv("LD_matrix.csv", header=FALSE, as.is=TRUE)
## the number of statistics
Z <- d$Z
L <- length(Z)
## transform Z-scores to chi-squares
y <- Z^2
## calculate TQ-statistic
yy <- y %*% rep(1, L)
## eigen decomposition of the LD matrix
ee <- eigen(S); eivec <- ee$vectors; eigva <- ee$values
pc <- eivec %*% diag(sqrt(1/eigva)) %*% t(eivec)
## calculate decorreated statistics
x <- (Z %*% pc)^2
k <- round(0.5*L)
#k <- L
px <- 1-pchisq(x, df=1)
P <- sort(px)
P.rtp <- RTP(sum(-log(P[1:k])), k, L)
P.art <- ART(sum(log(P[1:(k-1)])), P[k], k, L)
P.arta <- ART.A(P, k, L)
cat(P.art, P.rtp, P.arta[1], P.arta[2], "\n")