1
#' @title Indirect relations in a network
2
#' @description Derive indirect relations for a given network.
3
#' Observed relations, like presents or absence of a relation, are commonly not the center
4
#' of analysis, but are transformed in a new set of indirect relation like shortest path
5
#' distances among nodes. These transformations are usually an implicit step when centrality
6
#' indices are used. Making this step explicit gives more possibilities, for example
7
#' calculating partial centrality rankings with [positional_dominance].
8
#' @param g igraph object. The network for which relations should be derived.
9
#' @param type String giving the relation to be calculated. See Details for options.
10
#' @param lfparam Numeric parameter. Only used if type = "dist_lf".
11
#' @param dwparam Numeric parameter. Only used if type = "dist_walk".
12
#' @param netflowmode String, one of raw, frac, or norm. Only used if type = "depend_netflow".
13
#' @param rspxparam Numeric parameter. Only used if type = "depend_rsps" or type = "depend_rspn".
14
#' @param FUN A function that allows the transformation of relations. See Details.
15
#' @param ... Additional arguments passed to FUN.
16
#' @details The `type` parameter has the following options.
17
#'
18
#' \emph{'adjacency'} returns the adjacency matrix of the network.
19
#' 
20
#' \emph{'weights'} returns the weighted adjacency matrix of the network if an edge
21
#' attribute 'weight' is present.
22
#'
23
#' \emph{'dist_sp'} returns shortest path distances between all pairs of nodes.
24
#'
25
#' \emph{'depend_sp'} returns dyadic dependencies
26
#' \deqn{\delta(u,s) = \sum_{t \in V} \frac{\sigma(s,t|u)}{\sigma(s,t)}}
27
#' where \eqn{\sigma(s,t|u)} is the number of shortest paths from s to t that include u and
28
#' \eqn{\sigma(s,t)} is the total number of shortest (s,t)-paths. This relation is used
29
#' for betweenness-like centrality indices.
30
#'
31
#' \emph{'walks'} returns walk counts between pairs of nodes, usually they are
32
#' weighted decreasingly in their lengths or other properties which can be done by adding
33
#' a function in \code{FUN}.  See [transform_relations] for options.
34
#'
35
#' \emph{'dist_resist'} returns the resistance distance between all pairs of nodes.
36
#'
37
#' \emph{'dist_lf'} returns a logarithmic forest distance \eqn{d_\alpha(s,t)}. The logarithmic forest
38
#' distances form a one-parametric family of distances, converging to shortest path distances as \eqn{\alpha -> 0}
39
#' and to the resistance distance as \eqn{\alpha -> \infty}. See (Chebotarev, 2011) for more details.
40
#' The parameter `lfparam` can be used to tune \eqn{\alpha}.
41
#'
42
#'\emph{'dist_walk'} returns the walk distance \eqn{d_\alpha^W(s,t)} between nodes. The walk distances form a one-parametric
43
#'family of distances, converging to shortest path distances as \eqn{\alpha -> 0} and to longest
44
#'walk distances for \eqn{\alpha -> \infty}. Walk distances contain the logarithmic forest
45
#'distances as a special case. See (Chebotarev, 2012) for more details.
46
#'
47
#' \emph{'dist_rwalk'} returns the expected length of a random walk between two nodes. For more
48
#' details see (Noh and Rieger, 2004)
49
#'
50
#' \emph{'depend_netflow'} returns dependencies based on network flow (See Freeman et al.,1991). 
51
#' If `netflowmode="raw"`, the function returns 
52
#' \deqn{\delta(u,s) = \sum_{t \in V} f(s,t,G)-f(s,t,G-v)} 
53
#' where f(s,t,G) is the maximum flow from s to t in G and f(s,t,G-v) in G without the node v.
54
#' For `netflowmode="frac"` it returns dependencies in the form, similar to shortest path dependencies:
55
#'\deqn{\delta(u,s) = \sum_{t \in V} \frac{f(s,t,G)-f(s,t,G-v)}{f(s,t,G)}} 
56
#'
57
#' \emph{'depend_curflow'} returns pairwise dependencies based on current flow. The relation is
58
#' based on the same idea as 'depend_sp' and 'depend_netflow'. However, instead of considering
59
#' shortest paths or network flow, the current flow (or equivalent: random walks) between nodes
60
#' are of interest. See (Newman, 2005) for details.
61
#'
62
#' \emph{'depend_exp'} returns pairwise dependencies based on 'communicability':
63
#' \deqn{\delta(u,s)=\sum_{t \in V} \frac{exp(A)_{st}-exp(A+E(u))_{st}}{exp(A)_{st}},}
64
#' where E(u) has nonzeros only in row and column u, and 
65
#' in this row and column has -1 if A has +1. See (Estrada et al., 2009) for additional details.
66
#'
67
#' \emph{'depend_rsps'}. Simple randomized shortest path dependencies. 
68
#' The simple RSP dependency of a node u with respect to absorbing paths from s to t, 
69
#' is defined as the expected number of visits through u over all s-t-walks. The
70
#' parameter `rspxparam` is the "inverse temperature parameter". 
71
#' If it converges to infinity, only shortest paths are considered and the expected
72
#' number of visits to a node on a shortest path approaches the probability of
73
#' following that particular path. When the parameter converges to zero, then the 
74
#' dependencies converge to the expected number of visits to a node over all absorbing
75
#' walks with respect to the unbiased random walk probabilities. This means for undirected networks,
76
#' that the relations converge to adjacency. See (Kivimäki et al., 2016) for details.
77
#'
78
#' \emph{'depend_rspn'} Net randomized shortest path dependencies. 
79
#' The parameter `rspxparam` is the "inverse temperature parameter". The asymptotic 
80
#' for the infinity case are the same as for 'depend_rsps'. If the parameter approaches zero, then
81
#' it converges to 'depend_curflow'. The net randomized shortest path dependencies
82
#' are closely related to the random walk interpretation of current flows.
83
#'  See (Kivimäki et al., 2016) for technical details.
84
#'
85
#'
86
#' The function \code{FUN} is used to transform the indirect
87
#' relation. See [transform_relations] for predefined functions and additional help.
88
#'
89
#' @return A matrix containing indirect relations in a network.
90
#' @author David Schoch
91
#' @seealso [aggregate_positions] to build centrality indices, [positional_dominance] to derive dominance relations
92
#' @references Chebotarev, P., 2012. The walk distances in graphs. *Discrete Applied Mathematics*, 160(10), pp.1484-1500.  
93
#' 
94
#' Chebotarev, P., 2011. A class of graph-geodetic distances generalizing the shortest-path and
95
#' the resistance distances. *Discrete Applied Mathematics* 159,295-302.
96
#'
97
#' Noh, J.D. and Rieger, H., 2004. Random walks on complex networks. *Physical Review Letters*, 92(11), p.118701.
98
#'
99
#' Freeman, L.C., Borgatti, S.P., and White, D.R., 1991.
100
#' Centrality in Valued Graphs: A Measure of Betweenness Based on Network Flow. *Social Networks* 13(2), 141-154.
101
#'
102
#' Newman, M.E., 2005. A measure of betweenness centrality based on random walks. *Social Networks*, 27(1), pp.39-54.
103
#'
104
#' Estrada, E., Higham, D.J., and Hatano, N., 2009.
105
#' Communicability betweenness in complex networks. *Physica A* 388,764-774.
106
#'
107
#' Kivimäki, I., Lebichot, B., Saramäki, J., and Saerens, M., 2016.
108
#' Two betweenness centrality measures based on Randomized Shortest Paths
109
#' *Scientific Reports* 6: 19668
110
#' @examples
111
#' library(igraph)
112
#' g <- graph.empty(n=11,directed = FALSE)
113
#' g <- add_edges(g,c(1,11,2,4,3,5,3,11,4,8,5,9,5,11,6,7,6,8,
114
#'                    6,10,6,11,7,9,7,10,7,11,8,9,8,10,9,10))
115
#'
116
#' #shortest path distances
117
#' D <- indirect_relations(g,type = "dist_sp")
118
#'
119
#' #inverted shortest path distances
120
#' D <- indirect_relations(g,type = "dist_sp", FUN = dist_inv)
121

