center2 <- function(a){ mean.ak. <- rowMeans(a) mean.a.l <- colMeans(a) mean.a <- mean(a) A <- sweep(a, MARGIN = 1, STATS = mean.ak., FUN = "-") A <- sweep(A, MARGIN = 2, STATS = mean.a.l, FUN = "-") A <- A + mean.a return(A) } hsic.kernel <- function(a,b,index){ a <- exp(-a*b^index) return(a) } csorgo.stat <- function(a,p,n){ s1 <- 1; s2 <- 1; s3 <- 1 for(j in 1:p){ s1 <- s1*a[[j]] s2 <- s2*rowSums(a[[j]]) s3 <- s3*sum(a[[j]]) } Jn <- sum(s1)/n^2 - 2*sum(s2)/n^(p+1) + s3/n^(2*p) return(Jn) } csorgo.test <- function(X,vecd,index=1,cte=rep(1,length(vecd)),N=1000){ X <- as.matrix(X) n <- nrow(X) p <- length(vecd) In <- diag(n) J <- col(In)= csorgo.0) )/(N+1) list(csorgo.statistic = csorgo.0, pvalue = pvalue) } csorgoSerial.test <- function(X,p,index=1,cte=1,N=1000){ X <- as.matrix(X) m <- nrow(X) Y <- X d <- ncol(Y) vecd <- rep(d,p) n <- m-p+1 Im <- diag(m) J <- col(Im)= csorgo.0))/(N+1) list(csorgo.statistic = csorgo.0, pvalue = pvalue) } WnB.stat <- function(A.p,taille,ord,p,n,RES) { # Initialisation vect.WnB <- rep(0,taille) pvalue.WnB <- rep(0,taille) # Loop for W_nB outputs in the order of the cardinality of subsets B. nb <- 0 # Counter for subsets B, |B|>1. for(cardB in 2:ord){ for(j in 1:(choose(p,cardB))){ nb <- nb+1 temp <- 1 for(ensB in RES[[cardB]][,j]) temp <- temp*A.p[[ensB]] num.WnB <- sum(temp) #computation of W_nB vect.WnB[nb] <- num.WnB/n^2 } } return(vect.WnB) } dCov.test <- function(X,vecd,ord=length(vecd),index=1,method="dcov",cte=rep(1,length(vecd)),alpha=0.05, N=1000, dependogram=TRUE){ if(method!="hsic") method="dcov" X <- as.matrix(X) n <- nrow(X) p <- length(vecd) In <- diag(n) J <- col(In)= WnB.perm[,1]))/(N+1) if(taille>1){ pvalue.perm <- matrix(0,nrow=taille,ncol=(N+1)) for(k in 1:(N+1)) pvalue.perm[,k] <- (1/2+rowSums(WnB.perm[,-1] >= WnB.perm[,k]))/(N+1) Fisher.stat <- -2*colSums(log(pvalue.perm)) Tippett.stat <- apply(pvalue.perm,2,min) Fisher.pvalue <- mean(Fisher.stat[-1] >= Fisher.stat[1]) Tippett.pvalue <- mean(Tippett.stat[-1] <= Tippett.stat[1]) pvalues <- pvalue.perm[,1] #Labels for subsets nb <- 0 for(cardB in 2:ord){ for(j in 1:(choose(p,cardB))){ nb <- nb+1 names(WnB.0)[nb] <- paste("{",paste(RES[[cardB]][,j], sep = ",",collapse = ","),"}") } } names(pvalues) <- names(WnB.0) if(dependogram){ critical.values <- apply(WnB.perm[,-1], 1 ,quantile, probs = beta) names(critical.values) <- names(WnB.0) par(mar=c(7, 4, 4, 2) + 0.1) matplot(WnB.0,type="h",ylim=c(0,max(c(max(critical.values),max(WnB.0)))),xlim=c(0,taille+1),main=ifelse(method=="hsic",expression(HSIC),expression(DCOV)),ylab="Statistic",xaxt="n") points((1:taille),critical.values,pch="-") axis(side=1,at=1:taille,labels=labels(WnB.0),las = 2) par(mar=c(5, 4, 4, 2) + 0.1) } } if(taille==1) list(method = method, statistic = WnB.0, pvalue = pvalue) else list(method = method, statistics = WnB.0, pvalues = pvalues, Fisher.statistic = Fisher.stat[1], Fisher.pvalue = Fisher.pvalue, Tippett.statistic = Tippett.stat[1], Tippett.pvalue = Tippett.pvalue) } WnBSerial.stat <- function(A.p,taille,ord,p,n,RES) { # Initialisation vect.WnB <- rep(0,taille) pvalue.WnB <- rep(0,taille) # Loop for W_nB outputs in the order of the cardinality of subsets B. nb <- 0 # Counter for subsets B, |B|>1. for(cardB in 2:ord){ for(j in 1:(choose(p-1,cardB-1))){ nb <- nb+1 temp <- 1 for(ensB in RES[[cardB]][,j]) temp <- temp*A.p[[ensB]] num.WnB <- sum(temp) #computation of W_nB vect.WnB[nb] <- num.WnB/n^2 } } return(vect.WnB) } dCovSerial.test <- function(X,p,ord=p,index=1,method="dcov",cte=1,alpha=0.05, N=1000, dependogram=TRUE){ if(method!="hsic") method="dcov" X <- as.matrix(X) m <- nrow(X) Y <- X d <- ncol(Y) vecd <- rep(d,p) n <- m-p+1 Im <- diag(m) J <- col(Im)= WnB.perm[,1]))/(N+1) if(taille>1){ pvalue.perm <- matrix(0,nrow=taille,ncol=(N+1)) for(k in 1:(N+1)){ pvalue.perm[,k] <- (1/2+rowSums(WnB.perm[,-1] >= WnB.perm[,k]))/(N+1) } Fisher.stat <- -2*colSums(log(pvalue.perm)) Tippett.stat <- apply(pvalue.perm,2,min) Fisher.pvalue <- mean(Fisher.stat[-1] >= Fisher.stat[1]) Tippett.pvalue <- mean(Tippett.stat[-1] <= Tippett.stat[1]) pvalues <- pvalue.perm[,1] #Labels for subsets nb <- 0 for(cardB in 2:ord){ for(j in 1:(choose(p-1,cardB-1))){ nb <- nb+1 names(WnB.0)[nb] <- paste("{",paste(RES[[cardB]][,j], sep = ",",collapse = ","),"}") } } names(pvalues) <- names(WnB.0) if(dependogram){ #critical values per same cardinality critical.values <- rep(0, ord - 1) begin <- 1 end <- 0 for (cardB in 2:ord) { end <- end + choose(p-1, cardB-1) vecB <- as.vector(WnB.perm[,-1][begin:end, ]) vecB <- sort(vecB) critical.values[cardB - 1] <- quantile(vecB, probs=beta) begin <- end + 1 } par(mar=c(7, 4, 4, 2) + 0.1) matplot(WnB.0,type="h",ylim=c(0,max(c(max(critical.values),max(WnB.0)))),xlim=c(0,taille+1),main=ifelse(method=="hsic",expression(HSIC),expression(DCOV)),ylab="Statistic",xaxt="n") begin <- 1 end <- 0 for (cardB in 2:ord) { end <- end + choose(p - 1, cardB - 1) segments(begin - 0.5, critical.values[cardB - 1], end + 0.5, critical.values[cardB - 1], lty = 4) begin <- end + 1 } axis(side=1,at=1:taille,labels=labels(WnB.0),las = 2) par(mar=c(5, 4, 4, 2) + 0.1) } } if(taille==1) list(method = method, statistics = WnB.0, pvalue = pvalue) else list(method = method, statistics = WnB.0, pvalues = pvalues, Fisher.statistic = Fisher.stat[1], Fisher.pvalue = Fisher.pvalue, Tippett.statistic = Tippett.stat[1], Tippett.pvalue = Tippett.pvalue) }