R/test2graph.R

In kyoustat/PNAS: 'PNAS' is for Network Analysis and Statistical methods

#' test 2 graph in mm-way
#' 
#' @export
test2graph <- function(graph1, graph2, anchor=c("mean","median"), method=c("ad","ks","kw","ns","steel")){
  ##############################################################
  # preliminary : data names
  DNAME = paste(deparse(substitute(graph1)),"and",deparse(substitute(graph2)))
  
  ##############################################################
  # Check Input & Transform into 'igraph' object if matrix
  input1 = graph1
  if (!check_network(input1)){
    stop("* test2graph : an input graph should be either an 'igraph' object or adjacency matrix.")
  }
  if (is.matrix(input1)){
    if (isSymmetric(input1)){
      input1 = igraph::graph_from_adjacency_matrix(input1, mode="undirected")  
    } else {
      input1 = igraph::graph_from_adjacency_matrix(input1, mode="directed")  
    }
  }
  input2 = graph2
  if (!check_network(input2)){
    stop("* test2graph : an input graph should be either an 'igraph' object or adjacency matrix.")
  }
  if (is.matrix(input2)){
    if (isSymmetric(input2)){
      input2 = igraph::graph_from_adjacency_matrix(input2, mode="undirected")  
    } else {
      input2 = igraph::graph_from_adjacency_matrix(input2, mode="directed")  
    }
  }
  
  ##############################################################
  # check condition - (1) both should be connected, and (2) set the anchor selection
  if ((!igraph::is_connected(input1))||(!igraph::is_connected(input2))){
    stop("* test2graph : we need both input graphs to be connected.")
  }
  anchor = match.arg(anchor)
  
  ##############################################################
  # branching to use 'testKgraph'
  method = match.arg(method)
  if (all(method=="ks")){
    # computation : effective resistance whose sqrt is metric
    metric1 = sqrt(PNAS::effective(input1))
    metric2 = sqrt(PNAS::effective(input2))
    
    
    # computation : two-sample kolmogorov-smirnov test statistic (DAS)
    ids = return_id(metric1, metric2, anchor)
    id1 = ids[1]; tgt1 = as.vector(metric1[id1,-id1])
    id2 = ids[2]; tgt2 = as.vector(metric2[id2,-id2])
  
    # compute p-value using different methods from twosamples  
    ksout = suppressWarnings(stats::ks.test(tgt1, tgt2, alternative = "two.sided"))
    thestat = as.double(ksout$statistic)
    pvalue  = as.double(ksout$p.value)
    names(thestat) = "Tks"
  } else {
    listgraph = list()
    listgraph[[1]] = input1
    listgraph[[2]] = input2
    result = testKgraph(listgraph, method=method)
    thestat = result$statistic
    pvalue  = result$p.value
  }

  
  ##############################################################
  # REPORT
  hname   = "Test for Equality of Two Networks"
  Ha    = "two networks are not equal"
  res   = list(statistic=thestat, p.value=pvalue, alternative = Ha, method=hname, data.name = DNAME)
  class(res) = "htest"
  return(res)
}


#   -----------------------------------------------------------------------
#' @keywords internal
#' @noRd
return_id <- function(metric1, metric2, method){
  m = nrow(metric1)
  n = nrow(metric2)
  if (all(method=="random")){
    id1 = sample(1:m, 1)
    id2 = sample(1:n, 1)
  } else if (all(method=="mean")){
    id1 = which.min(rowSums(metric1^2))[1]
    id2 = which.min(rowSums(metric2^2))[1]
  } else if (all(method=="median")){
    id1 = which.min(rowSums(metric1))[1]
    id2 = which.min(rowSums(metric2))[1]
  }
  return(c(id1,id2))
}


# # ## test with some named graphs
# # ## scenario 1. star graph of sizes
# library(igraph)
# library(DAS)
# star10c1 = make_star(10)
# star20c1 = make_star(20)
# star20c2 = make_star(20, center=2)
# 
# 
# # test2graph(star20c1, star20c1)
# metric1 = sqrt(effective(star20c1)); embed1 = DAS::cmds(metric1)$embed
# metric2 = sqrt(effective(star20c2)); embed2 = DAS::cmds(metric2)$embed
# 
# id1 = which.min(rowSums(metric1^2))[1]; tgt1 = metric1[id1,-id1]
# id2 = which.min(rowSums(metric2^2))[1]; tgt2 = metric2[id2,-id2]
# 
# par(mfrow=c(3,2))
# plot(star20c1, main="star20c1")
# plot(star20c2, main="star20c2")
# plot(embed1[,1], embed1[,2], main="mds of star20c1")
# plot(embed2[,1], embed2[,2], main="mds of star20c2")
# plot(ecdf(tgt1), main="ecdf::star20c1")
# plot(ecdf(tgt2), main="ecdf::star20c2")
# test2graph(star20c1, star20c2, anchor = "mean", method="ad")
# test2graph(star20c1, star20c2, anchor = "mean", method="cvm")
# test2graph(star20c1, star20c2, anchor = "mean", method="kuiper")
# test2graph(star20c1, star20c2, anchor = "mean", method="wass")
# test2graph(star20c1, star20c2, anchor = "mean", method="ks")
# ks.test(tgt1, tgt2)
# 
# 
# niter = 12345
# pvals = array(0,c(2,niter))
# for (i in 1:niter){
#   star20c1 = make_star(20)
#   star20c2 = make_star(20, center=2)
#   pvals[1,i] = test2graph(star20c1, star20c2, method="ad")$p.value
#   pvals[2,i] = test2graph(star20c1, star20c2, method="ks")$p.value
#   print(paste("iteration ",i,"/",niter," complete...", sep=""))
# }
# 
# allnames = c("ad","cvm","kuiper","wass","ks")
# x11()
# par(mfrow=c(1,2), pty="s")
# for (i in 1:2){
#   tgtp = as.vector(pvals[i,])
#   alphas = seq(from=0.001, to=0.999, length.out=100)
#   errors = rep(0,100)
#   for (j in 1:100){
#     tgtalpha  = alphas[j]
#     errors[j] = sum(tgtp <= tgtalpha)/length(tgtp)
#   }
#   pm = paste("error by ",allnames[i],sep="")
#   plot(alphas, errors, main=pm)
# }