-
Notifications
You must be signed in to change notification settings - Fork 0
/
Int3.sampler
50 lines (45 loc) · 1.51 KB
/
Int3.sampler
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
Int3.sampler <- function(Sample = 1, database = Pop, nchr = 5, iter = 1000) {
P.values<<- list()
for (i in c(1:iter)) {
p.outcome<- c()
Draw <- DrawIncCSL(database, y = Sample)
for (p in 2:(nchr+1)) { #the Chr columns
for (q in p:(nchr+1)) {
for (r in q:(nchr+1)) {
if (p < q & q < r) {
summer <<- summary(aov(Draw[[trait]]~ as.factor(Draw[[p]]) *as.factor(Draw[[q]])as.factor(Draw[[r]])))
if (length(summer[[1]][["Pr(>F)"]]) < 4) {
p.outcome<- c(p.outcome, NA)
} else {
p.outcome<- c(p.outcome, -log(summer[[1]][["Pr(>F)"]][[3]],10))
}
}}
}
Val <- list(p.outcome)
P.values<<- c(P.values, Val)
Sigs <- as.data.frame((do.call(cbind, P.values)))
}
Sigs <- as.data.frame(t(Sigs))
nna<- c()
p95 <- c()
p95_BF <- c()
meanp<- c()
stdev<- c()
for (i in 1:ncol(Sigs)) {
nna<- c(nna, sum(is.na(Sigs[[i]])))
stdev<- c(stdev, sd(Sigs[[i]], na.rm=TRUE))
meanp<- c(meanp, mean(Sigs[[i]], na.rm=TRUE))
Sigs[[i]][is.na(Sigs[[i]])] <- 0
p95 <- c(p95, sum(Sigs[[i]]>-log(0.05,10))/nrow(Sigs)*100)
p95_BF <- <- c(p95, sum(Sigs[[i]]>-log(0.05/ncol(Sigs),10))/nrow(Sigs)*100)
}
Results <- data.frame(t(Sigs))
Results$STDEV<- stdev
Results$P95 <- p95
Results$P95_BF<- p95_BF
Results$MEANP<- meanp
Results$NNA<- nna
Results <- data.frame(t(Results))
return(Results[(nrow(Results)-4):(nrow(Results)),])
}
# returns a file with the averaged LOD-values, the standard deviation, the number of ANOVA-outcomes being significant (p < 0.05, uncorrected), the number of outcomes being significant ((p < 0.05, Bonferroni).