# TODO: Add comment # # Author: E.Korsching Jul 16, 2012, 2021 ############################################################################### perm.testing.1 <- function(x, ref, test, method="pearson", dataname="path/file", ylim=c(-1,1), colorG="#993300"){ ### testing of a given grouping + best in data grouping # x: data.frame # ref: c(4,5,6,7,8,10) (e.g. CK1 CK10 CK19 CK8.18 CK5 CK14) # test: c(1:3,9,11,14,17:20) # method: pearson, spearman or euclidean # dataname: relative "path/file" name # ylim: y axis range in regression plot # ini anzSamplings <- 1000 len.ref <- length(ref) len.test <- length(test) x <- x[,c(test,ref),drop=F] #order and filter x.names <- names(x) #new names test <- c(1:len.test) #new test ref <- c((len.test+1):(len.test+len.ref)) #new ref cat("\n 1") ## input order # 3s (6/10 in total approx. 5h) , check var=0 ! x.hy1.f90 <- enumeration.search.dep.boots.f90(x=x, anzSamplings=0, refset=ref, testset=test, method=method, sampleonly=F, output="boot.raw.col", ssqflag=1) cat("\n 2") # bootstrap - 28min x.hy1.b1000f90 <- enumeration.search.dep.boots.f90(x=x, anzSamplings=anzSamplings, refset=ref, testset=test, method=method, sampleonly=F, output="boot.raw.col", ssqflag=1) cat("\n 3") # extrema ssq.given.situation.min <- x.hy1.b1000f90[[2]][1,1] sample.min.row <- which.min(x.hy1.b1000f90[[2]][,1]) #642 sample.min.value <- min(x.hy1.b1000f90[[2]][,1], na.rm=T) #151.7538 sample.max.value <- max(x.hy1.b1000f90[[2]][,1], na.rm=T) ssq.sample <- x.hy1.b1000f90[[2]][ ,1] number.better <- sum(ssq.given.situation.min>ssq.sample, na.rm=T) sample.p.value <- number.better/anzSamplings #45/1000 = 0.045 cat("\n 4") # regression plot pdf(file=paste("results/",dataname,".regression.",format(Sys.time(), "%Y%m%d%H%M"),".pdf", sep=""), width=11, height=7, onefile=T, title=paste(dataname,"regression"), pointsize=12) plot.sample.search.dependencies(x=list( list(x.hy1.f90[[3]][1:len.test,,drop=F], x.hy1.f90[[1]][1,2:(len.test+1)], colorG) ), raw=NULL, ref.size.range=c(1,4), ylim=ylim, y.at=NULL, ylab="correlation", text=T, description="best", cex.t=0.6, cex.rt=1, mfrow=c(2,3), mar.b.f=0, ssqflag=1 ) plot.sample.search.dependencies(x=list( list(x.hy1.b1000f90[[3]][((sample.min.row*len.test)-(len.test-1)):(sample.min.row*len.test),,drop=F], x.hy1.b1000f90[[1]][sample.min.row,2:(len.test+1)], #11:2 colorG) ), raw=NULL, ref.size.range=c(1,4), ylim=ylim, y.at=NULL, ylab="correlation", text=T, description="bootstrap best", cex.t=0.6, cex.rt=1, mfrow=c(2,3), mar.b.f=0, ssqflag=1 ) dev.off() cat("\n 5") ## base distribution # 1s x.hy1.f90dist <- enumeration.search.dep.boots.f90.base.dist(x=x, refset=ref, testset=test, method=method, zz=NULL, prox.only=F, output="boot.raw.col", ssqflag=1) # all optimal orders - 3h55min x.hy1.extf90 <- enumeration.search.dep.boots.GOnly.extf90(x=x, gRef=ref, gTest=test, method=method, ssqflag=1) cat("\n 6") # extrema ssq.given.situation.max <- max(x.hy1.f90dist[[1]], na.rm=T) comb.min.row <- which.min(x.hy1.extf90[,1]) #406 comb.min.value <- min(x.hy1.extf90[,1], na.rm=T) #210.2306 comb.max.value <- max(x.hy1.extf90[,1], na.rm=T) orders.min <- x.hy1.extf90[comb.min.row, ] #order: best ssq, best test order of 'selected test', selected test, selected ref cat("\n 7") # histogram plots pdf(file=paste("results/",dataname,".histo.",format(Sys.time(), "%Y%m%d%H%M"),".pdf", sep=""), width=11, height=7, onefile=T, title=paste(dataname,"histo"), pointsize=12) hist.plot.simple(x=x.hy1.f90dist[[1]], bin.num=F, bin.size=F, xlab="ssqg range", bar.width=1, offset=0, digits=1, col=c("red","green","blue"), cex=1, h.title=paste(dataname,"distribution of all ssqg in the selected set") ) hist.plot.simple(x=x.hy1.b1000f90[[2]][,1], bin.num=F, bin.size=F, xlab="", bar.width=1, offset=0, digits=1, col=c("red","green","blue"), cex=1, h.title=paste(dataname,"distribution of bootstrap ssqg in the selected set, p=", round( sum(x.hy1.b1000f90[[2]][,1]ssq2.sample, na.rm=T) sample2.p.value <- number2.better/anzSamplings cat("\n 10") # regression plots zwei Mal identisch ? was passiert ? pdf(file=paste(getwd(),"/results/",dataname,".regression.best.",format(Sys.time(), "%Y%m%d%H%M"),".pdf", sep=""), width=11, height=7, onefile=T, title=paste(dataname,"regression.best"), pointsize=12) plot.sample.search.dependencies(x=list( list(x.hy1a.f90[[3]][1:len.test,,drop=F], x.hy1a.f90[[1]][1,2:(len.test+1)], colorG) ), raw=NULL, ref.size.range=c(1,4), ylim=ylim, y.at=NULL, ylab="correlation", text=T, description="best of best combination", cex.t=0.6, cex.rt=1, mfrow=c(2,3), mar.b.f=0, ssqflag=1) plot.sample.search.dependencies(x=list( list(x.hy1a.b1000f90[[3]][((sample2.min.row*len.test)-(len.test-1)):(sample2.min.row*len.test),,drop=F], x.hy1a.b1000f90[[1]][sample2.min.row,2:(len.test+1)], colorG) ), #11:2 raw=NULL, ref.size.range=c(1,4), ylim=ylim, y.at=NULL, ylab="correlation", text=T, description="bootstrap best of best combination", cex.t=0.6, cex.rt=1, mfrow=c(2,3), mar.b.f=0, ssqflag=1) dev.off() cat("\n 11") # return all information cat("\n\n") cat("\n #### data: ", dataname, " #### date: ", date(), " ####") cat("\n number of test permutations: ", factorial(len.test), " , number of ref combinations: ", choose(len.ref+len.test, len.ref)) cat("\n number of samplings: ", anzSamplings) cat("\n 1) given grouping") cat("\n test: ", x.names[test]," , ref: ", x.names[ref]) cat("\n min ssqg given grouping: ", ssq.given.situation.min, " , max (permutation) ssqg: ", ssq.given.situation.max) cat("\n min sampling ssqg: ", sample.min.value, " , max sampling ssqg: ", sample.max.value, " - row with min: ", sample.min.row) cat("\n smaller solutions (", number.better, " / ", anzSamplings, ") , sample p value ", sample.p.value) cat("\n 2) min ssqg of all combinations: ", comb.min.value, " , max ssqg: ", comb.max.value, " - row with min: ", comb.min.row) cat("\n 3) selected min combination") cat("\n test: ", x.names[test.2]," , ref: ", x.names[ref.2]) cat("\n min ssqg selected combination: ", ssq2.given.situation.min, " , max (permutation) ssqg: ", ssq2.given.situation.max) cat("\n min sampling ssqg: ", sample2.min.value, " , max sampling ssqg: ", sample2.max.value, " - row with min: ", sample2.min.row) cat("\n smaller solutions (", number2.better, " / ", anzSamplings, ") , sample p value: ", sample2.p.value) cat("\n\n") # return results - given grouping + best in data grouping return(list(given.f90=x.hy1.f90, bestInData.f90=x.hy1a.f90)) } #horst.test1 <- perm.testing.1(x=horstLCA.2,ref=c(1,2,8,9),test=c(3,4,5,6,7,10,11),dataname="horstLCA.2")