1
#' @title Probabilistic centrality rankings 
2
#' @description  Performs a complete and exact rank analysis of a given partial ranking.
3
#' This includes rank probabilities, relative rank probabilities and expected ranks.
4
#' 
5
#' @importFrom Rcpp evalCpp
6
#' @useDynLib netrankr
7
#' 
8
#' @param P A partial ranking as matrix object calculated with [neighborhood_inclusion]
9
#'    or [positional_dominance].
10
#' @param only.results Logical. return only results (default) or additionally 
11
#'     the ideal tree and lattice if \code{FALSE}.
12
#' @param verbose Logical. should diagnostics be printed. Defaults to \code{FALSE}.
13
#' @param force Logical. If \code{FALSE} (default), stops the analysis if the partial
14
#'     ranking has more than 40 elements and less than 0.4 comparable pairs. 
15
#'     Only change if you know what you are doing. 
16
#' @details The function derives rank probabilities from a given partial ranking 
17
#' (for instance returned by [neighborhood_inclusion] or [positional_dominance]). This includes the
18
#' calculation of expected ranks, (relative) rank probabilities and the number of possible rankings.
19
#'  Note that the set of rankings grows exponentially in the number of elements and the exact 
20
#'  calculation becomes infeasible quite quickly and approximations need to be used.
21
#'  See `vignette("benchmarks")` for guidelines and [approx_rank_relative], 
22
#'  [approx_rank_expected], and [mcmc_rank_prob] for approximative methods.
23
#' @return 
24
#' \item{lin.ext}{Number of possible rankings that extend `P`.}
25
#' \item{mse}{Array giving the equivalence classes of `P`.}
26
#' \item{rank.prob}{Matrix containing rank probabilities: \code{rank.prob[u,k]} is the probability that u has rank k.}
27
#' \item{relative.rank}{Matrix containing relative rank probabilities: \code{relative.rank[u,v]} is the probability that u is ranked lower than v.}
28
#' \item{expected.rank}{Expected ranks of nodes in any centrality ranking.}
29
#' \item{rank.spread}{Standard deviation of the ranking probabilities.}
30
#' \item{topo.order}{Random ranking used to build the lattice of ideals (if \code{only.results = FALSE}).}
31
#' \item{tree}{Adjacency list (incoming) of the tree of ideals (if \code{only.results = FALSE}).}
32
#' \item{lattice}{Adjacency list (incoming) of the lattice of ideals (if \code{only.results = FALSE}).}
33
#' \item{ideals}{List of order ideals (if \code{only.results = FALSE}).}
34
#' In all cases, higher numerical ranks imply a higher position in the ranking. That is,
35
#' the lowest ranked node has rank 1.
36
#' @author David Schoch, Julian Müller  
37
#' @references De Loof, K. 2009. Efficient computation of rank probabilities in posets. 
38
#' *Phd thesis*, Ghent University.
39
#' 
40
#' De Loof, K., De Meyer, H. and De Baets, B., 2006. Exploiting the
41
#'lattice of ideals representation of a poset. *Fundamenta Informaticae*, 71(2,3):309-321.
42
#'
43
#' @seealso [approx_rank_relative], [approx_rank_expected], [mcmc_rank_prob]
44
#' @examples
45
#' 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)
46
#' P
47
#' res <- exact_rank_prob(P)
48
#' 
49
#' #a warning is displayed if only one ranking is possible
50
#' tg <- threshold_graph(20,0.2)
51
#' P <- neighborhood_inclusion(tg)
52
#' res <- exact_rank_prob(P)
53
#' @export
54
exact_rank_prob <- function(P, only.results = T, verbose = F, force = F) {
55
    # Check for names ------------------------------------------------
56 2
    if (is.null(rownames(P)) & is.null(colnames(P))) {
57 2
        rownames(P) <- colnames(P) <- paste0("V",1:nrow(P))
58
    } 
59 2
    n_full <- nrow(P)
60 2
    P_full <- P
61
    # Equivalence classes ------------------------------------------------
62 2
    MSE <- which((P + t(P)) == 2, arr.ind = T)
63 2
    if (length(MSE) >= 1) {
64 2
        MSE <- t(apply(MSE, 1, sort))
65 2
        MSE <- MSE[!duplicated(MSE), ]
66 2
        g <- igraph::graph.empty()
67 2
        g <- igraph::add.vertices(g, nrow(P))
68 2
        g <- igraph::add.edges(g, c(t(MSE)))
69 2
        g <- igraph::as.undirected(g)
70 2
        MSE <- igraph::clusters(g)$membership
71 2
        equi <- which(duplicated(MSE))
72 2
        P <- P[-equi, -equi]
73
    } else {
74 2
        MSE <- 1:nrow(P)
75
    }
76 2
    if (is.null(nrow(P))) {
77 2
        warning("all elements are structurally equivalent and have the same rank")
78 2
        return()
79
    }
80
    # names <- rownames(P)
81
    # number of Elements
82 2
    nElem <- nrow(P)
83
    
84
    # check for linear order ---------------------------------------------
85 2
    if (comparable_pairs(P) == 1) {
86 2
        warning("P is already a ranking.\nExpected Ranks correspond to the only possible ranking.")
87 2
        expected_full <- rank(colSums(P_full), ties.method = "max")
88 2
        rank.spread_full <- rep(0, nrow(P_full))
89 2
        mrp_full <- P_full
90 2
        mrp_full[mrp_full == t(mrp_full)] <- 0
91 2
        rp_full  <-  matrix(0, nrow(P_full), nrow(P_full))
92 2
        for (i in 1:nrow(P_full)) {
93 2
            rp_full[i, expected_full[i]] <- 1
94
        }
95
        # add names
96 2
        rownames(rp_full) <- rownames(mrp_full) <- colnames(mrp_full) <- rownames(P_full)
97 2
        names(expected_full) <- names(rank.spread_full) <- rownames(P_full)
98 2
        colnames(rp_full) <- 1:ncol(rp_full)
99 2
        return(list(lin.ext = 1, 
100 2
                    mse = MSE, 
101 2
                    rank.prob = rp_full, 
102 2
                    relative.rank = mrp_full, 
103 2
                    expected.rank = expected_full, 
104 2
                    rank.spread = rank.spread_full))
105
    }
106
    
107
    # sanity check if applicable ------------------------------------------------
108 2
    if (nrow(P) > 40 & comparable_pairs(P) < 0.4 & force == F) {
109 0
        stop("Input data too big. Use approximations or set force=T if you know what you are doing")
110
    }
111
    # Prepare Data structures---------------------
112 2
    topo.order <- as.vector(igraph::topological.sort(igraph::graph_from_adjacency_matrix(P, "directed")))
113
    
114 2
    P <- P[topo.order, topo.order]
115 2
    ImPred <- igraph::get.adjlist(igraph::graph_from_adjacency_matrix(P, "directed"), "in")
116 2
    ImPred <- lapply(ImPred, function(x) as.vector(x) - 1)
117
    
118 2
    ImSucc <- igraph::get.adjlist(igraph::graph_from_adjacency_matrix(P, "directed"), "out")
119 2
    ImSucc <- lapply(ImSucc, function(x) as.vector(x) - 1)
120
    # TREEOFIDEALS ----------------------------------------------------
121 2
    if (verbose == TRUE) {
122 0
        print("building tree of ideals")
123
    }
124 2
    tree <- treeOfIdeals(ImPred)
125 2
    nIdeals <-  length(tree$label)
126 2
    if (verbose == TRUE) {
127 0
        print("tree of ideals built")
128
    }
129 2
    Ek <-  sapply(0:(nElem - 1), function(x) {
130 2
        which(tree$label == x) - 1
131
    })
132
    # tree$child=lapply(tree$child,function(x) {idx=order(tree$label[x+1],decreasing=T);x[idx]})
133 2
    if (verbose == TRUE) {
134 0
        print("building lattice of Ideals")
135
    }
136 2
    latofI <-  LatticeOfIdeals(tree$child, tree$parent, Ek, nElem, nIdeals)
137
    
138 2
    if (verbose == TRUE) {
139 0
        print("lattice of ideals built")
140
    }
141 2
    ideallist <-  listingIdeals(ImSucc, nElem, nIdeals)
142
    # ideallist=lapply(ideallist,sort)
143
    
144 2
    if (verbose == TRUE) {
145 0
        print("ideals listed")
146
    }
147
    
148 2
    if (verbose == TRUE) {
149 0
        print(paste("No of ideals:", nIdeals))
150 0
        print("Calculating Rank Probabilities")
151
    }
152
    
153 2
    res <- rankprobs(latofI, ideallist, nElem, nIdeals)
154 2
    if (verbose == TRUE) {
155 0
        print(paste("No. of possible Rankings: ", res$linext))
156
    }
157 2
    res$rp <- res$rp[order(topo.order), ]
158 2
    res$mrp <- res$mrp[order(topo.order), order(topo.order)]
159
    ############################### END
160 2
    expected <- res$rp %*% 1:nElem
161 2
    rank.spread <- rowSums((matrix(rep(1:nElem, each = nElem), nElem, nElem) - c(expected))^2 * res$rp)
162 2
    expected <- c(expected)
163
    ############################### Insert equivalent nodes again
164 2
    rp_full <- matrix(0, n_full, ncol(res$rp))
165 2
    mrp_full <- matrix(0, n_full, n_full)
166 2
    expected_full <- c(0, n_full)
167 2
    rank.spread_full <- rep(0, n_full)
168 2
    for (i in sort(unique(MSE))) {
169 2
        idx <- which(MSE == i)
170 2
        if (length(idx) > 1) {
171 0
            group.head <- i
172 0
            rp_full[idx, ] <- do.call(rbind, replicate(length(idx), res$rp[group.head, ], simplify = FALSE))
173 0
            mrp_full[idx, ] <- do.call(rbind, replicate(length(idx), res$mrp[group.head, MSE], simplify = FALSE))
174 0
            rank.spread_full[idx] <- rank.spread[group.head]
175 2
        } else if (length(idx) == 1) {
176 2
            rp_full[idx, ] <- res$rp[i, ]
177 2
            mrp_full[idx, ] <- res$mrp[i, MSE]
178 2
            rank.spread_full[idx] <- rank.spread[i]
179
        }
180
    }
181 2
    expected_full <- expected[MSE]
182 2
    for (val in sort(unique(expected_full), decreasing = T)) {
183 2
        idx <- which(expected_full == val)
184 2
        expected_full[idx] <- expected_full[idx] + sum(duplicated(MSE[expected_full <= val]))
185
    }
186
    # add names
187 2
    rownames(rp_full) <- rownames(mrp_full) <- colnames(mrp_full) <- rownames(P_full)
188 2
    names(expected_full) <- names(rank.spread_full) <- rownames(P_full)
189 2
    colnames(rp_full) <- 1:ncol(rp_full)
190
    ############################### 
191 2
    if (only.results) {
192 2
        return(list(lin.ext = res$linext, 
193 2
                    mse = MSE, 
194 2
                    rank.prob = rp_full, 
195 2
                    relative.rank = t(mrp_full), 
196 2
                    expected.rank = expected_full, 
197 2
                    rank.spread = sqrt(rank.spread_full)))
198
    } else {
199 2
        return(list(lin.ext = res$linext, 
200 2
                    mse = MSE, 
201 2
                    rank.prob = rp_full, 
202 2
                    relative.rank = t(mrp_full), 
203 2
                    expected.rank = expected_full, 
204 2
                    rank.spread = sqrt(rank.spread_full), 
205 2
                    topo.order = topo.order, 
206 2
                    tree = tree, 
207 2
                    lattice = latofI, 
208 2
                    ideals = ideallist))
209
    }
210
}
211

Read our documentation on viewing source code .

Loading