# R/mcmc.rank.R In netrankr: Analyzing Partial Rankings in Networks

#### Documented in mcmc_rank_prob

#' @title Estimate rank probabilities with Markov Chains
#' @description Performs a probabilistic rank analysis based on an almost uniform
#' sample of possible rankings that preserve a partial ranking.
#' @param P P A partial ranking as matrix object calculated with [neighborhood_inclusion]
#'    or [positional_dominance].
#' @param rp Integer indicating the number of samples to be drawn.
#' @details This function can be used instead of [exact_rank_prob]
#' if the number of elements in P is too large for an exact computation. As a rule of thumb,
#' the number of samples should be at least cubic in the number of elements in P.
#' See vignette("benchmarks",package="netrankr") for guidelines and benchmark results.
#' @return
#' \item{expected.rank}{Estimated expected ranks of nodes}
#' \item{relative.rank}{Matrix containing estimated relative rank probabilities:
#' \code{relative.rank[u,v]} is the probability that u is ranked lower than v.}
#'
#' @references Bubley, R. and Dyer, M., 1999. Faster random generation of linear extensions.
#' *Discrete Mathematics*, **201**(1):81-88
#'
#' @seealso [exact_rank_prob], [approx_rank_relative], [approx_rank_expected]
#' @author David Schoch
#' @examples
#' \dontrun{
#' data("florentine_m")
#' P <- neighborhood_inclusion(florentine_m)
#' res <- exact_rank_prob(P)
#' mcmc <- mcmc_rank_prob(P, rp = vcount(g)^3)
#'
#' # mean absolute error (expected ranks)
#' mean(abs(res$expected.rank - mcmc$expected.rank))
#' }
#' @export
mcmc_rank_prob <- function(P, rp = nrow(P)^3) {
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")
}

if (is.null(rownames(P)) & is.null(colnames(P))) {
name_vec <- rownames(P) <- colnames(P) <- paste0("V", seq_len(nrow(P)))
} else {
name_vec <- rownames(P)
}
n.full <- nrow(P)
MSE <- Matrix::which(P == Matrix::t(P) & P == 1, arr.ind = TRUE)
if (length(MSE) >= 1) {
MSE <- t(apply(MSE, 1, sort))
MSE <- MSE[!duplicated(MSE), ]
g <- igraph::graph.empty()
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") } init.rank <- as.vector(igraph::topological.sort(igraph::graph_from_adjacency_matrix(P, "directed"))) if (inherits(P, "Matrix")) { res <- mcmc_rank_sparse(P, init.rank - 1, rp) } else { res <- mcmc_rank_dense(P, init.rank - 1, rp) } res$expected <- res$expected + 1 expected.full <- c(0, n.full) rrp.full <- matrix(0, n.full, n.full) for (i in sort(unique(MSE))) { idx <- which(MSE == i) if (length(idx) > 1) { group.head <- i rrp.full[idx, ] <- do.call(rbind, replicate(length(idx), res$rrp[group.head, MSE], simplify = FALSE))
} else if (length(idx) == 1) {
rrp.full[idx, ] <- res$rrp[i, MSE] } } expected.full <- res$expected[MSE]
for (val in sort(unique(expected.full), decreasing = TRUE)) {
idx <- which(expected.full == val)
expected.full[idx] <- expected.full[idx] + sum(duplicated(MSE[expected.full <= val]))
}
rownames(rrp.full) <- colnames(rrp.full) <- names(expected.full) <- name_vec
res <- list(relative.rank = rrp.full, expected.rank = expected.full)
class(res) <- "netrankr_mcmc"
return(res)
}


## Try the netrankr package in your browser

Any scripts or data that you put into this service are public.

netrankr documentation built on Sept. 27, 2022, 1:07 a.m.