R/approximate.ranks.R

In netrankr: Analyzing Partial Rankings in Networks

#' @title Approximation of expected ranks
#' @description  Implements a variety of functions to approximate expected ranks
#' for partial rankings.
#'
#' @param P A partial ranking as matrix object calculated with [neighborhood_inclusion]
#'    or [positional_dominance].
#' @param method String indicating which method to be used. see Details.
#' @details The \emph{method} parameter can be set to
#' \describe{
#' \item{lpom}{local partial order model}
#' \item{glpom}{extension of the local partial order model.}
#' \item{loof1}{based on a connection with relative rank probabilities.}
#' \item{loof2}{extension of the previous method.}
#' }
#' Which of the above methods performs best depends on the structure and size of the partial
#' ranking. See `vignette("benchmarks",package="netrankr")` for more details.
#' @return A vector containing approximated expected ranks.
#' @author David Schoch
#' @references Brüggemann R., Simon, U., and Mey,S, 2005. Estimation of averaged
#' ranks by extended local partial order models. *MATCH Commun. Math.
#' Comput. Chem.*, 54:489-518.
#'
#' Brüggemann, R. and Carlsen, L., 2011. An improved estimation of averaged ranks
#' of partial orders. *MATCH Commun. Math. Comput. Chem.*,
#' 65(2):383-414.
#'
#' De Loof, L., De Baets, B., and De Meyer, H., 2011. Approximation of Average
#' Ranks in Posets. *MATCH Commun. Math. Comput. Chem.*, 66:219-229.
#'
#' @seealso [approx_rank_relative], [exact_rank_prob], [mcmc_rank_prob]
#' @examples
#' P <- matrix(c(0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, rep(0, 10)), 5, 5, byrow = TRUE)
#' # Exact result
#' exact_rank_prob(P)$expected.rank
#'
#' approx_rank_expected(P, method = "lpom")
#' approx_rank_expected(P, method = "glpom")
#' @export
approx_rank_expected <- function(P, method = "lpom") {
  if (!inherits(P, "Matrix") & !is.matrix(P)) {
    stop("P must be a dense or spare matrix")
  }
  if (!is.binary(P)) {
    stop("P is not a binary matrix")
  }

  # Equivalence classes ------------------------------------------------
  MSE <- Matrix::which((P + Matrix::t(P)) == 2, arr.ind = TRUE)
  if (length(MSE) >= 1) {
    MSE <- t(apply(MSE, 1, sort))
    MSE <- MSE[!duplicated(MSE), ]
    g <- igraph::graph.empty()
    g <- igraph::add.vertices(g, nrow(P))
    g <- igraph::add.edges(g, c(t(MSE)))
    g <- igraph::as.undirected(g)
    MSE <- igraph::clusters(g)$membership
    equi <- which(duplicated(MSE))
    P <- P[-equi, -equi]
  } else {
    MSE <- seq_len(nrow(P))
  }
  if (length(unique(MSE)) == 1) {
    stop("all elements are structurally equivalent and have the same rank")
  }

  # number of Elements
  n <- length(names)

  g <- igraph::graph_from_adjacency_matrix(P, "directed")
  n <- nrow(P)
  if (method == "lpom") {
    sx <- igraph::degree(g, mode = "in")
    ix <- (n - 1) - igraph::degree(g, mode = "all")
    r.approx <- (sx + 1) * (n + 1) / (n + 1 - ix)
    r.approx <- unname(r.approx)
  } else if (method == "glpom") {
    r.approx <- approx_glpom(P)
  } else if (method == "loof1") {
    P <- P + diag(1, n)
    s <- igraph::degree(g, mode = "in")
    l <- igraph::degree(g, mode = "out")
    r.approx <- s + 1
    for (x in 1:n) {
      Ix <- which(P[x, ] == 0 & P[, x] == 0)
      for (y in Ix) {
        approx.rank <- ((s[x] + 1) * (l[y] + 1))
        approx.num.ranks <- ((s[x] + 1) * (l[y] + 1) + (s[y] + 1) * (l[x] + 1))
        r.approx[x] <- r.approx[x] + approx.rank / approx.num.ranks
      }
    }
  } else if (method == "loof2") {
    P <- P + diag(1, n)
    s <- igraph::degree(g, mode = "in")
    l <- igraph::degree(g, mode = "out")
    s.approx <- s
    l.approx <- l
    for (x in 1:n) {
      Ix <- which(P[x, ] == 0 & P[, x] == 0)
      for (y in Ix) {
        s.approx[x] <- s.approx[x] + .sl.approx(s[x], s[y], l[x], l[y])
        l.approx[x] <- l.approx[x] + .sl.approx(s[y], s[x], l[y], l[x])
      }
    }
    r.approx <- s + 1
    s <- s.approx
    l <- l.approx
    for (x in 1:n) {
      Ix <- which(P[x, ] == 0 & P[, x] == 0)
      for (y in Ix) {
        approx.rank <- ((s[x] + 1) * (l[y] + 1))
        approx.num.ranks <- ((s[x] + 1) * (l[y] + 1) + (s[y] + 1) * (l[x] + 1))
        r.approx[x] <- r.approx[x] + approx.rank / approx.num.ranks
      }
    }
  }
  expected.full <- unname(r.approx[MSE])
  for (val in sort(unique(expected.full), decreasing = T)) {
    idx <- which(expected.full == val)
    expected.full[idx] <- expected.full[idx] +
      sum(duplicated(MSE[expected.full <= val]))
  }
  return(expected.full)
}

