-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathAgra_deltaKplot_fun.R
114 lines (104 loc) · 4.62 KB
/
Agra_deltaKplot_fun.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
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
##############################################################################/
##############################################################################/
#Delta-K method plotting functions
##############################################################################/
##############################################################################/
##############################################################################/
#Computing Delta-k####
##############################################################################/
chooseK<-function(str_out,nb_K,nb_rep) {
#'str_out': the summary of simulations file exported from STRUCTURE,
#before importing the file to R, replace the white space in the column
#header by underscores, you might also have to replace "?1" by "alpha"
#and remove double white space
#'nb_K': number of different K values considered
#'nb_rep':number of repetition for each K
datatable<-data.frame("K"=c(rep(1:nb_K,each=nb_rep)),"Ln(Pd)"=str_out[,4])
Lprim<-c(rep("NA",nb_rep))
for (i in ((nb_rep+1):(nb_K*nb_rep))) {
Lprim<-c(Lprim,str_out[i,4]-str_out[i-nb_rep,4])
}
datatable<-data.frame(datatable,as.numeric(Lprim))
Lsecond<-c(rep("NA",nb_rep))
for (i in (((2*nb_rep)+1):(nb_K*nb_rep))) {
Lsecond<-c(Lsecond,abs(datatable[i,3]-datatable[i-nb_rep,3]))
}
Lsecond<-c(Lsecond,rep("NA",nb_rep))
datatable<-data.frame(datatable,as.numeric(Lsecond))
reztable<-data.frame("K"=c(1:nb_K))
meanL<-c()
sdL<-c()
for (i in (1:nb_K)) {
meanL<-c(meanL,mean(datatable[datatable$K==i,2]))
sdL<-c(sdL,sd(datatable[datatable$K==i,2]))
}
reztable<-data.frame(reztable,meanL,sdL)
meanLprime<-c()
sdLprime<-c()
for (i in (1:nb_K)) {
meanLprime<-c(meanLprime,mean(as.numeric(datatable[datatable$K==i,3])))
sdLprime<-c(sdLprime,sd(datatable[datatable$K==i,3]))
}
reztable<-data.frame(reztable,meanLprime,sdLprime)
meanLsecond<-c()
sdLsecond<-c()
for (i in (1:nb_K)) {
meanLsecond<-c(meanLsecond,mean(as.numeric(datatable[datatable$K==i,4])))
sdLsecond<-c(sdLsecond,sd(datatable[datatable$K==i,4]))
}
reztable<-data.frame(reztable,meanLsecond,sdLsecond)
deltaK<-c()
for (i in (1:nb_K)) {
deltaK<-c(deltaK,reztable[reztable$K==i,6]/reztable[reztable$K==i,3])
}
reztable<-data.frame(reztable,deltaK)
return(reztable)
}
##############################################################################/
#function to plot variation of Delta K and Ln(P(X|K)) as a function of K####
##############################################################################/
plotdeltaK<-function(datadeltaK,nb_K,titre){
#'datadeltak': the output file of 'chooseK' function
#'nb_K': the number of different K considered
#'titre': the title of the plot you want to be displayed
plot(datadeltaK[1:(nb_K-2),8],type="b",pch=24,cex=2.5,lwd=4,lty=1,
col="transparent",bg="white",bty="n",ann=FALSE,axes=FALSE)
par(new=TRUE)
plot(datadeltaK[1:(nb_K-2),8],type="b",pch=24,bty="n",xaxt="n",
yaxt="n",ann=F,cex=2.5,lwd=4,lty=1)
axis(side=1,at=seq(1,13,1),lwd=3,font.axis=2)
axis(side=2,lwd=3,font.axis=2,las=1)
title(ylab="Delta K",font.lab=2,cex.lab=1.5)
par(new=TRUE)
plot(datadeltaK[1:(nb_K-2),2],type="b",pch=22,cex=2.5,lwd=4,lty=2,
col="grey50",bg="white",bty="n",xaxt="n",yaxt="n",ann=F)
axis(side=4,lwd=3,font.axis=2,col="grey50",las=1)
mtext("Ln(P(X|K))",side=4,line=4,font=2,cex=1.5,col="grey50")
title(main=titre,xlab="K",font.lab=2,cex.lab=1.5,cex.main=2)
}
#the same function using log(deltaK), just in order to see smaller variation
#of deltaK
plotlogdeltaK<-function(datadeltaK,nb_K,titre){
#'datadeltak': the output file of 'chooseK' function
#'nb_K': the number of different K considered
#'titre': the title of the plot you want to be displayed
op<-par(pty="s")
plot(log(datadeltaK[1:(nb_K-2),8]+1),type="b",pch=24,cex=2.5,lwd=4,
lty=1,col="transparent",bg="white",bty="n",ann=F)
par(new=TRUE)
plot(log(datadeltaK[1:(nb_K-2),8]+1),type="b",pch=24,bty="n",xaxt="n",
yaxt="n",ann=F,cex=2.5,lwd=4,lty=1)
axis(side=1,at=seq(1,13,1),lwd=3,font.axis=2)
axis(side=2,lwd=3,font.axis=2)
title(ylab="Ln(Delta K+1)",font.lab=2,cex.lab=1.5)
par(new=TRUE)
plot(datadeltaK[1:(nb_K-2),2],type="b",pch=22,cex=2.5,lwd=4,lty=2,
col="grey50",bg="white",bty="n",xaxt="n",yaxt="n",ann=F)
axis(side=4,lwd=3,font.axis=2,col="grey50")
mtext("Ln(P(X|K))", side=4, line=4,font=2,cex=1,col="grey50")
title(main=titre,xlab="K",font.lab=2,cex.lab=1.5,cex.main=2)
par(op)
}
##############################################################################/
#END
##############################################################################/