# TODO: Add comment # # Author: E.Korsching Mar 7, 2013 ############################################################################### estimation.ref.size.impact <- function(x, refset=c(1,5,2,3,12,13), testset=c(4,6:11,14:16), method="pearson", anzSamplings=1000, exp.name="froh.ck", plot1=F, fittype="both", bass=6, subfolder=""){ # Scale down the number of reference factors (without touching the test factors) to estimate the distance to noise # x: data set # j: number of reference factors # i: actual ref set size - linear downsizing from the end # l: loop/result counter # anzSamplings: 0 : only original, 1..many: original + bootstrap examples # exp.name: a suffix to the x data name, will be numbered (for each shortening) # plot1: plot all regression lines # fittype: "loess" [more exact] or "supsmu" [more smooth] (with bass parameter) or "both" # bass: smoothing parameter : less smooth - 0..10 - very smooth # method: pearson, spearman or euclidean # subfolder: e.g. "/results/e1/" will be appended to work folder path #ini j <- length(refset) # number of reference factors i <- j #start condition x.name <- deparse(substitute(x)) cat("\n x.name",x.name) erg1 <- data.frame(matrix(0, j-1, 9)) names(erg1) <- c("n.ref","real.ssqg","real.sum.ssqg","boot.ssqg.min","boot.ssqg.max","boot.ssqg.median","boot.min.sum.ssqg","boot.max.sum.ssqg","boot.median.sum.ssqg") #production run l <- 0 #counter start value while(i > 1){ l <- l+1 #loop counter, number of results txt.0 <- paste(x.name,".",exp.name,l, sep="") assign(txt.0, enumeration.search.dep.boots.f90(x, anzSamplings=anzSamplings, refset=refset[1:i], testset=testset, method=method, sampleonly=F), envir=.GlobalEnv) txt.1 <- get(txt.0) #make var #collect some results tmp0 <- anzSamplings+1 tmp1 <- which.min(txt.1[[2]][2:tmp0,1]) +1 #+1 because we start with 2 and therefore get a relative position ! tmp2 <- which.max(txt.1[[2]][2:tmp0,1]) +1 erg1[l,1] <- i #n.ref erg1[l,2] <- txt.1[[2]][1,1]/i #real.ssqg all values - erg1[,2:9] - normalized by i : number of factors erg1[l,3] <- txt.1[[2]][1,5]/i #real.sum.ssqg erg1[l,4] <- txt.1[[2]][tmp1,1]/i #boot.ssqg.min erg1[l,5] <- txt.1[[2]][tmp2,1]/i #boot.ssqg.max erg1[l,6] <- median(txt.1[[2]][2:tmp0,1])/i #boot.ssqg.median erg1[l,7] <- txt.1[[2]][tmp1,5]/i #boot.min.sum.ssqg erg1[l,8] <- txt.1[[2]][tmp2,5]/i #boot.max.sum.ssqg erg1[l,9] <- median(txt.1[[2]][2:tmp0,5])/i #boot.median.sum.ssqg, approximation here ! check if appropriate for normalisation i <- i-1 #shortening of ref } #plot -- regression line graphs if(plot1){ pdf(file=paste(getwd(),subfolder,x.name,".",exp.name,".qualityRef",".",format(Sys.time(), "%Y%m%d%H%M"),".pdf", sep=""), width=11, height=7, onefile=T, title=paste(x.name,".",exp.name,"qualityRef", sep=""), pointsize=12) for(ii in 1:l){ txt.0 <- paste(x.name,".",exp.name,ii,sep="") txt.1 <- get( txt.0 ) #make var tmp1 <- length(testset) tmp0 <- nrow(txt.1[[2]]) # plot original graph plot.sample.search.dependencies(x=txt.1[[3]][1:tmp1,,drop=F], o=txt.1[[1]][1,(tmp1+1):2], ylim=c(-0.7,0.7), y.at=NULL, ylab="correlation", text=T, description=txt.0, cex=0.8, mfrow=c(2,3) ) # supplement with bootstrap values plot(x=1,y=1, type="n", xlim=c(0,10), ylim=c(0,10), axes=F, xlab="", ylab="") text(x=1,y=9, labels=paste("min. : ", round( min(txt.1[[2]][2:tmp0,1]) ,digits=1), ",\n", "max. : ", round( max(txt.1[[2]][2:tmp0,1]) ,digits=1), ",\n", "median : ", round( median(txt.1[[2]][2:tmp0,1]) ,digits=1), sep=""), adj=0, cex=1.2) title(main="bootstrap ssqg") } dev.off() } # plot normalized ssqg pdf(file=paste(getwd(),subfolder,x.name,".",exp.name,".main.qualityRef",".",format(Sys.time(), "%Y%m%d%H%M"),".pdf", sep=""), width=11, height=7, onefile=T, title=paste(x.name,".",exp.name,"main.qualityRef", sep=""), pointsize=12) tmp0 <- c(0,j) tmp1 <- c(min(erg1[,c(2,4:6)]),max(erg1[,c(2,4:6)])) plot(x=erg1[,1],y=erg1[,2], type="p", pch=19, ylim=tmp1, xlab="number of reference factors", ylab="gobal ssq") #ssqg points(x=erg1[,1],y=erg1[,4], pch=2) #min sampling points(x=erg1[,1],y=erg1[,5], pch=6) #max sampling points(x=erg1[,1],y=erg1[,6], pch=5) #median sampling if(fittype=="supsmu"){ lines(supsmu(x=erg1[,1],y=erg1[,2], bass=bass), lty=2, col="blue") lines(supsmu(x=erg1[,1],y=erg1[,6], bass=bass), lty=2, col="red") } if(fittype=="loess"){ lines(erg1[,1], predict(loess(erg1[,2]~erg1[,1])), lty=2, col="blue") #lwd=2 lines(erg1[,1], predict(loess(erg1[,6]~erg1[,1])), lty=2, col="red") } if(fittype=="both"){ lines(supsmu(x=erg1[,1],y=erg1[,2], bass=bass), lty=2, col="blue") lines(supsmu(x=erg1[,1],y=erg1[,6], bass=bass), lty=2, col="red") lines(erg1[,1], predict(loess(erg1[,2]~erg1[,1])), lty=2, col="blue") lines(erg1[,1], predict(loess(erg1[,6]~erg1[,1])), lty=2, col="red") } legend(x=(tmp0[2]-((tmp0[2]-tmp0[1])*0.65)), y=(tmp1[2]-((tmp1[2]-tmp1[1])*0.1)), legend=c("ssqg","sampling, max ssqg","sampling, median ssqg","sampling, min ssqg"), col=c("black","black","black","black"), pch=c(19,6,5,2), cex=0.7) legend(x=(tmp0[2]-((tmp0[2]-tmp0[1])*0.65)), y=(tmp1[2]-((tmp1[2]-tmp1[1])*0.005)), legend=c("fit of ssqg","fit of sampling median"), col=c("blue","red"), lty=c(2,2), cex=0.7) title(main=paste("Reference factor downsizing \n Basic global ssq versus sampling global ssq (normalized by n ref)", sep="")) dev.off() return(erg1) } #save.image() #a <- date() #gfroRef <- estimation.ref.size.impact(x=fro.orig.1a, refset=c(1,2,16,17,18,19), testset=c(11,12,13,10,4,20,5,15), # anzSamplings=1000, exp.name="testref", plot1=T, fittype="both", bass=4, subfolder="/results/frohwitter/") #date() # 1.5min