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 1
  if (type == "dependencies") {
137 0
    warning('type="dependencies" is deprecated. Using "depend_sp" instead.\n')
138 0
    type <- "depend_sp"
139
  }
140 1
  if (type == "geodesic") {
141 0
    warning('type="geodesic" is deprecated. Using "dist_sp" instead.\n')
142 0
    type <- "dist_sp"
143
  }
144 1
  if (type == "resistance") {
145 0
    warning('type="resistance" is deprecated. Using "dist_resist" instead.\n')
146 0
    type <- "dist_resist"
147
  }
148 1
  if (type == "identity") {
149 0
    warning('type="identity" is deprecated. Using "adjacency" instead.\n')
150 0
    type <- "adjacency"
151
  }
152 1
  if (type == "dist_sp") {
153 1
    rel <- igraph::distances(g, mode = "all")
154 1
    rel <- FUN(rel, ...)
155 1
  } 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 1
  } else if (type == "adjacency") {
167 1
      rel <- igraph::get.adjacency(g, type = "both", sparse = FALSE, attr = NULL)
168 1
      rel <- FUN(rel, ...)
169 1
      diag(rel) <- 0
170 1
  } else if (type == "depend_sp") {
171 1
      adj <- lapply(igraph::get.adjlist(g), function(x) x - 1)
172 1
      rel <- dependency(adj)
173 1
  } else if (type == "walks") {
174 1
      eigen.decomp <- eigen(igraph::get.adjacency(g, type = "both"))
175 1
      lambda <- eigen.decomp$values
176 1
      X <- eigen.decomp$vectors
177 1
      rel <- X %*% diag(FUN(lambda, ...)) %*% t(X)
178 1
  } else if (type == "dist_resist") {
179 1
      L <- igraph::graph.laplacian(g, sparse = FALSE)
180 1
      n <- igraph::vcount(g)
181 1
      A <- L + matrix(1 / n, n, n)
182 1
      C <- solve(A)
183 1
      rel <- resistanceDistance(C, n)
184 1
      rel <- FUN(rel, ...)
185 1
  } else if (type == "dist_lf") {
186 1
      if (is.null(lfparam)) {
187 1
        stop('argument "lfparam" is missing for "dist_lf", with no default')
188
      }
189 1
      rel <- log_forest_fct(g, lfparam)
190 1
      rel <- FUN(rel, ...)
191 1
  } else if (type == "dist_walk") {
192 1
      if (is.null(dwparam)) {
193 1
        stop('argument "dwparam" is missing for "dist_walk", with no default')
194
      }
195 1
      rel <- dist_walk_fct(g, dwparam)
196 1
      rel <- FUN(rel, ...)
197 1
  } else if (type == "depend_netflow") {
198 1
      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 1
      rel <- depend_netflow_fct(g, netflowmode)
206 1
      rel <- FUN(rel, ...)
207 1
  } else if (type == "depend_exp") {
208 1
      rel <- depend_exp_fct(g)
209 1
      rel <- FUN(rel, ...)
210 1
  } else if (type == "depend_rsps") {
211 1
      if (is.null(rspxparam)) {
212 0
        stop('argument "rspxparam" is missing for "depend_rsps", with no default')
213
      }
214 1
      rel <- depend_rsps_fct(g, rspxparam)
215 1
      rel <- FUN(rel, ...)
216 1
  } else if (type == "depend_rspn") {
217 1
      if (is.null(rspxparam)) {
218 0
        stop('argument "rspxparam" is missing for "depend_rspn", with no default')
219
      }
220 1
      rel <- depend_rspn_fct(g, rspxparam)
221 1
      rel <- FUN(rel, ...)
222 1
  } else if (type == "depend_curflow") {
223 1
      rel <- depend_curflow_fct(g)
224 1
      rel <- FUN(rel, ...)
225 1
  } else if (type == "dist_rwalk") {
226 1
      rel <- dist_rwalk_fct(g)
227 1
      rel <- FUN(rel, ...)
228
  } else {
229 0
    stop(paste(type, "is not defined as indirect relation."))
230
  }
231
  #add names if present
232 1
  if(is.null(rownames(rel)) & !is.null(igraph::V(g)$name) ){
233 0
    rownames(rel) <- colnames(rel) <- igraph::V(g)$name
234
  }
235 1
  return(rel)
236
}
237

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

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

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

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

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

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

317 1
  return(combet)
318
}
319

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

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

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

337 1
  return(bet.mat)
338
}
339

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

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

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

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

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

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

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

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

390
dist_walk_fct <- function(g,dwparam) {
391 1
  n <- igraph::vcount(g)
392 1
  A <- igraph::get.adjacency(g, sparse = FALSE)
393 1
  bigeig <- eigen(A,only.values = TRUE)$values[1]
394 1
  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 1
  I <- diag(1,n)
399 1
  Rt <- solve((I-t*A))
400 1
  Ht <- log(Rt)
401 1
  Dt <- 0.5 * (diag(Ht) %*% t(rep(1, n)) + rep(1, n) %*% t(diag(Ht))) - Ht
402 1
  alpha <- dwparam#(1/dwparam - bigeig)^(-1)
403 1
  if(alpha!=1){
404 1
    gamma <- log(exp(1)+alpha^(2/n))*(alpha-1)/log(alpha)
405
  } else {
406 0
    gamma <- log(exp(1)+1)
407
  }
408 1
  Dt <- gamma * Dt
409 1
  return(Dt)
410
}

Read our documentation on viewing source code .

Loading