# TODO: Add comment # # Author: E.Korsching 01-2024 ############################################################################### #### new standardization of raw data - proof of concept # D : difference # 1) prove if the samples are homogenous enough, by D per group and D between groups, def. threshold # 2) optimize step by step for each factor the group variance so that <-> normalize (still pending) # var per factor is minimized (and whole D is optimized) # 3) raining noise on data to see when diff is away / D gets similar # # 4) compare classical with no dimensional reduction on several levels : pure expression and graphs/cluster # test - aad/mad rsp. MD and Wilcoxon signed-rank test a1 <- c(3,4,3,2) a2 <- c(3,3,4,5) wilcox.test(a1, a2, paired=F) # amdbi e.amdbi(a1, a2, p=2) # MD formula e.aadiff(a1, p=2) # library(diffeR) # Wilcoxon signed rank test, dependent, paired wilcox.test( x=a1, y=a2, # or not NULL alternative = "two.sided", # "less" "greater" mu=0, # null h. - distribution of x or of x - y (paired) is symmetric about mu paired=T, # T if x+y , F if only x exact=NULL, correct=T, #digits.rank=7,# if float conf.int=F, conf.level=0.95) #Wilcoxon signed rank test with continuity correction #data: a1 and a2 #V = 1.5, p-value = 0.5862 #alternative hypothesis: true location shift is not equal to 0 #Obs: Habt ihr gleiche Ränge, sog. „ties“, wird die wilcox.test()-Funktion keinen exakten p-Wert ermitteln. #Hierfür verwendet ihr wilcox_test() aus dem coin-Paket. #Die Syntax ist ähnlich, bis auf distribution = “exact“ anstatt exact = TRUE # library(coin) # wilcoxon_test() a.df <- data.frame(aa=c(a1,a2), g=factor( rep(c("a", "b"), c(4,4)) ) ) a.wt <- wilcox_test(aa ~ g, data=a.df, distribution="exact") a.wt #Exact Wilcoxon-Mann-Whitney Test #data: aa by g (a, b) #Z = -1.0838, p-value = 0.4857 #alternative hypothesis: true mu is not equal to 0 # Mann-Whitney-U, independent wilcox.test( x=a1, y=a2, alternative = "two.sided", # "less" "greater" mu=0, # null h. - distributions of x and y differ by a shift of mu - alternative - they differ by some other shift - "greater" x is shifted to the right of y paired=F, # F exact=NULL, correct=T, #digits.rank=7,# if float conf.int=F, conf.level=0.95) #Wilcoxon rank sum test with continuity correction #data: a1 and a2 #W = 4.5, p-value = 0.3529 #alternative hypothesis: true location shift is not equal to 0 # Mann-Whitney-U, independent wilcox.test( x=c(3,7,8,6,2,19), # c(1) c(1,1,1,1,1,1) round(runif(6,1,20),0) c(3,7,8,6,2,19) y=c(15,6,11,9,18,12), # c(2) c(2,2,2,2,2,2) round(runif(6,1,20),0) c(15,6,11,9,18,12) alternative = "two.sided", # "less" "greater" mu=0, # null h. - distributions of x and y differ by a shift of mu - alternative - they differ by some other shift - "greater" x is shifted to the right of y paired=F, # F exact=NULL, correct=T, #digits.rank=7,# if float conf.int=F, conf.level=0.95) # effect size a.wt <- wilcox.test( x=a1, y=a2, # or not NULL alternative = "two.sided", # "less" "greater" mu=0, # null h. - distribution of x or of x - y (paired) is symmetric about mu paired=T, # T if x+y , F if only x exact=NULL, correct=T) # Calculate the standardised z statistic Z and call it Zstat a.Zstat <- qnorm(a.wt$p.value/2) # -0.5443311 # Calculate the effect size using abs(Zstat)/sqrt(number of pairs) abs(a.Zstat)/sqrt(4) # 0.2721655 [0.1 small 0.3 moderate >0.5 large , Cohen’s classification] # Rank test wrong here: independent type: position dependency not considered #### better test for our purpose library(coin) a.y1 <- c(1.83, 0.50, 1.62, 2.48, 1.68, 1.88, 1.55, 3.06, 1.30) a.y2 <- c(0.878, 0.647, 0.598, 2.05, 1.06, 1.29, 1.06, 3.14, 1.29) a.df <- data.frame( y = c(a.y1, a.y2), x = gl(2, length(a.y1)), block = factor(rep(seq_along(a.y1), 2)) ) a.out <- symmetry_test(y ~ x | block, data=a.df, distribution="exact", alternative="two.sided") attributes(a.out) pvalue(a.out) statistic(a.out) a.y1 <- c(3,7,8,6,2,19) # a1 c(1,1,1,1,1,1) c(3,7,8,6,2,19) mean(c(3,7,8,6,2,19)) 7.5 a.y2 <- c(15,6,11,9,18,12) # a2 c(2,2,2,2,2,2) c(15,6,11,9,18,12) mean(c(15,6,11,9,18,12)) 11.8 a.df <- data.frame( y = c(a.y1, a.y2), x = gl(2, length(a.y1)), block = factor(rep(seq_along(a.y1), 2)) ) symmetry_test(y ~ x | block, data=a.df, distribution="exact", alternative="two.sided") ### alternatives - regression models # About the data, data$y contains count data and data$site_name is a factor with 9 levels. # There are 54 values in data$y, with 6 values per level of data$site_name. lrtest(glm(data$y ~ 1), glm(data$y ~ data$site_name, family="poisson")) # Likelihood ratio test #Model 1: data$y ~ 1 #Model 2: data$y ~ data$site_name ##Df LogLik Df Chisq Pr(>Chisq) #1 2 -89.808 #2 9 -31.625 7 116.37 < 2.2e-16 *** # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 lrtest(glm(data$y ~ 1, family="poisson"), glm(data$y ~ data$site_name, family="poisson")) # Likelihood ratio test #Model 1: data$y ~ 1 #Model 2: data$y ~ data$site_name ##Df LogLik Df Chisq Pr(>Chisq) #1 1 -54.959 #2 9 -31.625 8 46.667 1.774e-07 *** # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 waldtest(glm(data$y ~ data$site_name, family="poisson")) # Wald test #Model 1: data$y ~ data$site_name #Model 2: data$y ~ 1 #Res.Df Df F Pr(>F) #1 45 #2 53 -8 0.7398 0.6562 waldtest(glm(data$y ~ 1, family="poisson"), glm(data$y~data$site_name, family="poisson")) # Wald test #Model 1: data$y ~ 1 #Model 2: data$y ~ data$site_name #Res.Df Df F Pr(>F) #1 53 #2 45 8 0.7398 0.6562 #### permutation approaches # Vanhove blog - inspired by - https://janhove.github.io/posts/2015-02-26-explaining-key-concepts-using-permutation-tests/ # install.packages("data.table") see also (Steinhaus)-Johnson-Trotter-algorithm #library(data.table) #CJ() ## for variances # Define a function to compute difference in group variances a.var.diff <- function(data, Group1) { diff.var <- var(data[Group1]) - var(data[- Group1]) return(diff.var) } a.diffs <- mapply(a.var.diff, Group1=a.combinat, MoreArgs=list(data=actual.data) ) dotchart(sort(a.diffs), pch=16, xlab="Group variance difference", main="Group variance differences\nfor all possible combinations") abline(v=c(a.var.diff(actual.data, a1), -a.var.diff(actual.data, a1)), col="red", lty=2) mtext("ordered permutation results", side=2, line=1) ## for means - our case # read in actual data actual.data <- c( c(2, 7, 6), c(1, 1, 5) ) # define group a1 <- c(1, 2, 3) # define a function that computes the difference in means # between one part of a vector (indices in Group1) and the remaining part (indices NOT in Group1) a.mean.diff <- function(data, group) { diff.mean <- mean(data[group]) - mean(data[- group]) return(diff.mean) } a.mean.diff2 <- function(data, group) { diff.mean <- abs( mean(data[group]) - mean(data[- group]) ) return(diff.mean) } a.aadbi1 <- function(data, group, p) { a <- e.aadbi(x=data[group], y=data[- group], p=p) return(a) } a21 <- a.mean.diff(actual.data, a1) # [1] 2.666667 a22 <- a.mean.diff2(actual.data, a1) # [1] 2.666667 a23 <- a.aadbi1(actual.data, a1, p=1) # sample # [1] 3.75 a24 <- a.aadbi1(actual.data, a1, p=2) # population # [1] 3.333333 # combinations a.combinat <- combn(1:6, 3, simplify=F) # a <- combn(1:6, 2, simplify=F) # remove duplicates - ordered : therfore cut at half a.nc2 <- length(a.combinat) / 2 a.combinat2 <- a.combinat[1:a.nc2] pdf("results2025/perm_3_3_md_aadbi.pdf", width=7, height=9) par(mfrow=c(2,2)) # mean diff # apply function mean.diff - for every combination of indices in combinat - to actual.data a.diffs <- mapply(a.mean.diff, group=a.combinat, MoreArgs=list(data=actual.data) ); cat("\n diffs:",round(a.diffs,4),"\n") dotchart(sort(a.diffs), pch=16, xlab=paste("1) comb - md, red:+/-2.6,\nmin ",round(min(a.diffs), 4),", mean ",round(sum(a.diffs) / length(a.diffs), 4),", max ",round(max(a.diffs), 4),sep=""), main="Group mean differences\nfor all possible combinations") abline(v=c(a21, -a21), col="red", lty=2) mtext("ordered permutation results", side=2, line=1) a.diffs <- mapply(a.mean.diff, group=a.combinat2, MoreArgs=list(data=actual.data) ); cat("\n diffs:",round(a.diffs,4),"\n") dotchart(sort(a.diffs), pch=16, xlab=paste("2) comb1/2 - md, red:+/-2.6,\nmin ",round(min(a.diffs), 4),", mean ",round(sum(a.diffs) / length(a.diffs), 4),", max ",round(max(a.diffs), 4),sep=""), main="Group mean differences\nfor all unique combinations") abline(v=c(a21, -a21), col="red", lty=2) mtext("ordered permutation results", side=2, line=1) a.diffs <- mapply(a.mean.diff2, group=a.combinat, MoreArgs=list(data=actual.data) ); cat("\n diffs:",round(a.diffs,4),"\n") dotchart(sort(a.diffs), pch=16, xlab=paste("3) comb abs - md, red:2.6,\n mean ",round(sum(a.diffs) / length(a.diffs), 4),", max ",round(max(a.diffs), 4),sep=""), main="Group mean differences\nfor all possible combinations") abline(v=c(a21, -a21), col="red", lty=2) mtext("ordered permutation results", side=2, line=1) a.diffs <- mapply(a.mean.diff2, group=a.combinat2, MoreArgs=list(data=actual.data) ); cat("\n diffs:",round(a.diffs,4),"\n") dotchart(sort(a.diffs), pch=16, xlab=paste("4) comb1/2 abs - md, red:2.6,\n mean ",round(sum(a.diffs) / length(a.diffs), 4),", max ",round(max(a.diffs), 4),sep=""), main="Group mean differences\nfor all unique combinations") abline(v=c(a22, -a22), col="red", lty=2) mtext("ordered permutation results", side=2, line=1) # aadbi a.diffs <- mapply(a.aadbi2, group=a.combinat, MoreArgs=list(data=actual.data, p=2) ); cat("\n diffs:",round(a.diffs,4),"\n") dotchart(sort(a.diffs), pch=16, xlab=paste("5) comb abs - mad, red:3.75, blue:3.3,\n mean ",round(sum(a.diffs) / length(a.diffs), 4),", max ",round(max(a.diffs), 4),sep=""), main="Group (mad/aad) differences\nfor all possible combinations") abline(v=a23, col="red", lty=2) abline(v=a24, col="blue", lty=2) mtext("ordered permutation results", side=2, line=1) a.diffs <- mapply(a.aadbi2, group=a.combinat2, MoreArgs=list(data=actual.data, p=2) ); cat("\n diffs:",round(a.diffs,4),"\n") dotchart(sort(a.diffs), pch=16, xlab=paste("6) comb1/2 abs - mad, red:3.75, blue:3.3,\n mean ",round(sum(a.diffs) / length(a.diffs), 4),", max ",round(max(a.diffs), 4),sep=""), main="Group (mad/aad) differences\nfor all unique combinations") abline(v=a23, col="red", lty=2) abline(v=a24, col="blue", lty=2) mtext("ordered permutation results", side=2, line=1) par(mfrow=c(1,1)) dev.off() # # combinations a.combinat <- combn(1:6, 2, simplify=F) pdf("results2025/perm_2_4_md_aadbi.pdf", width=7, height=9) par(mfrow=c(2,2)) # mean diff # apply function mean.diff - for every combination of indices in combinat - to actual.data a.diffs <- mapply(a.mean.diff, group=a.combinat, MoreArgs=list(data=actual.data) ); cat("\n diffs:",round(a.diffs,4),"\n") dotchart(sort(a.diffs), pch=16, xlab=paste("1) comb - md,\nmin ",round(min(a.diffs), 4),", mean ",round(sum(a.diffs) / length(a.diffs), 4),", max ",round(max(a.diffs), 4),sep=""), main="Group mean differences\nfor all possible combinations") mtext("ordered permutation results", side=2, line=1) a.diffs <- mapply(a.mean.diff2, group=a.combinat, MoreArgs=list(data=actual.data) ); cat("\n diffs:",round(a.diffs,4),"\n") dotchart(sort(a.diffs), pch=16, xlab=paste("2) comb abs - md,\nmin ",round(min(a.diffs), 4),", mean ",round(sum(a.diffs) / length(a.diffs), 4),", max ",round(max(a.diffs), 4),sep=""), main="Group mean differences\nfor all possible combinations") mtext("ordered permutation results", side=2, line=1) # aadbi plot.new() a.diffs <- mapply(a.aadbi1, group=a.combinat, MoreArgs=list(data=actual.data, p=2) ); cat("\n diffs:",round(a.diffs,4),"\n") dotchart(sort(a.diffs), pch=16, xlab=paste("3) comb abs - mad,\nmin ",round(min(a.diffs), 4),", mean ",round(sum(a.diffs) / length(a.diffs), 4),", max ",round(max(a.diffs), 4),sep=""), main="Group (mad/aad) differences\nfor all possible combinations") mtext("ordered permutation results", side=2, line=1) par(mfrow=c(1,1)) dev.off() # # combinations a.combinat <- combn(1:6, 1, simplify=F) pdf("results2025/perm_1_5_md_aadbi.pdf", width=7, height=9) par(mfrow=c(2,2)) # mean diff # apply function mean.diff - for every combination of indices in combinat - to actual.data a.diffs <- mapply(a.mean.diff, group=a.combinat, MoreArgs=list(data=actual.data) ); cat("\n diffs:",round(a.diffs,4),"\n") dotchart(sort(a.diffs), pch=16, xlab=paste("1) comb - md,\nmin ",round(min(a.diffs), 4),", mean ",round(sum(a.diffs) / length(a.diffs), 4),", max ",round(max(a.diffs), 4),sep=""), main="Group mean differences\nfor all possible combinations") mtext("ordered permutation results", side=2, line=1) a.diffs <- mapply(a.mean.diff2, group=a.combinat, MoreArgs=list(data=actual.data) ); cat("\n diffs:",round(a.diffs,4),"\n") dotchart(sort(a.diffs), pch=16, xlab=paste("2) comb abs - md,\nmin ",round(min(a.diffs), 4),", mean ",round(sum(a.diffs) / length(a.diffs), 4),", max ",round(max(a.diffs), 4),sep=""), main="Group mean differences\nfor all possible combinations") mtext("ordered permutation results", side=2, line=1) # aadbi plot.new() a.diffs <- mapply(a.aadbi1, group=a.combinat, MoreArgs=list(data=actual.data, p=2) ); cat("\n diffs:",round(a.diffs,4),"\n") dotchart(sort(a.diffs), pch=16, xlab=paste("3) comb abs - mad,\nmin ",round(min(a.diffs), 4),", mean ",round(sum(a.diffs) / length(a.diffs), 4),", max ",round(max(a.diffs), 4),sep=""), main="Group (mad/aad) differences\nfor all possible combinations") mtext("ordered permutation results", side=2, line=1) par(mfrow=c(1,1)) dev.off() # Special case 1-1 # read in actual data actual.data <- c( 7, 5 ) # combinations a.combinat <- combn(1, 1, simplify=F) pdf("results2025/perm_1_1_md_aadbi_special_case.pdf", width=7, height=9) par(mfrow=c(2,2)) # mean diff # apply function mean.diff - for every combination of indices in combinat - to actual.data a.diffs <- mapply(a.mean.diff, group=a.combinat, MoreArgs=list(data=actual.data) ); cat("\n diffs:",round(a.diffs,4),"\n") dotchart(sort(a.diffs), pch=16, xlab=paste("1) comb - md,\nmin ",round(min(a.diffs), 4),", mean ",round(sum(a.diffs) / length(a.diffs), 4),", max ",round(max(a.diffs), 4),sep=""), main="Group mean differences\nfor all possible combinations") mtext("ordered permutation results", side=2, line=1) a.diffs <- mapply(a.mean.diff2, group=a.combinat, MoreArgs=list(data=actual.data) ); cat("\n diffs:",round(a.diffs,4),"\n") dotchart(sort(a.diffs), pch=16, xlab=paste("2) comb abs - md,\nmin ",round(min(a.diffs), 4),", mean ",round(sum(a.diffs) / length(a.diffs), 4),", max ",round(max(a.diffs), 4),sep=""), main="Group mean differences\nfor all possible combinations") mtext("ordered permutation results", side=2, line=1) # aadbi plot.new() a.diffs <- mapply(a.aadbi1, group=a.combinat, MoreArgs=list(data=actual.data, p=2) ); cat("\n diffs:",round(a.diffs,4),"\n") dotchart(sort(a.diffs), pch=16, xlab=paste("3) comb abs - mad,\nmin ",round(min(a.diffs), 4),", mean ",round(sum(a.diffs) / length(a.diffs), 4),", max ",round(max(a.diffs), 4),sep=""), main="Group (mad/aad) differences\nfor all possible combinations") mtext("ordered permutation results", side=2, line=1) par(mfrow=c(1,1)) dev.off() ## decision for all combinations and md as a perm scenario ## p values in the system # first, definition of a hypothesis # h0: it is all random # ha: it is non random -> presumably a systematic effect (not tested !) # definition of an alpha error -> smaller than alpha: ha # 5% ? # all exhaustively created combinatorial results define the population space # length(a.diffs) # the original result - if not populating a degenerate state - # is defining the p value: 1/n results * dg (degeneration grade) at that level # so the correction supports h0 and not ha # - dg likely if # the group distances are small (supportive) # each group variance is small to zero (ambiguous) #if(dg>0){ dg / length(a.diffs) }else{ 1 / length(a.diffs) } ## Two sample permutation test (exhaustive) over all group elements and their group mean values a <- perm.test.two.groups(x=c( c(2, 7, 6), c(1, 1, 5) ), length.g1=3, alternative="g", dg=F, file="test111.pdf", main="Group mean differences\nfor all possible combinations", sub="test", height=9) #### compare aadbi aadbi2 with classical group mean difference # a1 <- 500 a3 <- 6 a <- matrix(0,a1,2) dimnames(a)[[2]] <- c("a.mean.diff","aadbi") # aadbi2 identical to mean1-mean2 for(i in 1:a1){ a2 <- round(runif(a3, 1, 20),0) a[i,1] <- mean(a2[1:(a3/2)]) - mean(a2[(a3/2+1):a3]) a[i,2] <- e.aadbi(a2[1:(a3/2)], a2[(a3/2+1):a3], p=2) # e.aadbi2 } pdf("results2025/perm_cf_meanDiff_aadbi_aadbi2.pdf", width=11, height=5) # text aaa <- paste( "Comparison of meanDiff aadbi aadbi2\n\n", "mean diff : classical\n", "aadbi : all local differences between two groups - abs on diff (see graph)\n", "aadbi2 : all local differences between two groups - no abs on diff -> behaves identical to mean-mean\n", sep="") write2PDF(aaa, col="black", cex=0.7, ww=120, ll=25) par(mfrow=c(1,2)) a <- a[order(a[,1]),] plot(1:a1, a[,1], type="p", xlab="", ylab="", ylim=range(a), col="blue") points(1:a1, a[,2], pch=3, col="red") a <- a[order(a[,2]),] plot(1:a1, a[,1], type="p", xlab="", ylab="", ylim=range(a), col="blue") points(1:a1, a[,2], pch=3, col="red") par(mfrow=c(1,1)) dev.off() # a1 <- 500 a3 <- 6 a <- matrix(0,a1,2) dimnames(a)[[2]] <- c("a.mean.diff2","aadbi") for(i in 1:a1){ a2 <- round(runif(a3, 1, 20),0) a[i,1] <- abs(mean(a2[1:(a3/2)]) - mean(a2[(a3/2+1):a3])) a[i,2] <- e.aadbi(a2[1:(a3/2)], a2[(a3/2+1):a3], p=2) } pdf("results2025/perm_cf_meanDiff2_aadbi.pdf", width=11, height=5) # text aaa <- paste( "Comparison of abs mean diff - aadbi\n\n", "abs mean diff2 : abs(mean-mean)\n", "aadbi : all local differences between two groups - abs on diff (see graph)\n", sep="") write2PDF(aaa, col="black", cex=0.7, ww=120, ll=25) par(mfrow=c(1,2)) a <- a[order(a[,1]),] plot(1:a1, a[,1], type="p", xlab="", ylab="", ylim=range(a), col="blue") points(1:a1, a[,2], pch=3, col="red") a <- a[order(a[,2]),] plot(1:a1, a[,1], type="p", xlab="", ylab="", ylim=range(a), col="blue") points(1:a1, a[,2], pch=3, col="red") par(mfrow=c(1,1)) dev.off() save.image()