# pairwiseSERE
# compute pairwise SERE statistics

# input : counts
# output : matrix of SERE values

# created october 19th, 2012
# Marie-Agnes Dillies


pairwiseSERE <- function( counts ){
  
  sere <- matrix( NA, ncol=ncol(counts), nrow=ncol(counts) )
  for (i in 1:ncol(counts)){
    for (j in 1:ncol(counts)){
      sere[i,j] <- sigfun_Pearson( counts[,c(i,j)] )
    }
  }
  colnames(sere) <- rownames(sere) <- colnames(counts)
  return( formatC(sere, format="f", digits=2) )
}

sigfun_Pearson <- function(observed) {
  #calculate lambda and expected values
  laneTotals<- colSums(observed);
  total <- sum(laneTotals)
  fullObserved <- observed[rowSums(observed)>0,];
  fullLambda <- rowSums(fullObserved)/total;
  fullLhat <- fullLambda > 0;
  fullExpected<- outer(fullLambda, laneTotals);

  #keep values
  fullKeep <- which(fullExpected > 0);
  
  #calculate degrees of freedom (nrow*(ncol -1) >> number of parameters - calculated (just lamda is calculated >> thats why minus 1)
  #calculate pearson and deviance for all values
  oeFull <- (fullObserved[fullKeep] - fullExpected[fullKeep])^2/ fullExpected[fullKeep] # pearson chisq test
  dfFull <- length(fullKeep) - sum(fullLhat!=0);
  
  return(c(sqrt(sum(oeFull)/dfFull)));
}