122
#' #shortes path dependencies (used for betweenness)
123
#' D <- indirect_relations(g,type = "depend_sp")
124
#'
125
#' #walks attenuated exponentially by their length
126
#' W <- indirect_relations(g,type = "walks",FUN = walks_exp)
127
#'
128
#' @export
129
indirect_relations <- function(g, 
130
                               type = "dist_sp",
131
                               lfparam = NULL, 
132
                               dwparam = NULL,
133
                               netflowmode = "",
134
                               rspxparam = NULL,
135
                               FUN = identity, ...) {
136 2
  if (type == "dependencies") {
137 0
    warning('type="dependencies" is deprecated. Using "depend_sp" instead.\n')
138 0
    type <- "depend_sp"
139
  }
140 2
  if (type == "geodesic") {
141 0
    warning('type="geodesic" is deprecated. Using "dist_sp" instead.\n')
142 0
    type <- "dist_sp"
143
  }
144 2
  if (type == "resistance") {
145 0
    warning('type="resistance" is deprecated. Using "dist_resist" instead.\n')
146 0
    type <- "dist_resist"
147
  }
148 2
  if (type == "identity") {
149 0
    warning('type="identity" is deprecated. Using "adjacency" instead.\n')
150 0
    type <- "adjacency"
151
  }
152 2
  if (type == "dist_sp") {
153 2
    rel <- igraph::distances(g, mode = "all")
154 2
    rel <- FUN(rel, ...)
155 2
  } else if (type == "weights") {
156 0
    if(is.null(igraph::get.edge.attribute(g,"weight"))){
157 0
      warning('no weight attribute present. using "adjacency" instead.\n')
158 0
      rel <- igraph::get.adjacency(g, type = "both", sparse = FALSE, attr = NULL)
159 0
      rel <- FUN(rel, ...)
160 0
      diag(rel) <- 0
161
    } else {
162 0
      rel <- igraph::get.adjacency(g, type = "both", sparse = FALSE, attr = "weight")
163 0
      rel <- FUN(rel, ...)
164 0
      diag(rel) <- 0
165
    }
166 2
  } else if (type == "adjacency") {
167 2
      rel <- igraph::get.adjacency(g, type = "both", sparse = FALSE, attr = NULL)
168 2
      rel <- FUN(rel, ...)
169 2
      diag(rel) <- 0
170 2
  } else if (type == "depend_sp") {
171 2
      adj <- lapply(igraph::get.adjlist(g), function(x) x - 1)
172 2
      rel <- dependency(adj)
173 2
  } else if (type == "walks") {
174 2
      eigen.decomp <- eigen(igraph::get.adjacency(g, type = "both"))
175 2
      lambda <- eigen.decomp$values
176 2
      X <- eigen.decomp$vectors
177 2
      rel <- X %*% diag(FUN(lambda, ...)) %*% t(X)
178 2
  } else if (type == "dist_resist") {
179 2
      L <- igraph::graph.laplacian(g, sparse = FALSE)
180 2
      n <- igraph::vcount(g)
181 2
      A <- L + matrix(1 / n, n, n)
182 2
      C <- solve(A)
183 2
      rel <- resistanceDistance(C, n)
184 2
      rel <- FUN(rel, ...)
185 2
  } else if (type == "dist_lf") {
186 2
      if (is.null(lfparam)) {
187 2
        stop('argument "lfparam" is missing for "dist_lf", with no default')
188
      }
189 2
      rel <- log_forest_fct(g, lfparam)
190 2
      rel <- FUN(rel, ...)
191 2
  } else if (type == "dist_walk") {
192 2
      if (is.null(dwparam)) {
193 2
        stop('argument "dwparam" is missing for "dist_walk", with no default')
194
      }
195 2
      rel <- dist_walk_fct(g, dwparam)
196 2
      rel <- FUN(rel, ...)
197 2
  } else if (type == "depend_netflow") {
198 2
      if (netflowmode == "" | !netflowmode %in% c("raw", "frac", "norm")) {
199 0
        stop('netflowmode must be one of"raw","frac","norm"\n')
200
      }
201
      # if (netflowmode == "norm") {
202
      #   warning('"norm" not supported yet. Using "frac" instead.\n')
203
      #   netflowmode <- "frac"
204
      # }
205 2
      rel <- depend_netflow_fct(g, netflowmode)
206 2
      rel <- FUN(rel, ...)
207 2
  } else if (type == "depend_exp") {
208 2
      rel <- depend_exp_fct(g)
209 2
      rel <- FUN(rel, ...)
210 2
  } else if (type == "depend_rsps") {
211 2
      if (is.null(rspxparam)) {
212 0
        stop('argument "rspxparam" is missing for "depend_rsps", with no default')
213
      }
214 2
      rel <- depend_rsps_fct(g, rspxparam)
215 2
      rel <- FUN(rel, ...)
216 2
  } else if (type == "depend_rspn") {
217 2
      if (is.null(rspxparam)) {
218 0
        stop('argument "rspxparam" is missing for "depend_rspn", with no default')
219
      }
220 2
      rel <- depend_rspn_fct(g, rspxparam)
221 2
      rel <- FUN(rel, ...)
222 2
  } else if (type == "depend_curflow") {
223 2
      rel <- depend_curflow_fct(g)
224 2
      rel <- FUN(rel, ...)
225 2
  } else if (type == "dist_rwalk") {
226 2
      rel <- dist_rwalk_fct(g)
227 2
      rel <- FUN(rel, ...)
228
  } else {
229 0
    stop(paste(type, "is not defined as indirect relation."))
230
  }
231
  #add names if present
232 2
  if(is.null(rownames(rel)) & !is.null(igraph::V(g)$name) ){
233 0
    rownames(rel) <- colnames(rel) <- igraph::V(g)$name
234
  }
235 2
  return(rel)
236
}
237