.sl.approx <- function(sx, sy, lx, ly) {
  ((sx + 1) * (ly + 1)) / ((sx + 1) * (ly + 1) + (sy + 1) * (lx + 1))
}
#############################
#' @title Approximation of relative rank probabilities
#' @description Approximate relative rank probabilities \eqn{P(rk(u)<rk(v))}.
#' In a network context, \eqn{P(rk(u)<rk(v))} is the probability that u is
#' less central than v, given the partial ranking P.
#' @param P A partial ranking as matrix object calculated with [neighborhood_inclusion]
#'    or [positional_dominance].
#' @param iterative Logical scalar if iterative approximation should be used.
#' @param num.iter Number of iterations to be used. defaults to 10 (see Details).
#' @details The iterative approach generally gives better approximations
#' than the non iterative, if only slightly. The default number of iterations
#' is based on the observation, that the approximation does not improve
#' significantly beyond this value. This observation, however, is based on
#' very small networks such that increasing it for large network may yield
#' better results. See `vignette("benchmarks",package="netrankr")` for more details.
#' @author David Schoch
#' @references De Loof, K. and De Baets, B and De Meyer, H., 2008. Properties of mutual
#' rank probabilities in partially ordered sets. In *Multicriteria Ordering and
#' Ranking: Partial Orders, Ambiguities and Applied Issues*, 145-165.
#'
#' @return a matrix containing approximation of relative rank probabilities.
#' \code{relative.rank[i,j]} is the probability that i is ranked lower than j
#' @seealso [approx_rank_expected], [exact_rank_prob], [mcmc_rank_prob]
#' @examples
#' P <- matrix(c(0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, rep(0, 10)), 5, 5, byrow = TRUE)
#' P
#' approx_rank_relative(P, iterative = FALSE)
#' approx_rank_relative(P, iterative = TRUE)
#' @export
approx_rank_relative <- function(P, iterative = TRUE, num.iter = 10) {
  if (!inherits(P, "Matrix") & !is.matrix(P)) {
    stop("P must be a dense or spare matrix")
  }
  if (!is.binary(P)) {
    stop("P is not a binary matrix")
  }

  # Equivalence classes ------------------------------------------------
  MSE <- Matrix::which((P + Matrix::t(P)) == 2, arr.ind = T)

  if (length(MSE) >= 1) {
    MSE <- t(apply(MSE, 1, sort))
    MSE <- MSE[!duplicated(MSE), ]
    g <- igraph::graph.empty()
    g <- igraph::add.vertices(g, nrow(P))
    g <- igraph::add.edges(g, c(t(MSE)))
    g <- igraph::as.undirected(g)
    MSE <- igraph::clusters(g)$membership
    equi <- which(duplicated(MSE))
    P <- P[-equi, -equi]
  } else {
    MSE <- 1:nrow(P)
  }

  if (length(unique(MSE)) == 1) {
    stop("all elements are structurally equivalent and have the same rank")
  }

  relative.rank <- approx_relative(colSums(P), rowSums(P), P, iterative, num.iter)
  mrp.full <- matrix(0, length(MSE), length(MSE))
  for (i in sort(unique(MSE))) {
    idx <- which(MSE == i)
    if (length(idx) > 1) {
      group.head <- i
      mrp.full[idx, ] <- do.call(rbind, replicate(length(idx), relative.rank[group.head, MSE], simplify = FALSE))
    } else if (length(idx) == 1) {
      group.head <- idx
      mrp.full[group.head, ] <- relative.rank[i, MSE]
    }
  }

  diag(mrp.full) <- 0
  return(mrp.full)
}