-
Notifications
You must be signed in to change notification settings - Fork 0
/
Influence_of_Dispersal.R
65 lines (48 loc) · 2.43 KB
/
Influence_of_Dispersal.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
########
L <- mapply(rbind, bl1_combine_positions, bl2_combine_positions, bl3_combine_positions,
il1_combine_positions, il2_combine_positions, il3_combine_positions, SIMPLIFY=FALSE)
M <- mapply(rbind, bm1_combine_positions, bm2_combine_positions, bm3_combine_positions,
im1_combine_positions, im2_combine_positions, im3_combine_positions, SIMPLIFY=FALSE)
H <- mapply(rbind, bh1_combine_positions, bh2_combine_positions, bh3_combine_positions,
ih1_combine_positions, ih2_combine_positions, ih3_combine_positions, SIMPLIFY=FALSE)
plotL <- PlotRates(L, cl, ylab = "Rate")
plotM <- PlotRates(M, cm, xlab = "Threshold (p.TBI)")
plotH <- PlotRates(H, ch)
plotlegend <- PlotRates(L, cl, leg = "right")
legend <- get_legend(plotlegend)
pgrid <- plot_grid(plotL, plotM, plotH, ncol = 3)
pdisp <- plot_grid(pgrid, legend, rel_widths = c(1, .1))
pdisp
sapply(L, function(x) mean(x$FNR))[11]
sapply(M, function(x) mean(x$FNR))[11]
sapply(H, function(x) mean(x$FNR))[11]
sapply(L, function(x) mean(x$FPR))[11]
sapply(M, function(x) mean(x$FPR))[11]
sapply(H, function(x) mean(x$FPR))[11]
alpha <- c(0.0001, 0.00025, 0.0005, 0.00075, 0.001, 0.0025, 0.005, 0.0075, 0.01, 0.025, 0.05, 0.075, 0.1)
ml <- sapply(L, function(x) mean(x$FPR))
mm <- sapply(M, function(x) mean(x$FPR))
mh <- sapply(H, function(x) mean(x$FPR))
plot(alpha, ml, col = "red", pch=19, ylim = c(0,0.052))
points(alpha, mm, col = "blue", pch=19)
points(alpha, mh, col = "green", pch=19)
# Lval <- sapply(L, function(x) mean(x$FPR))
# Mval <- sapply(M, function(x) mean(x$FPR))
# Hval <- sapply(H, function(x) mean(x$FPR))
#
# kwtlist <- list(L[[9]]$FPR, M[[9]]$FPR, H[[9]]$FPR)
# kruskal.test(kwtlist)
# vals9 <- data.frame(FPR=c(L[[9]]$FPR, M[[9]]$FPR, H[[9]]$FPR), FNR=c(L[[9]]$FNR, M[[9]]$FNR, H[[9]]$FNR),
# group=rep(factor(c("L", "M", "H")), each=length(L[[9]]$FPR)))
# summary(aov(FPR~group, data=vals9))
#
# chisq.test(Lval, Mval, simulate.p.value=TRUE, B=1000000)
# do.call('rbind', L)
#
# vals9cont <- data.frame(FPR = unlist(L, cm, ch),
# disp = rep(factor(c("L", "M", "H")), each=length(cl[[1]])*13),
# threshold = rep(c(0.0001, 0.00025, 0.0005, 0.00075,
# 0.001, 0.0025, 0.005, 0.0075,
# 0.01, 0.025, 0.05, 0.075,
# 0.1), each=180))
# anova(lm(FPR~disp, data=vals9cont))