238
# -------------------------------------------------------------------------------
239

240
log_forest_fct <- function(g, lfparam) {
241 2
  n <- igraph::vcount(g)
242 2
  gamma <- log(exp(1) + lfparam ^ (2 / n))
243

244 2
  L <- igraph::graph.laplacian(g, sparse = FALSE)
245 2
  I <- diag(1, n)
246 2
  Q <- solve(I + lfparam * L)
247

248 2
  if (lfparam == 1) {
249 0
    H <- gamma * log(Q)
250
  } else {
251 2
    H <- gamma * (lfparam - 1) * logb(Q, lfparam)
252
  }
253 2
  rel <- 0.5 * (diag(H) %*% t(rep(1, n)) + rep(1, n) %*% t(diag(H))) - H
254 2
  return(rel)
255
}
256

257
depend_netflow_fct <- function(g, netflowmode) {
258 2
  n <- igraph::vcount(g)
259 2
  mflow <- matrix(0, n, n)
260
  # maxflow
261 2
  for (s in 1:n) {
262 2
    for (t in 1:n) {
263 2
      if (s != t) {
264 2
        mflow[s, t] <- igraph::graph.maxflow(g, s, t)$value
265
      }
266
    }
267
  }
268 2
  if (netflowmode == "norm") {
269 0
    flo <- mflow
270 0
    diag(flo) <- 0
271 0
    maxoflo <- rep(0, n)
272 0
    for (i in 1:n) maxoflo[i] <- sum(mflow[-i, -i])
273
  }
274 2
  flow_smat <- matrix(0, n, n)
275 2
  for (i in 1:n) {
276 2
    g_i <- igraph::delete.vertices(g, i)
277 2
    for (s in 1:n) {
278 2
      for (t in 1:n) {
279 2
        if (i != s & s != t & i != t) {
280 2
          flow <- igraph::graph.maxflow(g_i, s - (s > i), t - (t > i))$value
281 2
          flow_smat[i, s] <- switch(netflowmode,
282 2
            raw = flow_smat[i, s] + mflow[s, t] - flow,
283 2
            norm = flow_smat[i, s] + mflow[s, t] - flow,
284 2
            frac = flow_smat[i, s] + (mflow[s, t] - flow) / mflow[s, t]
285
          )
286
        }
287
      }
288
    }
289
  }
290 2
  if (netflowmode == "norm") {
291 0
    flow_smat <- flow_smat / maxoflo * 2
292
  }
293 2
  return(flow_smat)
294
}
295

