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
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
|
# read filename from the command line
args <- commandArgs(trailingOnly=TRUE)
filename <- args[1]
outfile1 <- args[2]
outfile2 <- args[3]
# read file
LOD <- as.matrix(read.table(filename,header=FALSE))
N <- ncol(LOD) # number of runs + 1
L <- nrow(LOD) # number of splits
LOD <- matrix(LOD,ncol=N)
# sort LOD by last column
O <- order(LOD[,N])
LOD <- matrix(LOD[O,],ncol=N)
# create PP table
PP <- (10**LOD)/(1+10**LOD)
# split last column out of tables
aveLOD <- as.vector(LOD[,N])
avePP <- as.vector( PP[,N])
LOD <- matrix(LOD[,1:N-1],ncol=N-1)
PP <- matrix(PP[,1:N-1],ncol=N-1)
lodToPP <- function(x) {y<-exp(x); y/(1+y)}
ppToLod <- function(x) {log10(x/(1-x))}
#-------------- Plot 1 -------------------
svg(file=outfile1,height=3,width=7)
par(mar=c(4, 4, 0, 4) + 0.1)
plot(avePP,xlab="Split",ylab="PP",type="n",ylim=c(0,1),xaxt="n")
axis(side=1,at=1:L,1:L)
# plot results from another analysis for comparison
# lines(avePP_other,col=hsv(0.63,0.4,1),lwd=2)
ptcolor <- hsv(0, 1, 0, 0.125)
for(i in 1:(N-1)) {
xs <- 1:L
xs <- xs + rnorm(L, 0, N/500)
points(xs, PP[,i], col=ptcolor, pch=".",cex=4)
}
# Plot the bars representing the range for each split
minPP <- apply(PP, 1, min, na.rm=TRUE)
maxPP <- apply(PP, 1, max, na.rm=TRUE)
for(i in 1:L) {
X <- sort(PP[i,])
LX <- length(X)
for(j in 1:(LX-1))
{
lines(c(i,i), c(X[j],X[j+1]), lwd=min(j,LX-j), col=hsv(1,0.5,0,0.5))
}
}
lodticks <- c(-3,-2,-1,0,1,2,3)
axis(side=4, at=lodToPP(lodticks),labels = lodticks)
mtext(side=4,line=3,expression(log[10](PP/(1-PP))))
ppmin <- min(minPP)
ppmax <- max(maxPP)
# Plot the estimate for each split, and connect the dots
lines(avePP,col=hsv(1,1,1),lwd=2)
#-------------- Plot 2 -------------------
svg(file=outfile2,height=3,width=7)
lodmax <- log10(ppmax/(1-ppmax));
lodmin <- log10(ppmin/(1-ppmin));
ymax <- lodmin + 1.02*(lodmax - lodmin)
ymin <- lodmax + 1.02*(lodmin - lodmax)
ymax <- max(2.1,ymax)
ymin <- min(-2.1,ymin)
par(mar=c(4, 4, 0, 4) + 0.1)
plot(aveLOD,xlab="Split",ylab=expression(log[10](PP/(1-PP))), ylim=c(ymin,ymax),type="n",xaxt="n")
axis(side=1,at=1:L,1:L)
lodToPP <- function(x) {y=exp(x);y/(1+y)}
ppToLod <- function(x) {log10(x/(1-x))}
#par(mar=c(4,1,1,1))
ppticks =c(0.001,0.01,0.1,0.5,0.9,0.99,0.999)
pplabels= c("0.001","0.01","0.1","0.5","0.9","0.99","0.999")
axis(side=4, at=ppToLod(ppticks),labels = pplabels)
mtext(side=4,line=3,'PP')
# plot results from another analysis for comparison
# lines(LOD,col=hsv(0.63,0.4,1),lwd=2)
for(i in 1:(N-1)) {
xs <- 1:L
xs <- xs + rnorm(L, 0, N/500)
points(xs, LOD[,i], col=ptcolor, pch=".",cex=4)
}
# Plot the bars representing the range for each split
minLOD <- apply(LOD, 1, min, na.rm=TRUE)
maxLOD <- apply(LOD, 1, max, na.rm=TRUE)
for(i in 1:L) {
X <- sort(LOD[i,])
LX <- length(X)
for(j in 1:(LX-1))
{
lines(c(i,i), c(X[j],X[j+1]), lwd=min(j,LX-j), col=hsv(1,0.5,0,0.5))
}
}
# Plot the estimate for each split, and connect the dots
lines(aveLOD,col=hsv(1,1,1),lwd=2)
|
