Post date: Oct 18, 2017 7:57:42 PM
I am completing several summary analyses for the wing mapping project, including some aimed primarily at making graphics. Here is what I have:
1. Summary plots of PVE, PGE, and n-gamma. These are image plots and are made with the script mkPePlots.R in /uufs/chpc.utah.edu/common/home/u6000989/projects/lycaeides_wings/Gemma/combined_output/.
2. Plots showing posterior distributions for PVE, PGE and n-gamma. Files are in /uufs/chpc.utah.edu/common/home/u6000989/projects/lycaeides_wings/Gemma/combined_output/params/. First I had to run get*Dat.pl (no arguments) to get the MCMC stuff together. Plots and quantitative summaries of posterior density widths then come from plotPveDensities.R, pasted below:
N<-69
M<-500000
cs<-c(rep("orange",17),rep("darkgray",52))
files<-c("pve_est_all.txt","pve_est_AN.txt","pve_est_GNP.txt","pve_est_ID.txt","pve_est_JA.txt","pve_est_ME.txt","pve_est_MW.txt","pve_est_RI.txt","pve_est_SIN.txt","pve_est_SN.txt","pve_est_WA.txt","pve_est_YBG.txt")
files_o<-files[c(1:2,4:8,10:11,3,9,12)]
nms<-c("All","AN","GNP","ID","JA","ME","MW","RI","SIN","SN","WA","YBG")
nms_o<-nms[c(1:2,4:8,10:11,3,9,12)]
x<-vector("list",12)
## read in data, only once
for(i in 1:12){
x[[i]]<-matrix(scan(files_o[i],n=N*M,sep=" "),nrow=M,ncol=N,byrow=TRUE)
}
ymax<-c(15,11,8,8,rep(15,4),12,8,8, 8)
## make plot
pdf("pveDens.pdf",width=8,height=11)
par(mar=c(4,5,4,1))
par(mfrow=c(4,3))
for(i in 1:12){
hdr<-paste("(",letters[i],") ",nms_o[i],sep="")
plot(density(x[[i]][,1],from=0,to=1,na.rm=TRUE),ylim=c(0,ymax[i]),xlab="PVE",ylab="density",main=hdr,col=cs[1],lwd=0.9,cex.lab=1.4,cex.main=1.2)
for(j in 2:69){
lines(density(x[[i]][,j],from=0,to=1,na.rm=TRUE),col=cs[j],lwd=0.9)
}
if(i==1){legend(0.0,15,c("size","position"),lty=1,col=unique(cs),bty='n')}
abline(h=1,lty=2,lwd=1.2)
}
dev.off()
## HPD 90 interval widths
library(coda)
hpd90w<-matrix(NA,nrow=12,ncol=N)
for(i in 1:12){
for(j in 1:N){
o<-HPDinterval(as.mcmc(x[[i]][,j]),prob=0.9)
hpd90w[i,j]<-o[[2]]-o[[1]]
}
}
mean(as.vector(hpd90w[1,])) ## all
#[1] 0.2022477
mean(as.vector(hpd90w[2:9,])) ## taxa
#[1] 0.5674118
mean(as.vector(hpd90w[10:12,])) ## pops
#[1] 0.7639302
## pge
files<-c("pge_est_all.txt","pge_est_AN.txt","pge_est_GNP.txt","pge_est_ID.txt","pge_est_JA.txt","pge_est_ME.txt","pge_est_MW.txt","pge_est_RI.txt","pge_est_SIN.txt","pge_est_SN.txt","pge_est_WA.txt","pge_est_YBG.txt")
files_o<-files[c(1:2,4:8,10:11,3,9,12)]
xpge<-vector("list",12)
## read in data, only once
for(i in 1:12){
xpge[[i]]<-matrix(scan(files_o[i],n=N*M,sep=" "),nrow=M,ncol=N,byrow=TRUE)
}
ymax<-c(18,12,12,9,18,18,9,18,9,6,6,9)
## make plot
pdf("pgeDens.pdf",width=8,height=11)
par(mar=c(4,5,4,1))
par(mfrow=c(4,3))
for(i in 1:12){
hdr<-paste("(",letters[i],") ",nms_o[i],sep="")
plot(density(xpge[[i]][,1],from=0,to=1,na.rm=TRUE),ylim=c(0,ymax[i]),xlab="PGE",ylab="density",main=hdr,col=cs[1],lwd=0.9,cex.lab=1.4,cex.main=1.2)
for(j in 2:69){
lines(density(xpge[[i]][,j],from=0,to=1,na.rm=TRUE),col=cs[j],lwd=0.9)
}
if(i==1){legend(0.48,18,c("size","position"),lty=1,col=unique(cs),bty='n')}
abline(h=1,lty=2,lwd=1.2)
}
dev.off()
## HPD 90 interval widths
library(coda)
hpd90w<-matrix(NA,nrow=12,ncol=N)
for(i in 1:12){
for(j in 1:N){
o<-HPDinterval(as.mcmc(xpge[[i]][,j]),prob=0.9)
hpd90w[i,j]<-o[[2]]-o[[1]]
}
}
mean(as.vector(hpd90w[1,])) ## all
# [1] 0.3413153
mean(as.vector(hpd90w[2:9,])) ## taxa
# [1] 0.7660171
mean(as.vector(hpd90w[10:12,])) ## pops
# [1] 0.7983621
files<-c("ngam_est_all.txt","ngam_est_AN.txt","ngam_est_GNP.txt","ngam_est_ID.txt","ngam_est_JA.txt","ngam_est_ME.txt","ngam_est_MW.txt","ngam_est_RI.txt","ngam_est_SIN.txt","ngam_est_SN.txt","ngam_est_WA.txt","ngam_est_YBG.txt")
files_o<-files[c(1:2,4:8,10:11,3,9,12)]
## n-gamma
xng<-vector("list",12)
## read in data, only once
for(i in 1:12){
xng[[i]]<-matrix(scan(files_o[i],n=N*M,sep=" "),nrow=M,ncol=N,byrow=TRUE)
}
ymax<-c(1,0.8,0.15,0.55,1,.25,0.2,.25,.25,.1,.25,.75)
## make plot
pdf("ngDens.pdf",width=8,height=11)
par(mar=c(4,5,4,1))
par(mfrow=c(4,3))
for(i in 1:12){
hdr<-paste("(",letters[i],") ",nms_o[i],sep="")
plot(density(xng[[i]][,1],from=0,to=100,na.rm=TRUE),ylim=c(0,ymax[i]),xlab="no. loci",ylab="density",main=hdr,col=cs[1],lwd=0.9,cex.lab=1.4,cex.main=1.2)
for(j in 2:69){
lines(density(xng[[i]][,j],from=0,to=100,na.rm=TRUE),col=cs[j],lwd=0.9)
}
if(i==1){legend(45,0.95,c("size","position"),lty=1,col=unique(cs),bty='n')}
# abline(h=1,lty=2,lwd=1.2)
}
dev.off()
## HPD 90 interval widths
library(coda)
hpd90w<-matrix(NA,nrow=12,ncol=N)
for(i in 1:12){
for(j in 1:N){
o<-HPDinterval(as.mcmc(xng[[i]][,j]),prob=0.9)
hpd90w[i,j]<-o[[2]]-o[[1]]
}
}
mean(as.vector(hpd90w[1,])) ## all
#[1] 108.8406
mean(as.vector(hpd90w[2:9,])) ## taxa
#[1] 132.5851
mean(as.vector(hpd90w[10:12,])) ## pops
#[1] 155.3043
save(list=ls(),file="paramSum.rdat")
3.