296
depend_exp_fct <- function(g) {
297 2
  A <- igraph::get.adjacency(g, "both", sparse = F)
298 2
  eigen_A <- eigen(A)
299 2
  n <- nrow(A)
300 2
  expA <- eigen_A$vectors %*% diag(exp(eigen_A$values)) %*% t(eigen_A$vectors)
301 2
  C <- (n - 1) ^ 2 - (n - 1)
302 2
  combet <- matrix(0, n, n)
303 2
  for (i in 1:n) {
304 2
    E <- matrix(0, n, n)
305 2
    E[which(A[, i] == 1), i] <- -1
306 2
    E[i, which(A[i, ] == 1)] <- -1
307 2
    E <- A + E
308 2
    eigen_E <- eigen(E)
309 2
    expE <- eigen_E$vectors %*% diag(exp(eigen_E$values)) %*% t(eigen_E$vectors)
310 2
    expE <- (expA - expE) / expA
311 2
    expE[i, ] <- 0
312 2
    expE[, i] <- 0
313 2
    diag(expE) <- 0
314 2
    combet[i, ] <- 1 / C * rowSums(expE)
315
  }
316

317 2
  return(combet)
318
}
319

320
depend_rsps_fct <- function(g, rspxparam) {
321 2
  n <- igraph::vcount(g)
322 2
  I <- diag(1, n)
323

324 2
  A <- igraph::get.adjacency(g, sparse = F)
325 2
  D <- diag(1 / igraph::degree(g))
326 2
  P_ref <- D %*% A
327 2
  C <- 1 / A
328 2
  C[is.infinite(C)] <- 0
329 2
  W <- P_ref * exp(-rspxparam * C)
330

331 2
  Z <- solve((I - W), I)
332 2
  Zdiv <- 1 / Z
333 2
  bet.mat <- t(sapply(1:n, function(x)
334 2
    (Z[x, ] * t(Zdiv)) %*% Z[, x] - sum(Z[x, ] * diag(Zdiv) * Z[, x])))
335 2
  diag(bet.mat) <- 0
336

337 2
  return(bet.mat)
338
}
339

