R/testKgraph.R

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

#' test K graph in mm-way
#' 
#' @export
testKgraph <- function(glist, anchor=c("mean","median"), method=c("ad","kw","ns","steel")){
  ##############################################################
  # preliminary : data names and glist
  DNAME = deparse(substitute(glist))
  if ((!is.list(glist))||(length(glist)<2)){
    stop("* testKgraph : input 'glist' should be a list of length at least 2.")
  }
  anchor = match.arg(anchor)
  method = match.arg(method)
  N = length(glist)
  
  ##############################################################
  # extract ecdf for each graph
  listecdf <- list()
  for (n in 1:N){
    # 1. convert it to the graph
    now.graph = testgraph_trf(glist[[n]])
    # 2. compute metric space & minimal index
    now.metric = sqrt(PNAS::effective(now.graph))
    now.id = return_idsingle(now.metric, anchor)
    # 3. extract the distance vector
    listecdf[[n]] = as.vector(now.metric[now.id,-now.id])
  }
  
  ##############################################################
  # computation
  if (all(method=="ad")){
    adout = kSamples::ad.test(listecdf, method="asymptotic")$ad
    thestat = as.double(adout[1,1])
    pvalue  = as.double(adout[1,3])
    names(thestat) = "Tad"
  } else if (all(method=="kw")){
    kwout   = kSamples::qn.test(listecdf, method="asymptotic", test="KW")$qn
    thestat = as.double(kwout[1])
    pvalue  = as.double(kwout[2])
    names(thestat) = "Tkw"
  } else if (all(method=="vdW")){
    vdwout = kSamples::qn.test(listecdf, method="asymptotic", test="vdW")$qn
    thestat = as.double(vdwout[1])
    pvalue  = as.double(vdwout[2])
    names(thestat) = "TvdW"
  } else if (all(method=="ns")){
    nsout = kSamples::qn.test(listecdf, method="asymptotic", test="NS")$qn
    thestat = as.double(nsout[1])
    pvalue  = as.double(nsout[2])
    names(thestat) = "Tns"
  } else if (all(method=="steel")){
    stout = kSamples::Steel.test(listecdf, method="asymptotic")$st
    thestat = as.double(stout[1])
    pvalue  = as.double(stout[2])
    names(thestat) = "Tsteel"
  }
  
  ##############################################################
  # REPORT
  hname   = "Test for Homogeneity of K Networks"
  Ha    = "not all K networks are equal"
  res   = list(statistic=thestat, p.value=pvalue, alternative = Ha, method=hname, data.name = DNAME)
  class(res) = "htest"
  return(res)
}


#   -----------------------------------------------------------------------
#' @keywords internal
#' @noRd
testgraph_trf <- function(graph){
  input2 = graph
  if (!check_network(input2)){
    stop("* testKgraph : 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")  
    }
  }  
  if (!igraph::is_connected(input2)){
    stop("* testKgraph : we need every input graph to be connected.")
  }
  return(input2)
}
#' @keywords internal
#' @noRd
return_idsingle <- function(metric, method){
  if (all(method=="mean")){
    idnow = which.min(rowSums(metric^2))[1]
  } else if (all(method=="median")){
    idnow = which.min(rowSums(metric))[1]
  }
  return(idnow)
}

# graph5 = list()
# for (i in 1:5){
#   graph5[[i]] = make_star(20, center=sample(1:20, 1))
# }
# testKgraph(graph5, method="ad")
# testKgraph(graph5, method="kw")
# testKgraph(graph5, method="ns")
# testKgraph(graph5, method="steel")