# TODO: Add comment # # Author: E.Korsching 11-2025 ############################################################################### #### standardization / normalization # sometimes overlapping with proximity approaches which contain normalizing functionality #### core e.var <- function(x, p=1){ # variance (more general: dispersion) xlen <- length(x) if(p==1){ n <- xlen-1 }else{ n <- xlen } # sample or population var v <- sum( (x-mean(x))^2 ) / n return(v) } e.sd <- function(x, p=1){ # sd standard deviation var <- e.var(x, p=p) sd <- sqrt( var ) return(sd) } e.aadiff <- function(x, p=2){ # univariat # mean absolute difference or absolut mean difference (maesure of dispersion) xlen <- length(x) if(p==1){ n <- xlen-1 }else{ n <- xlen } # sample or population z <- vector("numeric",(xlen^2-xlen)/2) k <- 1 for(i in 1:(xlen-1)){ for(j in (i+1):xlen){ z[k] <- abs(x[i]-x[j]) k <- k+1 } } aadiff <- sum(z) / (n^2) return(aadiff) } e.aadbi <- function(x, y, p=2){ # bivariate, no direction # arithmetic mean of absolute co-differences xlen <- length(x) ylen <- length(y) n <- xlen * ylen if(p==1){ n <- n-1 } # sample or population n z <- vector("numeric",n) k <- 1 for(i in 1:xlen){ for(j in 1:ylen){ z[k] <- abs(x[i]-y[j]) k <- k+1 } } aadbi <- sum(z) / n return(aadbi) } e.aadbi2 <- function(x, y, p=2){ # bivariate, direction => identical final outcome to classical: mean() - mean() # arithmetic mean of co-differences xlen <- length(x) ylen <- length(y) n <- xlen * ylen if(p==1){ n <- n-1 } # sample or population n z <- vector("numeric",n) k <- 1 for(i in 1:xlen){ for(j in 1:ylen){ z[k] <- x[i]-y[j] k <- k+1 } } aadbi2 <- sum(z) / n return(aadbi2) } # show details exp.sd <- function(x, p=1){ # test xlen <- length(x) if(p==1){ n <- xlen-1; n1 <- "(n-1)" }else{ n <- xlen; n1 <- "n" } # sample or population cat("\n used denominator ",n1) cat("\n x ",x) a <- mean(x) cat("\n mean(x) ",a) b <- x-a cat("\n x-mean(x) ",b) c <- (x-mean(x))^2 cat("\n (x-mean(x))^2 ",c) d <- sum( (x-mean(x))^2 ) cat("\n sum(x-mean(x))^2) ",d) e <- sum( (x-mean(x))^2 ) / n cat("\n",paste("sum(x-mean(x))^2)/",n1,sep=""),e) cat("\n sqrt()") f <- sqrt(e) cat("\n sd ",f,"\n") return() } exp.aad <- function(x, p=1){ # test xlen <- length(x) if(p==1){ n <- xlen-1; n1 <- "(n-1)" }else{ n <- xlen; n1 <- "n" } # sample or population var cat("\n used denominator ",n1) cat("\n x ",x) a <- mean(x) cat("\n mean(x) ",a) b <- x-a cat("\n x-mean(x) ",b) c <- abs(x-mean(x)) cat("\n abs(x-mean(x)) ",c) d <- sum( abs(x-mean(x)) ) cat("\n sum(abs(x-mean(x))) ",d) e <- sum( abs(x-mean(x)) ) / n cat("\n",paste("sum(abs(x-mean(x))/",n1,sep="")) cat("\n mad ",e,"\n") return() } #### norm intra column/sample (in general also inter - so by row - see below) # a1 <- c(1,3,2,5,3,5) # div by max norm.xi.max.x <- function(x){ # xi / max(x) maxx <- max(x) norm <- x/maxx # 0..1 return(norm) } # sum( norm.xi.max.x(a1) ) # rescale norm # see adaptScale and adaptScale.E in 0functions/0general # mean norm norm.xi.mean.x <- function(x){ # xi - mean(x) mx <- mean(x) norm <- x-mx # sum --> ~ 0 return(norm) } # sum( norm.xi.mean.x(a1) ) # mean range norm norm.xi.mr.x <- function(x){ # xi - mean(x) / max(x)-min(x) mx <- mean(x) norm <- (x-mx) / (max(x)-min(x)) # sum --> ~ 0 return(norm) } # sum( norm.xi.mr.x(a1) ) # unit sum norm.xi.unit.sum <- function(x){ # xi / sum(xi) sx <- sum(x) if(sx==0){ stop("vector sum is 0\n") } norm <- x/sx return(norm) # sum =1 } # sum( norm.xi.unit.sum(a1) ) # unit length / L2norm norm.xi.L2 <- function(x){ # L2norm - denominator is the Euclidean magnitude of the vector d <- sqrt( sum( x^2 ) ) if(d==0){ stop("magnitude/unit length is 0") } L2n <- x / d # sum --> ~ 2.223782 return(L2n) } # sum( norm.xi.L2(a1) ) # Z score sd norm.xi.Zsd <- function(x){ # xi - mean(x) / sd typical Z score n <- length(x) if(n<=1){ stop("n or n-1 is 0") } me <- mean(x) sd <- sqrt( sum((x-me)^2) / (n-1) ) if(sd==0){ stop("z score: SD is 0") } norm <- (x-me) / sd # sum --> ~ 0 [2.4e-16] return(norm) } # sum( norm.xi.Zsd(a1) ) # Z score aad (absolut average deviation / mean absolut deviation) norm.xi.Zaad <- function(x){ # xi - mean(x) / aad n <- length(x) if(n==0){ stop("n is 0") } me <- mean(x) aad <- sum( abs(x-me) ) / n if(aad==0){ stop("z: aad is 0") } norm <- (x-me) / aad # sum --> ~ 0 [7.7e-16] return(norm) } # sum( norm.xi.Zaad(a1) ) # Z score mad (median absolut deviation) norm.xi.Zmad <- function(x){ # xi - median(x) / mad mx <- median(x) mad <- median( abs(x-mx) ) if(mad==0){ stop("z: mad is 0") } norm <- (x-mx) / mad # sum --> ~ 0.6666667 return(norm) } # sum( norm.xi.Zmad(a1) ) # div by log xprod norm.xi.log.xprod <- function(x){ # ln(xi) / ln(psum(xi)) xlen <- length(x) if(sum(x<=1)>=1){ x[which(x<=1)] <- 1.000001; warning("numerical problem") } # acceptable for us lnx <- log(x, base=2) ps <- prod(x) # small increase if 1.000001,1.000001, ... lnps <- log(ps, base=2) norm <- lnx / lnps # sum --> 1 return(norm) } # 2^0.668 =1.588 ; log(2, 2) =1 ; log(1, 2) =0 ; log(0.5, 2) =-1 ; log(0.0000001, 2) =-23.25 # sum( norm.xi.log.xprod(a1) ) # normalize only between samples (inter column) scale.iqrM <- function( x, # matrix, data.frame plot.it=T, # plot normalization graph? t.name="", # name for title logscale="", # y axis with log scale: "y" save.graph=F, # save graph in pdf file file.name="", # save data in file (needs name) outlier=F, # show outlier mfrow=c(2,1), srt=90, cex=0.6, line=0, adj=1, col.f=c("blue","green"), col.b="black" ){ # normalization for columns # ini nr <- dim(x)[1] nc <- dim(x)[2] cnames <- names(x) mar1 <- c(1,5,4,1) mar2 <- c(8,5,1,1) if(plot.it){ # before scaling #layout layout( matrix(c(1,2), nrow=2, ncol=1, byrow=TRUE), widths=c(34), #relative proportions heights=c(9,13), #relative proportions respect=TRUE ) # layout.show(2) #function to check the layout if(save.graph){ # png(filename=paste(getwd(),"/results/",file.name,".png", sep=""), # width=800, height=800, units="px", pointsize=12, bg="white") pdf(file=paste(getwd(),"/",file.name,".",format(Sys.time(), "%Y%m%d%H%M"),".pdf", sep=""), width=11, height=7, onefile=T, title="normalization", pointsize=12) } cex.sav <- par()$cex adj.sav <- par()$adj par(cex=cex,adj=adj) par(mar=mar1) if(logscale=="y"){ xx <- x+1 }else{ xx <- x } boxplot(x=xx, # wisker : 1.5 * IQR range top/bottom outline=outlier, notch=F, ylim=range(boxplot(xx,plot=F)$stats), ylab=if(logscale=="y"){"intensity (+1)"}else{"intensity"}, cex.lab=1.1, boxwex= .7, log=logscale, plot=T, axes=F, col=col.f[1], border=col.b ) axis( side=2, cex.axis=1.1, cex.lab=1.1 ) tmp <- paste(t.name," - IQR median normalization\n") mtext(text=tmp, side=3, line=line, outer=F, adj=1) } # normalization new version 2023-01 # Adjust IQ ranges to be the same as max of IQRs colIQR <- apply(x, 2, IQR, na.rm=T) # use standard fn -> see IQR help divisor <- colIQR/max(colIQR) # for each col the appropriate constant for(i in 1:nc){ x[,i] <- x[,i]/divisor[i] #col ranges werden durch den Teiler angeglichen } # Adjust medians to be the same as max of medians colMed <- apply(x, 2, median, na.rm=T) adjustment <- max(colMed) - colMed for(i in 1:nc){ x[,i] <- x[,i]+adjustment[i] } xlabel <- paste("scaling: to interquartile range and median: X.adj=Xi * max(IQR(Xi))/IQR(Xi) , X.norm=X.adj * max(median(X.adj.i))/median(X.adj.i)") cat("\n Cols are normalized\n") cat("\n the new values are in the range: ", round(range(x, na.rm=T)[1],2), " - ", round(range(x, na.rm=T)[2],2), "\n") cat("\n NA values ? : ",sum(is.na(x))," [0:none,>0:yes] \n") if(plot.it){ # after scaling par(mar=mar2) if(logscale=="y"){ xx <- x+1 }else{ xx <- x } boxplot(x=xx, outline=outlier, notch=F, ylim=range(boxplot(xx,plot=F)$stats), ylab=if(logscale=="y"){"intensity (+1)"}else{"intensity"}, cex.lab=1.1, boxwex= .7, log=logscale, plot=T, axes=F, col=col.f[2], border=col.b ) if(is.null(cnames)){ axis( side=1 , at=c(1:nc) , labels=paste(1:nc,sep=","), tick=F, las=2, srt=srt ) }else{ axis( side=1 , at=c(1:nc) , labels=cnames(x), tick=F, las=2, srt=srt ) } axis( side=2 , cex.axis=1.1, cex.lab=1.1) if(save.graph){ dev.off() } } # restore par(mfrow=c(1,1), mar=c(5,4,4,2), cex=cex.sav, adj=adj.sav) # save file if(file.name!=""){ write.table(file.name, x, append=F, sep="\t", dec=".", row.names=T, col.names=NA) } return(x) } # a <- scale.iqrM( cbind(rnorm(20),rnorm(20),rnorm(20),rnorm(20),rnorm(20)), t.name="Test") # sigmoid / logistic function norm.sigmoid <- function(x){ 1 / (1 + exp(-x)) } # sigmoid / logistic function [0..1] base e #a <- rnorm(20); range(a); quantile(x=a, probs=0.999) #plot(a, rep(-0.1,length(a)), col="gray", xlim=c(-2,2), ylim=c(-0.2,1), xlab="y_in", ylab="y_out", pch=3) #points(a, norm.sigmoid(a), col="blue") # tanh (hyperbolic tangent) norm.tanhyp <- function(x){ (2 / (1 + exp(-2*x))) -1 } # tanh (hyperbolic tangent) [-1..1] base e #points(a, adaptScale(norm.tanhyp(a), minS=-1,maxS=1,minT=0,maxT=1), col="red") range adjusted [0..1] # adjusted tanh (hyperbolic tangent) norm.adj.tanhyp <- function(x){ (2 / (1 + exp(-2*x/quantile(x, probs=0.999)))) -1 } # adjusted tanh (hyperbolic tangent) [-1..1] base e #points(a, adaptScale(norm.adj.tanhyp(a), minS=-1,maxS=1,minT=0,maxT=1), col="pink") range adjusted [0..1] #### statistical tests # Two sample permutation test (exhaustive) over all group elements and their group mean values perm.test.two.groups <- function(x, length.g1, alternative="t", dg=F, weight=NULL, file=NULL, main="", sub="", height=9){ # conservative permutation test # on a significant difference (d) between two groups # x: vector with all measurements: first group first and then second group # length.g1: length of group 1 (integer) # alternative: two.sided "t", greater "g", less "l" based on initial reference # dg: T: consider degeneration grade (dg) for the original state as weakness - correct by multiplying p with dg # weight: value: apply a weight to the p values based on all possible pairwise single differences (d) and the # consistency of their sign - also considering d=0 cases - value: 0..1 e.g. 0.1 # file: not NULL: string: "path/name.extension" is saving a graph if length.g1 <= 4 tlen <- length(x) group <- 1:length.g1 ylen <- tlen - length.g1 mean.diff <- function(data, group){ md <- mean(data[group]) - mean(data[- group]) return(md) } # direct result dr <- mean.diff(data=x, group=group) dr1 <- mean(x[group]) dr2 <- mean(x[- group]) # combinations combinat <- combn(1:tlen, length.g1, simplify=F) # apply diffs <- mapply(mean.diff, group=combinat, MoreArgs=list(data=x) ) l.diffs <- length(diffs) gr.dr <- sum(diffs >= dr) le.dr <- sum(diffs <= -dr) tw.dr <- gr.dr + le.dr if(alternative=="t"){ a.dr <- tw.dr }else if(alternative=="g"){ a.dr <- gr.dr }else if(alternative=="l"){ a.dr <- le.dr } s.diffs <- sort(diffs) cat("\n diffs:",round(s.diffs,4),"\n") # dg if(dg){ n.dr <- sum(diffs == dr) + sum(diffs == -dr) p.dr <- (a.dr * n.dr) / l.diffs }else{ p.dr <- a.dr / l.diffs } # weight if(!is.null(weight)){ p.dr.adj <- perm.test.two.groups.weight(p.dr, weight, x, length.g1, ylen) }else{ p.dr.adj <- NA } # save graph if(!is.null(file) & length.g1<=4){ pdf(file=file, width=7, height=height) par(mfrow=c(1,1)) dotchart(s.diffs, pch=16, xlab=paste(sub,",\nmin ",round(min(diffs), 4),", mean ",round(sum(diffs) / length(diffs), 4),", max ",round(max(diffs), 4),sep=""), main=main) abline(v=c(dr, -dr), col=c("blue","gray"), lty=2) mtext("ordered permutation results", side=2, line=1) dev.off() } erg <- list( statistic=dr, p.value=p.dr, estimate=c(dr1,dr2), p.adj=p.dr.adj, method="Two sample permutation test (exhaustive) over all group elements and their group mean values", data.name=deparse(substitute(x)) ) attr(erg, "statistic") <- "permutation test" attr(erg, "estimate") <- c("mean of x","mean of y") return(erg) } #a <- perm.test.two.groups(x=c(8,3,7, 4,9,6), length.g1=3, alternative="l", dg=F, weight=0.1, # 0.1 c(2,7,6, 1,1,5) # file="test111.pdf", main="Group mean differences\nfor all possible combinations", sub="test", height=9) perm.test.two.groups.weight <- function(p.dr, weight, x, xlen, ylen){ # subfunction - conservative permutation test # weighting based on a sign score of the pairwise differences # p.dr: non weighted p value # x: basic data, xlen: length group one, ylen: length group two # weight: weight factor aadbi <- function(x, y, xlen, ylen, n){ # cf. => e.aadbi2, identical final outcome to classical: mean() - mean() # co-differences z <- vector("numeric",n) k <- 1 for(i in 1:xlen){ for(j in 1:ylen){ z[k] <- sign( x[i]-y[j] ) # 1,0,-1 k <- k+1 } } return(z) } n <- xlen * ylen z <- aadbi(x=x[1:xlen], y=x[(xlen+1):(xlen+ylen)], xlen, ylen, n) cat("\nz ",z," n ",n) zlen <- length(z) za <- sum(z == 1) # greater zb <- sum(z == 0) # d=0 zc <- sum(z == -1) # less if((za==zlen | zc==zlen) & zb==0){ #cat("\na",za,zc) w <- weight } else if(za=4 & za<(n-1) & zc<(n-1)){ #cat("\nb",za,zc) w <- -weight/2 }else{ #cat("\nb",0) w <- 0 } } else if(zc>0 & zc<=zlen){ # check data #cat("\nc",zc) w <- -weight } p.adj <- p.dr * (1 - w) # factor: >1 p value --> 1, <1 p value --> 0 cat("\n") return(p.adj) } #perm.test.two.groups.weight(0.2, weight=0.1, x=c(2,9,3, 4,3,6), xlen=3, ylen=3) #c #perm.test.two.groups.weight(0.2, weight=0.1, x=c(8,9,7, 4,3,6), xlen=3, ylen=3) #a #perm.test.two.groups.weight(0.2, weight=0.1, x=c(8,3,7, 4,9,6), xlen=3, ylen=3) #b #perm.test.two.groups.weight(0.2, weight=0.1, x=c(8,3, 4,9), xlen=2, ylen=2) #b else #perm.test.two.groups.weight(0.2, weight=0.1, x=c(8,8, 9,7), xlen=2, ylen=2) #b #### applications ## norm. functions # uniform distribution sample.unif.test <- function(a, b, c, fnIn="", samplN=10, main=""){ # runif() - according to sampling number, runif a,b,c parameters # fnIn : norm function (see above) if(fnIn==""){ stop("Give a function name") } coloR <- c("blue","darkgoldenrod","cyan","coral2","cadetblue1","chartreuse2","darkorchid1","deeppink2","gray0","gold4") fnOp <- get(fnIn) xlen <- a erg <- matrix(0, samplN, xlen) erg2 <- erg for(i in 1:samplN){ tmp <- runif(a, b, c) # runif(n, min=0, max=1) erg[i,] <- tmp erg2[i,] <- fnOp(x=tmp) } r.range <- t(apply(erg, 1, range)) r.range2 <- t(apply(erg2, 1, range)) r.means <- apply(erg, 1, mean) r.means2 <- apply(erg2, 1, mean) r.sums <- apply(erg, 1, sum) r.sums2 <- apply(erg2, 1, sum) sample.stat <- cbind(r.range, r.means, r.sums, r.range2, r.means2, r.sums2) colnames(sample.stat) <- c("c.min", "c.max", "c.means", "c.sums", "c.min2", "c.max2", "c.means2", "c.sums2") c.range <- t(apply(erg, 2, range)) c.range2 <- t(apply(erg2, 2, range)) c.means <- apply(erg, 2, mean) c.means2 <- apply(erg2, 2, mean) c.sums <- apply(erg, 2, sum) c.sums2 <- apply(erg2, 2, sum) cross.sample.stat <- cbind(c.range, c.means, c.sums, c.range2, c.means2, c.sums2) colnames(cross.sample.stat) <- c("c.min", "c.max", "c.means", "c.sums", "c.min2", "c.max2", "c.means2", "c.sums2") #print(erg) par(mfrow=c(2,1)) r.e <- range(erg) r.x <- 1:xlen plot(0,0, type="n", xlim=c(0, xlen+1), ylim=c(r.e[1], r.e[2]), xlab="erg", ylab="strength", main="raw-runif") for(i in 1:samplN){ lines(x=r.x, y=erg[i,], col=coloR[i]) } r.e <- range(erg2) r.x <- 1:xlen plot(0,0, type="n", xlim=c(0, xlen+1), ylim=c(r.e[1], r.e[2]), xlab="erg2", ylab="norm strength", main=main) for(i in 1:samplN){ lines(x=r.x, y=erg2[i,], col=coloR[i]) } #hist(jittD[2:samplN], breaks=15, xlab="range of mean(jitter)", main="") par(mfrow=c(1,1)) return(list(sample=erg, score=erg2, sample.stat=sample.stat, cross.sample.stat=cross.sample.stat)) } #pdf("results02pr/02norm_tests1.pdf", width=7,height=7) # a <- sample.unif.test(a=10, b=0, c=150, fnIn="norm.xi.max.x", samplN=10, main="div by max") # a <- sample.unif.test(a=10, b=0, c=150, fnIn="norm.xi.rescale", samplN=10, main="rescale norm") # a <- sample.unif.test(a=10, b=0, c=150, fnIn="norm.xi.mean.x", samplN=10, main="mean norm") # a <- sample.unif.test(a=10, b=0, c=150, fnIn="norm.xi.mr.x", samplN=10, main="mean range norm") # a <- sample.unif.test(a=10, b=0, c=150, fnIn="norm.xi.unit.sum", samplN=10, main="unit sum") # a <- sample.unif.test(a=10, b=0, c=150, fnIn="norm.xi.unit.len", samplN=10, main="unit length") # a <- sample.unif.test(a=10, b=0, c=150, fnIn="norm.xi.Zsd", samplN=10, main="Zscore") # a <- sample.unif.test(a=10, b=0, c=150, fnIn="norm.xi.log.xprod", samplN=10, main="div by log xprod") #dev.off() # normal distribution sample.norm.test <- function(a, b1, b2, c1, c2, fnIn="", samplN=10, main=""){ # rnorm() - according to sampling number, rnorm a,b,c parameters # fnIn : norm function (see above) if(fnIn==""){ stop("Give a function name") } coloR <- c("blue","darkgoldenrod","cyan","coral2","cadetblue1","chartreuse2","darkorchid1","deeppink2","gray0","gold4") fnOp <- get(fnIn) xlen <- a erg <- matrix(0, samplN, xlen) erg2 <- erg for(i in 1:samplN){ mean.te <- runif(a, b1, b2) # runif(n, min, max) sd.te <- runif(a, c1, c2) tmp <- rnorm(a, mean.te, sd.te) # rnorm(n, mean, sd) erg[i,] <- tmp erg2[i,] <- fnOp(x=tmp) } r.range <- t(apply(erg, 1, range)) r.range2 <- t(apply(erg2, 1, range)) r.means <- apply(erg, 1, mean) r.means2 <- apply(erg2, 1, mean) r.sums <- apply(erg, 1, sum) r.sums2 <- apply(erg2, 1, sum) sample.stat <- cbind(r.range, r.means, r.sums, r.range2, r.means2, r.sums2) colnames(sample.stat) <- c("c.min", "c.max", "c.means", "c.sums", "c.min2", "c.max2", "c.means2", "c.sums2") c.range <- t(apply(erg, 2, range)) c.range2 <- t(apply(erg2, 2, range)) c.means <- apply(erg, 2, mean) c.means2 <- apply(erg2, 2, mean) c.sums <- apply(erg, 2, sum) c.sums2 <- apply(erg2, 2, sum) cross.sample.stat <- cbind(c.range, c.means, c.sums, c.range2, c.means2, c.sums2) colnames(cross.sample.stat) <- c("c.min", "c.max", "c.means", "c.sums", "c.min2", "c.max2", "c.means2", "c.sums2") #print(erg) par(mfrow=c(2,1)) r.e <- range(erg) r.x <- 1:xlen plot(0,0, type="n", xlim=c(0, xlen+1), ylim=c(r.e[1], r.e[2]), xlab="erg", ylab="strength", main="raw-rnorm") for(i in 1:samplN){ lines(x=r.x, y=erg[i,], col=coloR[i]) } r.e <- range(erg2) r.x <- 1:xlen plot(0,0, type="n", xlim=c(0, xlen+1), ylim=c(r.e[1], r.e[2]), xlab="erg2", ylab="norm strength", main=main) for(i in 1:samplN){ lines(x=r.x, y=erg2[i,], col=coloR[i]) } #hist(jittD[2:samplN], breaks=15, xlab="range of mean(jitter)", main="") par(mfrow=c(1,1)) return(list(sample=erg, score=erg2, sample.stat=sample.stat, cross.sample.stat=cross.sample.stat)) } #pdf("results02pr/02norm_tests2.pdf", width=7,height=7) # a <- sample.norm.test(a=10, b1=0, b2=200, c1=1, c2=2, fnIn="norm.xi.max.x", samplN=10, main="div by max") # a <- sample.norm.test(a=10, b1=0, b2=200, c1=1, c2=2, fnIn="norm.xi.rescale", samplN=10, main="rescale norm") # a <- sample.norm.test(a=10, b1=0, b2=200, c1=1, c2=2, fnIn="norm.xi.mean.x", samplN=10, main="mean norm") # a <- sample.norm.test(a=10, b1=0, b2=200, c1=1, c2=2, fnIn="norm.xi.mr.x", samplN=10, main="mean range norm") # a <- sample.norm.test(a=10, b1=0, b2=200, c1=1, c2=2, fnIn="norm.xi.unit.sum", samplN=10, main="unit sum") # a <- sample.norm.test(a=10, b1=0, b2=200, c1=1, c2=2, fnIn="norm.xi.unit.len", samplN=10, main="unit length") # a <- sample.norm.test(a=10, b1=0, b2=200, c1=1, c2=2, fnIn="norm.xi.Zscore", samplN=10, main="Zscore") # a <- sample.norm.test(a=10, b1=0, b2=200, c1=1, c2=2, fnIn="norm.xi.log.xprod", samplN=10, main="div by log xprod") #dev.off() ## sd mad Mad - test approach sample.norm.test.cf <- function(dataIn, fnIn="", p=c(1,1), rd=0, ab.l=NULL, pathPrefix="", main=""){ # sd mad Mad - test approach - take artificial data sets or real data # dataIn: df, matrix - working on rows with booth functions # fnIn : 2 functions to compare: e.g. c("e.sd","e.mad"), main: Text e.g. for test identification # p: p-parameter of the functions, rd:0: do nothing, 1: return uncommon data # ab.l: NULL: no adjustment: or: how close get both fn values by applying adjustment parameters, sqrt(2/pi) mad=ab.l*sd nr <- nrow(dataIn) nc <- ncol(dataIn) cat("\n dataIn nr ",nr, " nc ",nc,"\n") fnO1 <- get(fnIn[1]) fnO2 <- get(fnIn[2]) erg2 <- matrix(0, nr, 2) dimnames(erg2)[[2]] <- fnIn # for each full row for(i in 1:nr){ erg2[i,] <- c(fnO1(x=dataIn[i,], p=p[1]), fnO2(x=dataIn[i,], p=p[2])) } r.1.gr.eq.2 <- t(apply(erg2, 1, function(x){ x[1]>=x[2] } )) # vec logi r.1.eq.2 <- t(apply(erg2, 1, function(x){ x[1]==x[2] } )) r.1.lt.2 <- t(apply(erg2, 1, function(x){ x[1]0){ hist.plot.color(x=ucdiff, bin.num=19, norm=F, xlab="bins", ylab="frequency", xlim=range(ucdiff), ylog=F, digits=2, col.f="gray", col.sh=NULL, col.b=NULL, cex=1, h.title=paste("uncommon diff: ",fnIn[1]," < ",fnIn[2],sep=""), add=F) } # a <- length(r.1.gr.eq.2) b <- sum(r.1.gr.eq.2) c <- round((b/a)*100, 1) a1 <- length(r.1.eq.2) b1 <- sum(r.1.eq.2) c1 <- round((b1/a1)*100, 1) a2 <- length(r.1.lt.2) b2 <- sum(r.1.lt.2) c2 <- round((b2/a2)*100, 1) outxl <- 11 outx <- data.frame(tx=rep("",outxl),pos=rep("",outxl),ind=rep("",outxl), stringsAsFactors=F) outx[1,] <- c(paste("sample# ",a,sep=" "), 1, 0) outx[2,] <- c(paste(fnIn[1],">=",fnIn[2],":",b,", ",c,"%",sep=" "), 3, 2) outx[3,] <- c(paste(fnIn[1],"==",fnIn[2],":",b1,", ",c1,"%",sep=" "), 5, 2) outx[4,] <- c(paste(fnIn[1],"<",fnIn[2],":",b2,", ",c2,"%",sep=" "), 7, 2) c3 <- round(range(r.2.div.1),2) outx[5,] <- c(paste("div :",fnIn[2],"/",fnIn[1],": mean:",round(mean(r.2.div.1),2),", range:",paste(c3,collapse="-"),sep=" "), 9, 0) d <- round(range(r.1.diff.2), 1) outx[6,] <- c(paste("diff:",fnIn[1],"-",fnIn[2],": mean:",round(mean(r.1.diff.2),2),", range:",paste(d,collapse="-"),sep=" "), 11, 0) d1 <- round(range(erg2[r.1.gr.eq.2,1]), 1) d2 <- round(range(erg2[r.1.gr.eq.2,2]), 1) e <- round(mean(erg2[r.1.gr.eq.2,1]), 1) f <- round(mean(erg2[r.1.gr.eq.2,2]), 1) outx[7,] <- c(paste(fnIn[1],": range ge:",paste(d1,collapse="-"),", ",fnIn[2],": range gr:",paste(d2,collapse="-"),sep=" "), 13, 2) outx[8,] <- c(paste(fnIn[1],": mean ge:",e,", ",fnIn[2],": mean gr:",f,sep=" "), 15, 2) g <- round(range(erg2[!r.1.gr.eq.2,1]), 1) h <- round(range(erg2[!r.1.gr.eq.2,2]), 1) g1 <- round(mean(erg2[!r.1.gr.eq.2,1]), 1) h1 <- round(mean(erg2[!r.1.gr.eq.2,2]), 1) outx[9,] <- c(paste(fnIn[1],": range lt:",paste(g,collapse="-"),", ",fnIn[2],": range lt:",paste(h,collapse="-"),sep=" "), 17, 2) outx[10,] <- c(paste(fnIn[1],": mean lt:",g1,", ",fnIn[2],": mean lt:",h1,sep=" "), 19, 2) outx[11,] <- c(main, 22, 0) write2PDF(x=outx, add=F, col="black", cex=1) dev.off() if(rd==1){ tmp <- cbind(erg2[!r.1.gr.eq.2, ],erg[!r.1.gr.eq.2, ]) # return uncommon results where 1 < 2 # order diff decreasing t.1.diff.2 <- t(apply(tmp, 1, function(x){ x[1]-x[2] } )) tmp <- round(tmp[order(t.1.diff.2, decreasing=F),], 1) par(mfrow=c(1,1)) return(tmp) } par(mfrow=c(1,1)) return() } #old a <- sample.norm.test.cf(a=10, b1=1, b2=200, c1=0.5, c2=3, fnIn=c("e.sd","e.mad"), p=c(1,2,0), samplN=100000, rd=1, ab.l=sqrt(2/pi), main="test") #new: #sample.norm.test.cf(dataIn=data.prd.01, fnIn=c("e.sd","e.aadev"), p=c(1,1), rd=0, ab.l=sqrt(2/pi), pathPrefix="results_fn_tests", main="test02.pdf") # z score -- with sd mad Mad sample.z.test.cf <- function(dataIn, fnIn="", p=c(1,1,1), rd=0, ab.l=NULL, pathPrefix="", main=""){ # z score -- with sd mad Mad - take created data sets or real data # dataIn: df, matrix - working on rows with booth functions # fnIn : 3 functions: e.g. c("e.sd","e.mad","e.aadev"), main: Text e.g. for test identification # p: p-parameter of the functions, rd:0: do nothing, 1: return uncommon data # ab.l: NULL: no adjustment: or: how close get both fn values by applying adjustment parameters, sqrt(2/pi) mad=ab.l*sd nr <- nrow(dataIn) nc <- ncol(dataIn) cat("\n dataIn nr ",nr, " nc ",nc,"\n") fnO1 <- get(fnIn[1]) fnO2 <- get(fnIn[2]) fnO3 <- get(fnIn[3]) erg1 <- matrix(0, nr, nc) erg2 <- matrix(0, nr, nc) erg3 <- matrix(0, nr, nc) erg21 <- matrix(0, nr, nc) erg31 <- matrix(0, nr, nc) erg32 <- matrix(0, nr, nc) for(i in 1:nr){ erg1[i,] <- (dataIn[i,] - mean(dataIn[i,]))/fnO1(x=dataIn[i,], p=p[1]) # z erg2[i,] <- (dataIn[i,] - mean(dataIn[i,]))/fnO2(x=dataIn[i,], p=p[2]) erg3[i,] <- (dataIn[i,] - median(dataIn[i,]))/fnO3(x=dataIn[i,], p=p[3]) # erg21[i,] <- erg2[i,]-erg1[i,] erg31[i,] <- erg3[i,]-erg1[i,] erg32[i,] <- erg3[i,]-erg2[i,] } mean21 <- mean(colMeans(erg21)) # whole mean mean31 <- mean(colMeans(erg31)) mean32 <- mean(colMeans(erg32)) # stat pdf(paste(pathPrefix,main,sep="/"), width=7, height=11) par(mfrow=c(3,3)) numDens <- 50 # 1 hist hist.plot.p.l.curve(x=erg1, bin.num=10, norm=T, xlab="scale", ylab="density", xlim=NULL, digits=2, col.f="green", cex=1, h.title=paste("z all: ",fnIn[1],sep=""), add=F) abline(v=0) # a1.1 <- (a1[,3]+a1[,2])/2 # a1 <- hist.plot.p.l.curve...verbose=T # a1.2 <- find.peaks2(x=a1[,1], thresh=0.05, span=0.25, lspan=0.05, noisey=T) # a1.3 <- a1.1[a1.2] # abline(v=a1.3, col="blue") # a1.2 <- find.peaks3(x=a1.1, y=a1[,1], w=3, span=0.05) # a1.3 <- a1.1[a1.2$i] # abline(v=a1.2$x, col="blue") # 2 hist hist.plot.p.l.curve(x=erg2, bin.num=10, norm=T, xlab="scale", ylab="density", xlim=NULL, digits=2, col.f="green", cex=1, h.title=paste("z all: ",fnIn[2],sep=""), add=F) abline(v=0) # 3 hist hist.plot.p.l.curve(x=erg3, bin.num=10, norm=T, xlab="scale", ylab="density", xlim=NULL, digits=2, col.f="green", cex=1, h.title=paste("z all: ",fnIn[3],sep=""), add=F) abline(v=0) ## # 21 hist hist.plot.p.l.curve(x=erg21, bin.num=10, norm=T, xlab="scale", ylab="density", xlim=NULL, digits=2, col.f="green", cex=1, h.title=paste("z-diff all: ",fnIn[2]," - ",fnIn[1],sep=""), add=F) abline(v=0) # 31 hist hist.plot.p.l.curve(x=erg31, bin.num=10, norm=T, xlab="scale", ylab="density", xlim=NULL, digits=2, col.f="green", cex=1, h.title=paste("z-diff all: ",fnIn[3]," - ",fnIn[1],sep=""), add=F) abline(v=0) # 32 hist hist.plot.p.l.curve(x=erg32, bin.num=10, norm=T, xlab="scale", ylab="density", xlim=NULL, digits=2, col.f="green", cex=1, h.title=paste("z-diff all: ",fnIn[3]," - ",fnIn[2],sep=""), add=F) abline(v=0) # legend outxl <- 4 outx <- data.frame(tx=rep("",outxl),pos=rep("",outxl),ind=rep("",outxl), stringsAsFactors=F) outx[1,] <- c(paste("global mean",fnIn[2],"-",fnIn[1],": ",round(mean21,5),sep=" "), 1, 0) outx[2,] <- c(paste("global mean",fnIn[3],"-",fnIn[1],": ",round(mean31,5),sep=" "), 3, 0) outx[3,] <- c(paste("global mean",fnIn[3],"-",fnIn[2],": ",round(mean32,5),sep=" "), 5, 0) outx[4,] <- c(paste(main,sep=" "), 8, 0) write2PDF(x=outx, add=F, col="black", cex=1) par(mfrow=c(1,1)) dev.off() if(rd){ t.row <- apply(erg21,2,which.max) erg.21.max.d <- erg21[t.row,] erg.21.e1 <- erg1[t.row,] erg.21.e2 <- erg2[t.row,] t.row <- apply(erg31,2,which.max) erg.31.max.d <- erg31[t.row,] erg.31.e1 <- erg1[t.row,] erg.31.e2 <- erg3[t.row,] t.row <- apply(erg32,2,which.max) erg.32.max.d <- erg32[t.row,] erg.32.e1 <- erg2[t.row,] erg.32.e2 <- erg3[t.row,] xlen <- length(erg.21.max.d) e1 <- rbind(erg.21.max.d, rep(0,xlen), erg.31.max.d, rep(0,xlen), erg.32.max.d) # warnings?? e2 <- rbind(erg.21.e1, rep(0,xlen), erg.31.e1, rep(0,xlen), erg.32.e1) e3 <- rbind(erg.21.e2, rep(0,xlen), erg.31.e2, rep(0,xlen), erg.32.e2) return(list(diff=e1,prim1=e2,prim2=e3)) } return() } #sample.z.test.cf(dataIn=data.prd.01, fnIn=c("e.sd","e.mad","e.aadev"), p=c(1,2,0), rd=0, ab.l=sqrt(2/pi), pathPrefix="results_fn_tests", main="z-sd-mad-e.aadev-rd1.pdf")