R/A-dep-tests.R

In IndependenceTests: Non-Parametric Tests of Independence Between Random Vectors

A.dep.tests <- function(Xmat,choice=1,d=0,m=d,freqname="",type="text") {

  if (!(choice %in% 1:2)) stop("choice should be 1 or 2")
  
  Xd <- function(x,d) {
    n <- length(x)
    V1 <- NULL
    for (i in 1:d) V1 <- cbind(V1,x[i:(n-d+i)])
    V1
  }
  
  
  if (choice==1) {
    if (is.table(Xmat)) {
      d <- length(dim(Xmat))
      tmat <- sapply(dimnames(Xmat),length)
    } else {
      if (freqname=="") {
        d <- ncol(Xmat)
        tmat <- sapply(lapply(Xmat,MARGIN=2,FUN=unique),length)
        V1 <- Xmat
      } else {
        Xmat <- eval(parse(text=paste("xtabs(formula=",freqname,"~.,data=Xmat)")))
        d <- length(dim(Xmat))
        tmat <- sapply(dimnames(Xmat),length)
      }
    }
  }
  
  if (choice==2) {
    if (d==0) stop("You should enter a value for d")
    if (m>d) stop("You should enter a value for m which is <=d")
    V1 <- Xd(Xmat,d) # la fonction Xd est ci-dessus
    #Xmat <- table(as.data.frame(V1))
    tmat <- sapply(lapply(data.frame(V1),MARGIN=2,FUN=unique),length)
  }
  
                                        
  RES <- as.list(1:(m-1))
  TheTAs <- RES

  if (choice==1) {
  
    for (cardA in 2:m) {
      RES[[cardA-1]] <- as.matrix(combn(d,cardA))
      TheTAs[[cardA-1]] <- as.data.frame(t(as.matrix(RES[[cardA-1]][1,])))
      dimnames(RES[[cardA-1]]) <- list(NULL,apply(RES[[cardA-1]],2,f<-function(...){paste(...,collapse=",")}))
      names(TheTAs[[cardA-1]]) <- dimnames(RES[[cardA-1]])[[2]]
    }
  
    ThefAs <- TheTAs
    ThepvalAs <- TheTAs
  
    for (cardA in 2:m) {
      bs <- choose(d,cardA)
      for (j in 1:bs) {
        
        A <- RES[[cardA-1]][,j]
        if (is.table(Xmat)) X2A <- summary(margin.table(Xmat, A))$statistic else X2A <- summary(table(as.data.frame(V1[,A])))$statistic
        
        somme <- 0
        if (cardA>2) {
          for (cardB in 2:(cardA-1)) {
            somme <- somme + sum(TheTAs[[cardB-1]][list(NULL,apply(as.matrix(combn(A,cardB)),2,f<-function(...){paste(...,collapse=",")}))[[2]]])
          }
        }
        
        TA <- X2A-somme
        TheTAs[[cardA-1]][paste(as.character(A),collapse=",")] <- TA
        
        dA <- length(A)
        I <- rep(0,dA)
        
        for (k in 1:dA) I[k] <- tmat[A[k]]
        fA <- prod(as.vector(I-1))
        ThefAs[[cardA-1]][paste(as.character(A),collapse=",")] <- fA
        ThepvalAs[[cardA-1]][paste(as.character(A),collapse=",")] <- pchisq(TA,df=fA,lower.tail=FALSE)
        
      }
      
    }

  } else { # choice==2

    for (cardA in 2:m) {
      RES[[cardA-1]] <- rbind(rep(1,choose(d-1,cardA-1)),as.matrix(combn(2:d,cardA-1))) # tous les ensembles A avec 1
      TheTAs[[cardA-1]] <- as.data.frame(t(as.matrix(RES[[cardA-1]][1,])))
      dimnames(RES[[cardA-1]]) <- list(NULL,apply(RES[[cardA-1]],2,f<-function(...){paste(...,collapse=",")}))
      names(TheTAs[[cardA-1]]) <- dimnames(RES[[cardA-1]])[[2]]
    }
    
    ThefAs <- TheTAs
    ThepvalAs <- TheTAs
  
    for (cardA in 2:m) {
      bs <- choose(d-1,cardA-1)
      for (j in 1:bs) {
        
        A <- RES[[cardA-1]][,j]
        if (is.table(Xmat)) X2A <- summary(margin.table(Xmat, A))$statistic else X2A <- summary(table(as.data.frame(V1[,A])))$statistic
     
        somme <- 0
        X2Awo1 <- 0
        if (cardA>2) {
        if (is.table(Xmat)) X2Awo1 <- summary(margin.table(Xmat, A[-1]))$statistic else X2Awo1 <- summary(table(as.data.frame(V1[,A[-1]])))$statistic
          for (cardB in 2:(cardA-1)) {
            somme <- somme + sum(TheTAs[[cardB-1]][list(NULL,apply(rbind(rep(1,choose(cardA-1,cardB-1)),as.matrix(combn(A[-1],cardB-1))),2,f<-function(...){paste(...,collapse=",")}))[[2]]])
          }
        }

        TA <- X2A-X2Awo1-somme
        TheTAs[[cardA-1]][paste(as.character(A),collapse=",")] <- TA
        
        dA <- length(A)
        I <- rep(0,dA)
        
        for (k in 1:dA) I[k] <- tmat[A[k]]
        fA <- prod(as.vector(I-1))
        ThefAs[[cardA-1]][paste(as.character(A),collapse=",")] <- fA
        ThepvalAs[[cardA-1]][paste(as.character(A),collapse=",")] <- pchisq(TA,df=fA,lower.tail=FALSE)
        
      }
      
    }



  }
    

  names(TheTAs) <- paste("|A|=",2:m,sep="")
  names(ThefAs) <- names(TheTAs)
  names(ThepvalAs) <- names(TheTAs)
  
  res <- list(TA=TheTAs,fA=ThefAs,pvalA=ThepvalAs)
  
  
  X <- c()
  for (cardA in 2:m) {
    X <- rbind(X,eval(parse(text=paste("t(rbind(res$TA$`|A|=",cardA,"`,res$fA$`|A|=",cardA,"`,res$pvalA$`|A|=",cardA,"`))",sep=""))))
  }
  
  colnames(X) <- c("TA","fA","pvalA")
  
# Test of Mutual independence
  M <- sum(X[,1])
  fM <- sum(X[,2])
  pvalM <- pchisq(M,df=fM,lower.tail=FALSE)
  
  X2tmp <- rbind(c("","",""),round(X,2),c("","",""),round(c(M,fM,pvalM),2))
#  Anames <- sapply(strsplit(rownames(X),"",fixed=TRUE),FUN=function(x) paste(x,collapse=","))
  Anames <- rownames(X)
  if (choice==1) rownames(X2tmp) <- c("A",Anames,"","X2") else rownames(X2tmp) <- c("A",Anames,"","Y2")


  ans <- if (choice == 1) {
    c(res, list(X = X, X2 = M, f = fM, pval = pvalM))
  } else {
    c(res, list(X = X, Y2 = M, f = fM, pval = pvalM))
  }

  ans$print.table <- X2tmp
  ans$type <- type
  ans$choice <- choice
  ans$X.summary <- X

  class(ans) <- "A.dep.tests"
  ans
      
}



print.A.dep.tests <- function(x, ...) {
  if (identical(x$type, "html")) {
    sortie <- cbind(subset = rownames(x$X.summary), round(x$X.summary, 2))
    sortie <- rbind(sortie, c("", sortie[nrow(sortie), 2:4]))

    if (x$choice == 1) {
      sortie[nrow(sortie) - 1, ] <- c(
        "<B>all</B>",
        "<B>X<sup>2</sup></B>",
        "<B>total</B>",
        "<B>p.value</B>"
      )
    }

    if (x$choice == 2) {
      sortie[nrow(sortie) - 1, ] <- c(
        "<B>all</B>",
        "<B>Y<sup>2</sup></B>",
        "<B>total</B>",
        "<B>p.value</B>"
      )
    }

    print(xtable(data.frame(sortie)), type = "html", include.rownames = FALSE)
  } else {
    print(data.frame(x$print.table))
  }

  invisible(x)
}