340
depend_rspn_fct <- function(g, rspxparam) {
341 2
  n <- igraph::vcount(g)
342 2
  I <- diag(1, n)
343

344 2
  A <- igraph::get.adjacency(g, sparse = F)
345 2
  D <- diag(1 / igraph::degree(g))
346 2
  P_ref <- D %*% A
347 2
  C <- 1 / A
348 2
  C[is.infinite(C)] <- 0
349 2
  W <- P_ref * exp(-rspxparam * C)
350

351 2
  Z <- solve((I - W), I)
352 2
  Zdiv <- 1 / Z
353 2
  adj <- igraph::get.adjlist(g, "all")
354 2
  A <- lapply(adj, function(x) as.vector(x) - 1)
355 2
  bet.mat <- dependRspn(A, Z, Zdiv, W, n)
356 2
  return(bet.mat)
357
}
358

359
depend_curflow_fct <- function(g) {
360 2
  n <- igraph::vcount(g)
361 2
  D <- diag(igraph::degree(g))
362 2
  A <- igraph::get.adjacency(g, sparse = F)
363 2
  L <- D - A
364

365 2
  Tmat <- solve(L[1:(n - 1), 1:(n - 1)])
366 2
  Tmat <- rbind(cbind(Tmat, 0), 0)
367

368 2
  el <- igraph::get.edgelist(g, names = FALSE)
369 2
  m <- igraph::ecount(g)
370 2
  el <- el - 1L
371

372 2
  bet.mat <- dependCurFlow(Tmat, el, m, n)
373 2
  return(bet.mat)
374
}
375

376
dist_rwalk_fct <- function(g) {
377 2
  n <- igraph::vcount(g)
378 2
  A <- igraph::get.adjacency(g, sparse = FALSE)
379 2
  M <- A / rowSums(A)
380 2
  e <- rep(1, n - 1)
381 2
  H <- matrix(0, n, n)
382 2
  for (j in 1:n) {
383 2
    Mj <- M[-j, -j]
384 2
    Hij <- solve(diag(e) - Mj) %*% e
385 2
    H[j, -j] <- Hij # transposed to original (to fit framework with rowSums)
386
  }
387 2
  return(H)
388
}
389

390
dist_walk_fct <- function(g,dwparam) {
391 2
  n <- igraph::vcount(g)
392 2
  A <- igraph::get.adjacency(g, sparse = FALSE)
393 2
  bigeig <- eigen(A,only.values = TRUE)$values[1]
394 2
  t <- (1/dwparam+bigeig)^(-1)
395
  # if(dwparam > (1/bigeig)){
396
  #   stop(paste0("dwparam to large. To ensure convergence, use a value greater 0 and less than 1/",bigeig))
397
  # }
398 2
  I <- diag(1,n)
399 2
  Rt <- solve((I-t*A))
400 2
  Ht <- log(Rt)
401 2
  Dt <- 0.5 * (diag(Ht) %*% t(rep(1, n)) + rep(1, n) %*% t(diag(Ht))) - Ht
402 2
  alpha <- dwparam#(1/dwparam - bigeig)^(-1)
403 2
  if(alpha!=1){
404 2
    gamma <- log(exp(1)+alpha^(2/n))*(alpha-1)/log(alpha)
405
  } else {
406 0
    gamma <- log(exp(1)+1)
407
  }
408 2
  Dt <- gamma * Dt
409 2
  return(Dt)
410
}

Read our documentation on viewing source code .

Loading