Post date: May 12, 2014 8:55:14 PM
I decided to compute the area under the ROC curve to assess prediction performance of the BSLMM with cross-validation for the survival results (I will use root mean square error for the weight measurements). These are the same methods used by Zhou et al. 2013. The AUC calculations are done in R with the RORC package. AUC above 0.5 means that there is some improvement relative to a null model with just an intercept. We see some predictive ability for GLA on Ms and Ac, even with family cross-validation (see projects/lycaeides_hostplant/melGemma/auc.pdf, gray is family, white is individual). The R code I used to generate this is projects/lycaeides_hostplant/melGemma/scripts/plotAUC.R (pasted below). I want to make a four pane figure that includes this figure, adult weight (RMSE) and RMSE for survival and weight in the combined experiments. I still need to run the cross-validation for the latter.
library(RORC)
## 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("auc.pdf",width=5,height=5)
par(mar=c(5.5,5.5,0.5,0.5))
boxplot(auc,col=rep(c("white","gray"),4),ylim=c(0,1),axes=FALSE,ylab="area under the curve (auc)",xlab="experiment",cex.lab=1.6)
axis(1,at=seq(1.5,7.5,2),labels=c("GLAxAc","GLAxMs","SLAxAc","SLAxMs"),cex.axis=1.1)
axis(2,at=seq(0,1,0.2),cex.axis=1.1)
box()
abline(h=0.5,lty=3,lwd=2)
dev.off()