Post date: May 23, 2014 9:44:3 PM
The MCMC runs for cross-validation all finished. The results are summarized in a series of files in projects/lycaeides_hostplant/melGemma/results/, *Predict*.txt. These include family and individual cross-validation. I used two methods to summarize cross-validation. For survival in the plant x population treatments I computed the area under the ROC curve (auc), and for all other cases I computed the cross-validated r-squared (q-squared). The code for these analyses and the predictive plot is in projects/lycaeides_hostplant/melGemma/results/plotPredictive.R and pasted below. The main result is that we have little to know predictive power in most cases and particularly for the combined treatments (for survival this could in part reflect working with residuals). The main exception is survival in GLA on Ac or Ms, but still our predictive power is limited.
library(ROCR)
################ survival probit #################################
## read cross-validataion results
x<-vector('list',4)
x[[1]]<-read.table("results/survPredictglaTrtAc.txt",header=TRUE)
x[[2]]<-read.table("results/survPredictglaTrtMs.txt",header=TRUE)
x[[3]]<-read.table("results/survPredictslaTrtAc.txt",header=TRUE)
x[[4]]<-read.table("results/survPredictslaTrtMs.txt",header=TRUE)
## family cross-validataion
xfam<-vector('list',4)
xfam[[1]]<-read.table("results/survPredictFamglaTrtAc.txt",header=TRUE)
xfam[[2]]<-read.table("results/survPredictFamglaTrtMs.txt",header=TRUE)
xfam[[3]]<-read.table("results/survPredictFamslaTrtAc.txt",header=TRUE)
xfam[[4]]<-read.table("results/survPredictFamslaTrtMs.txt",header=TRUE)
N<-numeric(4)
for(i in 1:4){
N[i]<-dim(x[[i]])[1]
}
nr<-dim(x[[1]])[2]-1
auc<-matrix(NA,nrow=10,ncol=8) ## area under ROC curve, reps by experiments
for(ex in 1:4){
for(i in 1:10){
j<-which(is.na(x[[ex]][,(i+1)])==FALSE)
obs<-x[[ex]][j,1]
est<-x[[ex]][j,(i+1)]
pred<-prediction(est,obs)
stats_auc<-performance(pred,'auc')
auc[i,(ex - 1) * 2 + 1]<-as.numeric(stats_auc@y.values)
j<-which(is.na(xfam[[ex]][,(i+1)])==FALSE)
obs<-xfam[[ex]][j,1]
est<-xfam[[ex]][j,(i+1)]
pred<-prediction(est,obs)
stats_auc<-performance(pred,'auc')
auc[i,(ex - 1) * 2 + 2]<-as.numeric(stats_auc@y.values)
}
}
pdf("predictive.pdf",width=10,height=10)
par(mfrow=c(2,2))
par(mar=c(4,5.5,2,0.5))
boxplot(auc,col=rep(c("white","gray"),4),ylim=c(0,1),axes=FALSE,ylab="area under the curve (auc)",xlab="",cex.lab=1.9)
axis(1,at=seq(1.5,7.5,2),labels=c("GLAxAc","GLAxMs","SLAxAc","SLAxMs"),cex.axis=1.3)
axis(2,at=seq(0,1,0.2),cex.axis=1.3)
box()
abline(h=0.5,lty=3,lwd=2)
par(xpd=TRUE)
text(-0.7,1.1,"(a) survival [plant x pop.]",adj=0,cex=1.8)
par(xpd=FALSE)
####################### weight ############################
## read cross-validataion results
x<-vector('list',4)
x[[1]]<-read.table("results/wgtPredictglaTrtAc.txt",header=TRUE)
x[[2]]<-read.table("results/wgtPredictglaTrtMs.txt",header=TRUE)
x[[3]]<-read.table("results/wgtPredictslaTrtAc.txt",header=TRUE)
x[[4]]<-read.table("results/wgtPredictslaTrtMs.txt",header=TRUE)
## family cross-validataion
xfam<-vector('list',4)
xfam[[1]]<-read.table("results/wgtFamPredictglaTrtAc.txt",header=TRUE)
xfam[[2]]<-read.table("results/wgtFamPredictglaTrtMs.txt",header=TRUE)
xfam[[3]]<-read.table("results/wgtFamPredictslaTrtAc.txt",header=TRUE)
xfam[[4]]<-read.table("results/wgtFamPredictslaTrtMs.txt",header=TRUE)
N<-numeric(4)
for(i in 1:4){
N[i]<-dim(x[[i]])[1]
}
nr<-dim(x[[1]])[2]-1
q2<-matrix(NA,nrow=10,ncol=8) ## cross-validation r2, reps by experiments
for(ex in 1:4){
for(i in 1:10){
j<-which(is.na(x[[ex]][,(i+1)])==FALSE)
obs<-x[[ex]][j,1]
est<-x[[ex]][j,(i+1)]
q2[i,(ex - 1) * 2 + 1]<-1-sum((obs-est)^2,na.rm=T)/sum((obs-mean(obs,na.rm=T))^2,na.rm=T)
j<-which(is.na(xfam[[ex]][,(i+1)])==FALSE)
obs<-xfam[[ex]][j,1]
est<-xfam[[ex]][j,(i+1)]
q2[i,(ex - 1) * 2 + 2]<-1-sum((obs-est)^2,na.rm=T)/sum((obs-mean(obs,na.rm=T))^2,na.rm=T)
}
}
boxplot(q2,col=rep(c("white","gray"),4),ylim=c(-0.5,1),axes=FALSE,ylab=expression(paste("cross-validated ",r^2," (",q^2,")")),xlab="",cex.lab=1.9)
axis(1,at=seq(1.5,7.5,2),labels=c("GLAxAc","GLAxMs","SLAxAc","SLAxMs"),cex.axis=1.3)
axis(2,at=seq(-0.5,1,0.5),cex.axis=1.3)
abline(h=0,lty=3,lwd=2)
box()
par(xpd=TRUE)
text(-0.7,1.15,"(b) weight [plant x pop.]",adj=0,cex=1.8)
par(xpd=FALSE)
####################### survival combined ############################
## read cross-validataion results
x<-vector('list',5)
x[[1]]<-read.table("results/survPredictCombcombGla.txt",header=TRUE)
x[[2]]<-read.table("results/survPredictCombcombSla.txt",header=TRUE)
x[[3]]<-read.table("results/survPredictCombcombAc.txt",header=TRUE)
x[[4]]<-read.table("results/survPredictCombcombMs.txt",header=TRUE)
x[[5]]<-read.table("results/survPredictCombcombAll.txt",header=TRUE)
## family cross-validataion
xfam<-vector('list',5)
xfam[[1]]<-read.table("results/survFamPredictCombcombGla.txt",header=TRUE)
xfam[[2]]<-read.table("results/survFamPredictCombcombSla.txt",header=TRUE)
xfam[[3]]<-read.table("results/survFamPredictCombcombAc.txt",header=TRUE)
xfam[[4]]<-read.table("results/survFamPredictCombcombMs.txt",header=TRUE)
xfam[[5]]<-read.table("results/survFamPredictCombcombAll.txt",header=TRUE)
N<-numeric(5)
for(i in 1:5){
N[i]<-dim(x[[i]])[1]
}
nr<-dim(x[[1]])[2]-1
q2<-matrix(NA,nrow=10,ncol=10) ## cross-validation r2, reps by experiments
for(ex in 1:5){
for(i in 1:10){
j<-which(is.na(x[[ex]][,(i+1)])==FALSE)
obs<-x[[ex]][j,1]
est<-x[[ex]][j,(i+1)]
q2[i,(ex - 1) * 2 + 1]<-1-sum((obs-est)^2,na.rm=T)/sum((obs-mean(obs,na.rm=T))^2,na.rm=T)
j<-which(is.na(xfam[[ex]][,(i+1)])==FALSE)
obs<-xfam[[ex]][j,1]
est<-xfam[[ex]][j,(i+1)]
q2[i,(ex - 1) * 2 + 2]<-1-sum((obs-est)^2,na.rm=T)/sum((obs-mean(obs,na.rm=T))^2,na.rm=T)
}
}
boxplot(q2,col=rep(c("white","gray"),5),ylim=c(-0.5,1),axes=FALSE,ylab=expression(paste("cross-validated ",r^2," (",q^2,")")),xlab="treatment",cex.lab=1.9)
axis(1,at=seq(1.5,9.5,2),labels=c("GLA","SLA","Ac","Ms","all"),cex.axis=1.3)
axis(2,at=seq(-0.5,1,0.5),cex.axis=1.3)
abline(h=0,lty=3,lwd=2)
box()
par(xpd=TRUE)
text(-0.9,1.15,"(c) survival [combined]",adj=0,cex=1.8)
par(xpd=FALSE)
####################### weight combined ############################
## read cross-validataion results
x<-vector('list',5)
x[[1]]<-read.table("results/wgtPredictCombcombGla.txt",header=TRUE)
x[[2]]<-read.table("results/wgtPredictCombcombSla.txt",header=TRUE)
x[[3]]<-read.table("results/wgtPredictCombcombAc.txt",header=TRUE)
x[[4]]<-read.table("results/wgtPredictCombcombMs.txt",header=TRUE)
x[[5]]<-read.table("results/wgtPredictCombcombAll.txt",header=TRUE)
## family cross-validataion
xfam<-vector('list',5)
xfam[[1]]<-read.table("results/wgtFamPredictCombcombGla.txt",header=TRUE)
xfam[[2]]<-read.table("results/wgtFamPredictCombcombSla.txt",header=TRUE)
xfam[[3]]<-read.table("results/wgtFamPredictCombcombAc.txt",header=TRUE)
xfam[[4]]<-read.table("results/wgtFamPredictCombcombMs.txt",header=TRUE)
xfam[[5]]<-read.table("results/wgtFamPredictCombcombAll.txt",header=TRUE)
N<-numeric(5)
for(i in 1:5){
N[i]<-dim(x[[i]])[1]
}
nr<-dim(x[[1]])[2]-1
q2<-matrix(NA,nrow=10,ncol=10) ## cross-validation r2, reps by experiments
for(ex in 1:5){
for(i in 1:10){
j<-which(is.na(x[[ex]][,(i+1)])==FALSE)
obs<-x[[ex]][j,1]
est<-x[[ex]][j,(i+1)]
q2[i,(ex - 1) * 2 + 1]<-1-sum((obs-est)^2,na.rm=T)/sum((obs-mean(obs,na.rm=T))^2,na.rm=T)
j<-which(is.na(xfam[[ex]][,(i+1)])==FALSE)
obs<-xfam[[ex]][j,1]
est<-xfam[[ex]][j,(i+1)]
q2[i,(ex - 1) * 2 + 2]<-1-sum((obs-est)^2,na.rm=T)/sum((obs-mean(obs,na.rm=T))^2,na.rm=T)
}
}
boxplot(q2,col=rep(c("white","gray"),5),ylim=c(-0.5,1),axes=FALSE,ylab=expression(paste("cross-validated ",r^2," (",q^2,")")),xlab="treatment",cex.lab=1.9)
axis(1,at=seq(1.5,9.5,2),labels=c("GLA","SLA","Ac","Ms","all"),cex.axis=1.3)
axis(2,at=seq(-0.5,1,0.5),cex.axis=1.3)
abline(h=0,lty=3,lwd=2)
box()
par(xpd=TRUE)
text(-0.9,1.15,"(d) weight [combined]",adj=0,cex=1.8)
par(xpd=FALSE)
dev.off()