# TODO: Add comment # # Author: E.Korsching 1.12.2011 refset, testset ############################################################################### proximity.bootstrap.sampling <- function(x, refset, testset, anzSamplings=100, method="Pearson", output=""){ # Create a proximity matrix and if selected bootstrap / permutation results of the proximity table # proximity table: rows: test group factors, columns: reference group factors # x: matrix, data.frame : will be processed column wise # testset: factor columns which should form the order # refset: factor columns / targets which rule the order of the test columns # anzSamplings: number of samplings (shuffling/bootstrap), 0: ssq best result, >0: ssq best sampling result # method: L2 (Euclidean), Pearson, Spearman -- special cor=NA handling: see below Pearson # output: "": return only proximity matrix # "boot.raw","boot.prox","boot.raw.col","boot.prox.col", "perm.raw","perm.prox","perm.raw.col","perm.prox.col" # return stacked matrix of proximity tables with the original unsampled object first # ini if(!(method %in% c("L2","Pearson","Spearman"))){ stop("\n wrong similarity method") } if(anzSamplings<0){ stop("\n anzSamplings has to be 0, 1, 2, ... (integer)") } anzSamplings <- anzSamplings+1 # add the (unsampled) proximity table of the original data set - at least that will be produced nr <- length(testset) nc <- length(refset) nrx <- nrow(x) ncx <- ncol(x) if(nr<2 | nc<2 | nrx<2){ cat("\n check: testset >=2 , refset >=2 , data rows >=2"); break } var.list <- apply(x,2,var) var.list <- var.list[order(var.list)] # calculate all ordered (0 ascending) variances cat("\n var ", var.list, "\n") # output var values if(method=="Pearson" & sum(var.list==0)>=1){ cat("\n warning -- var=0 -- special handling for 'cor' activated, check results \n") } # give warning (especially for cor measure) ## ct sub-functions -- one ref versus the test group if(method=="L2"){ #ct <- function(x,y) apply(y,2,function(x,y) sqrt(sum((x-y)^2)) ,x) ct <- function(x,y){ # L2 norm l2r <- sqrt(sum(x^2)) ncy <- ncol(y) l2t <- vector("numeric",ncy) for(i in 1:ncy){ l2t[i] <- sqrt(sum(y[,i]^2)) } # normalization x <- x/l2r for(i in 1:ncy){ y[,i] <- y[,i]/l2t[i] } # Euclidean distance erg <- vector("numeric",ncy) for(i in 1:ncy){ erg[i] <- sqrt(sum((x-y[,i])^2)) } return(erg) } } if(method=="Pearson"){ ct <- function(x,y){ cor.fn <- function(x,y){ tmp <- cov(x,y) # means: test if var(x) or var(y) =0 if(tmp==0){ # if cov 0 : -> set cor by definition to a value nearly 0 (randomly chosen of +0.00000001 or -0.00000001 ) tmp <- sample(x=c(0.00000001,-0.00000001), size=1) }else{ tmp <- cor(x,y) #calculate regular Pearson cor coefficient } return(tmp) } apply(y,2,cor.fn,x) } } if(method=="Spearman"){ ct <- function(x,y){ cor.fn <- function(x,y){ # tmp <- cov(x,y) # means: test if var(x) or var(y) =0 ########## für Spearman sinnvoll? ######### # if(tmp==0){ # if cov 0 : -> set cor by definition to a value nearly 0 (randomly chosen of +0.00000001 or -0.00000001 ) # tmp <- sample(x=c(0.00000001,-0.00000001), size=1) # }else{ tmp <- cor(x,y, method="Spearman") # calculate regular Spearman rank cor coefficient # } return(tmp) } apply(y,2,cor.fn,x) } } ## work -- output always a matrix #def z <- matrix(0, nr*(anzSamplings), nc) name.x <- names(x) dimnames(z)[[2]] <- name.x[refset] dimnames(z)[[1]] <- rep(name.x[testset], times=(anzSamplings)) y <- matrix(0, nrx, ncx) # dim of raw data #-- apply all ref versus a test group and join all ref-test vectors of ct() to a matrix #-- ct(): apply all test vs. the selected ref, return vec of results for this ref vs. test group ## normal proximity table -- in the case of additional sampling approaches it is the first section in the container matrix, # all the following blocks in the container matrix are samplings j <- 1 z[((j*nr)-(nr-1)):((j+1)*nr-nr),] <- apply(x[ ,refset,drop=F], 2, ct, x[ ,testset,drop=F]) # create a full prox. table (col-ref, row-test) if(anzSamplings==1){ return(z) } ## sampling (bootstrap) -- in 4 flavours if(output=="boot.raw"){ # bootstrap sampling on the raw data (all values) for(j in 2:(anzSamplings)){ y <- matrix(sample(as.vector(unlist(x)),replace=T), nrx, ncx) # basic 0-hypothesis (with replacement (bootstrap)) z[((j*nr)-(nr-1)):((j+1)*nr-nr), ] <- apply(y[ ,refset,drop=F], 2, ct, y[ ,testset,drop=F]) } return(z) } if(output=="boot.prox"){ # bootstrap sampling on the proximity table (all values) for(j in 2:(anzSamplings)){ z[((j*nr)-(nr-1)):((j+1)*nr-nr), ] <- sample(z[((1*nr)-(nr-1)):((1+1)*nr-nr), ],replace=T) } return(z) } if(output=="boot.raw.col"){ # bootstrap sampling on the raw data (per factor column) for(j in 2:(anzSamplings)){ for(i in 1:ncx){ y[,i] <- sample(x[,i],replace=T) } z[((j*nr)-(nr-1)):((j+1)*nr-nr), ] <- apply(y[ ,refset,drop=F], 2, ct, y[ ,testset,drop=F]) } return(z) } if(output=="boot.prox.col"){ # bootstrap sampling on the proximity table (per reference column) for(j in 2:(anzSamplings)){ for(i in 1:nc){ z[((j*nr)-(nr-1)):((j+1)*nr-nr), i] <- sample(z[((1*nr)-(nr-1)):((1+1)*nr-nr), i],replace=T) } } return(z) } ## sampling (permutation) -- in 4 flavours if(output=="perm.raw"){ # permutation sampling on the raw data (all values) for(j in 2:(anzSamplings)){ y <- matrix(sample(as.vector(unlist(x)),replace=F), nrx, ncx) # basic 0-hypothesis (without replacement (permutation)) [x is data.frame ergo a list obj] z[((j*nr)-(nr-1)):((j+1)*nr-nr), ] <- apply(y[ ,refset,drop=F], 2, ct, y[ ,testset,drop=F]) } return(z) } if(output=="perm.prox"){ # permutation sampling on the proximity table (all values) for(j in 2:(anzSamplings)){ z[((j*nr)-(nr-1)):((j+1)*nr-nr), ] <- sample(z[((1*nr)-(nr-1)):((1+1)*nr-nr), ],replace=F) } return(z) } if(output=="perm.raw.col"){ # permutation sampling on the raw data (per factor column) for(j in 2:(anzSamplings)){ for(i in 1:ncx){ y[,i] <- sample(x[,i],replace=F) } z[((j*nr)-(nr-1)):((j+1)*nr-nr), ] <- apply(y[ ,refset,drop=F], 2, ct, y[ ,testset,drop=F]) } return(z) } if(output=="perm.prox.col"){ # permutation sampling on the proximity table (per reference column) for(j in 2:(anzSamplings)){ for(i in 1:nc){ z[((j*nr)-(nr-1)):((j+1)*nr-nr), i] <- sample(z[((1*nr)-(nr-1)):((1+1)*nr-nr), i],replace=F) } } return(z) } } #aa2 <- proximity.bootstrap.sampling(He589.c, refset=c(1,5,3,2,12,13), testset=c(4,6:11,14:16), anzSamplings=4, method="cor", output="")