1 Motivación

En muchos casos la existencia o no de un enlace entre dos nodos en un grafo tiene poca o ninguna incertidumbre para decidir sobre la misma, en estos casos se busca una manera natural de predecir un enlace usando la información observada acerca de la red, de manera tal que se puede plantear la inferencia de los enlaces en tres casos:

  1. Predicción de enlaces dados los vértices y el conocimiento de algunas aristas.

  2. Predicción de enlaces dados los vértices y sus características nodales.

  3. Predicción de enlaces dados algunos vertices y sus características nodales.

Tipos de Inferencia Sobre Topología de la Red

De manera visual se puede identificar los tres tipos de problemas de inferencia de topología de red, para un grafo \(\mathbf{G}\) de juguete. Las aristas existentes se muestran en lineas solidas; las no aristas, en lineas punteadas. Los vértices y aristas observados se muestran en obscuro (es decir, rojo y azul, respectivamente); vértices y aristas no observados, en claro (es decir, rosa y azul claro).

2 Topología de la Red

Sea \(\mathbf{Y}=\left[y_{ij}\right]\) la matriz de adyacencia asociada al grafo \(\mathbf{G}=(\mathbf{V},\mathbf{E})\) de tamaño \(N_v\times N_v\) con \(N_v=\mid\mathbf{V}\mid\) y \(\mathbf{X}=\left(x_1,x_2,\dots,x_{N_v}\right)^T\) el conjunto de características nodales.

Supongase que de las entradas de la matríz \(\mathbf{Y}\) sólo son observadas algunas de ellas denotadas por \(\mathbf{Y^{obs}}\) y se desconocen las restantes denotadas por \(\mathbf{Y^{miss}}\).

El problema principal rádica en predecir las entradas \(\mathbf{Y^{miss}}\) no observadas de \(\mathbf{Y}\) dados \(\mathbf{Y^{obs}=y^{obs}}\) y posiblemente algunos atributos nodales \(\mathbf{X}=\mathbf{x}=\left(x_1,x_2,\dots,x_{N_v}\right)^T\)

Cómo supuesto para la información faltante acerca de la existencia o no de los enlaces es que son por perdida al azar, lo que implica que está información faltante se puede predecir con la información obersvada que está completa, para la cual se proponen dos alternativas:

3 Puntajes

Sean \(i, j\) pares de vértices con enlace existente o no, desconocido y sea \(s(i,j)\) un puntaje que rankea a las diadas cómo medida de potenciales enlaces o no.

3.1 Tipos de Puntaje

  • Negativo de distancia geodésica:

\[ s(i,j)=-\texttt{dist}_G^{obs}(i,j) \]

  • Comparación de las vecindades \(\mathcal{N}_i\) y \(\mathcal{N}_j\) de \(i\) y \(j\):

\[ s(i,j)=\mid\mathcal{N}_i^{obs} \cap \mathcal{N}_j^{obs}\mid \]

  • Coeficiente de Jaccard:

\[ s(i,j)=\frac{\mid\mathcal{N}_i^{obs} \cap \mathcal{N}_j^{obs}\mid}{\mid\mathcal{N}_i^{obs} \cup \mathcal{N}_j^{obs}\mid} \]

  • Coeficiente de Liben-Nowell y Kleinberg:

\[ s(i,j)=\sum_{k\in\left(\mathcal{N}_i^{obs} \cap \mathcal{N}_j^{obs}\right)}\frac{1}{\log\mid\mathcal{N}_k^{obs} \mid} \]

3.1.1 Ejemplo: Red de Juguete

Consideremos el siguiente grafo de juguete \(\mathbf{G}=(\mathbf{V},\mathbf{E})\):

con \(V=\left\lbrace 1,2,3,4,5\right\rbrace\) y \(E=\left\lbrace\lbrace1,2\rbrace,\lbrace1,3\rbrace,\lbrace1,4\rbrace,\lbrace2,3\rbrace,\lbrace2,4\rbrace,\lbrace4,5\rbrace\right\rbrace\)

  • Negativo de distancia geodésica:

\[ s(1,5)=-\texttt{dist}_G^{obs}(1,5)=-2 \]

\[ s(1,2)=-\texttt{dist}_G^{obs}(1,2)=-1 \]

  • Comparación de las vecindades:

\[ s(1,5)=\mid\mathcal{N}_1^{obs} \cap \mathcal{N}_5^{obs}\mid=\mid\lbrace2,3,4\rbrace\cap\lbrace4\rbrace\mid=\mid\lbrace4\rbrace\mid=1 \] \[ s(1,2)=\mid\mathcal{N}_1^{obs} \cap \mathcal{N}_2^{obs}\mid=\mid\lbrace2,3,4\rbrace\cap\lbrace1,3,4\rbrace\mid=\mid\lbrace3,4\rbrace\mid=2 \] - Coeficiente de Jaccard:

\[ s(1,5)=\frac{\mid\mathcal{N}_1^{obs} \cap \mathcal{N}_5^{obs}\mid}{\mid\mathcal{N}_1^{obs} \cup \mathcal{N}_5^{obs}\mid}=\frac{\mid\lbrace2,3,4\rbrace \cap \lbrace4\rbrace\mid}{\mid\lbrace2,3,4\rbrace \cup \lbrace4\rbrace\mid}=\frac{\mid\lbrace4\rbrace\mid}{\mid\lbrace2,3,4\rbrace\mid}=\frac{1}{3} \] \[ s(1,2)=\frac{\mid\mathcal{N}_1^{obs} \cap \mathcal{N}_2^{obs}\mid}{\mid\mathcal{N}_1^{obs} \cup \mathcal{N}_2^{obs}\mid}=\frac{\mid\lbrace2,3,4\rbrace \cap \lbrace1,3,4\rbrace\mid}{\mid\lbrace2,3,4\rbrace \cup \lbrace1,3,4\rbrace\mid}=\frac{\mid\lbrace3,4\rbrace\mid}{\mid\lbrace1,2,3,4\rbrace\mid}=\frac{2}{4}=\frac{1}{2} \]

  • Coeficiente de Liben-Nowell y Kleinberg:

\[ s(1,5)=\sum_{k\in\left(\mathcal{N}_1^{obs} \cap \mathcal{N}_5^{obs}\right)}\frac{1}{\log\mid\mathcal{N}_k^{obs} \mid}=\frac{1}{\log\mid\mathcal{N}_4^{obs} \mid}=\frac{1}{\log\mid\lbrace1,2,5\rbrace\mid}=\frac{1}{\log 3} \] \[ s(1,2)=\sum_{k\in\left(\mathcal{N}_1^{obs} \cap \mathcal{N}_2^{obs}\right)}\frac{1}{\log\mid\mathcal{N}_k^{obs} \mid}=\frac{1}{\log\mid\mathcal{N}_3^{obs} \mid}+\frac{1}{\log\mid\mathcal{N}_4^{obs} \mid}=\frac{1}{\log\mid\lbrace1,2\rbrace\mid}+\frac{1}{\log\mid\lbrace1,2,5\rbrace\mid}=\frac{1}{\log 2}+\frac{1}{\log 3} \]

3.1.2 Ejemplo: Blogs de Política

Red de Blogs políticos franceses. Un enlace indica que al menos uno de los dos blogs pertenecientes a una diada hace referencia al otro en su página web.

Librerias:

library(sand)
library(igraph)
library(ROCR)
library(dplyr)
library(purrr)
library(vioplot)
library(ggplot2)

Número de nodos:

n <- vcount(fblog)
n
## [1] 192

Número de posibles relaciones o diadas:

diadas<-n*(n-1)/2
diadas
## [1] 18336

Matriz de Adyacencia y \(Vec(A_{\Delta_{inf}})\)

A <- as_adjacency_matrix(fblog)
Avec <- A[lower.tri(A)]
Avec %>% head(3)
## [1] 1 0 1

1. Comparación de vecindades:

ncn <- numeric()
for(i in (1:(n-1))){
  .ni <- ego(fblog, 1, i)
  .nj <- ego(fblog, 1, (i+1):n)
  .nbhd.ij <- mapply(intersect, .ni, .nj, SIMPLIFY=FALSE)
  .temp <- unlist(lapply(.nbhd.ij, length)) -2*A[i, (i+1):n]
  ncn <- c(ncn, .temp)
}
ncn %>% head(3)
##             bix.enix.org/      www.arnaudcaron.net/  dominiquevoynet.net/blog 
##                         1                         2                         1
pred_ncn <- prediction(ncn, Avec)
perf_ncn <- performance(pred_ncn, "auc")
paste0('AUC=',slot(perf_ncn, "y.values"),sep='')
## [1] "AUC=0.927517857323709"

2. Negativo distancia geodésica:

dist <- (-1)*igraph::distances(graph=fblog,v=V(fblog),to=V(fblog))
dist <- dist[lower.tri(dist)]
dist %>% head(3)
## [1] -1 -2 -1
pred_dist <- prediction(dist, Avec)
perf_dist  <- performance(pred_dist, "auc")
paste0('AUC=',slot(perf_dist, "y.values"),sep='')
## [1] "AUC=1"

3. Coeficiente de Jaccard

jac <- numeric()
for(i in (1:(n-1))){
  .ni <- ego(fblog, 1, i)
  .nj <- ego(fblog, 1, (i+1):n)
  .nbhd.ij1 <- mapply(base::intersect, .ni, .nj, SIMPLIFY=FALSE)
  .nbhd.ij2 <- mapply(base::union, .ni, .nj, SIMPLIFY=FALSE)
  .temp1 <- unlist(lapply(.nbhd.ij1, length)) -2*A[i, (i+1):n]
  .temp2 <- unlist(lapply(.nbhd.ij2, length)) -2*A[i, (i+1):n]
  jac <- c(jac, .temp1/.temp2)
}
jac %>% head(3)
##             bix.enix.org/      www.arnaudcaron.net/  dominiquevoynet.net/blog 
##                0.06666667                0.33333333                0.11111111
pred_jac <- prediction(jac, Avec)
perf_jac <- performance(pred_jac, "auc")
paste0('AUC=',slot(perf_jac, "y.values"),sep='')
## [1] "AUC=0.929047513636755"

4. Coeficiente de Liben-Nowell y Kleinberg:

LNK<-numeric()
for (i in 1:(n-1)) { 
    ni <- ego(fblog, 1, i)
    nj <- ego(fblog, 1, (i+1):n)
    nbhd.ij <- mapply(base::intersect, ni, nj, SIMPLIFY = FALSE)
    nbhd.ij2<- lapply(nbhd.ij, function(x)
      subset(x, !x %in% c(i, i+1)))
    prueba <- lapply(nbhd.ij2,FUN = 
                     function(x) { apply(as.data.frame(x), MARGIN = 1,FUN = 
                     function(y) {1 / log(length(unlist(ego(fblog, order = 1, nodes = y))) - 1)})})
    temp<-unlist(lapply(prueba,sum))
    LNK<-c(LNK,temp)
}
LNK %>% head(3)
## [1] 0.4342945 0.4342945 0.4342945
pred_LNK <- prediction(LNK, Avec)
perf_LNK <- performance(pred_LNK, "auc")
paste0('AUC=',slot(perf_LNK, "y.values"),sep='')
## [1] "AUC=0.962521353450686"

3.1.3 Resumen Scores Medios

Y Real Vecinos Comúnes Negativo Distancia Geodésica Coeficiente Jaccard Coeficiente Liben-Nowell y Kleinberg
0 1.079681 -2.668915 0.0354973 0.322302
1 8.010482 -1.000000 0.2712412 2.945697

3.1.4 Resumen AUC

Vecinos Comúnes Negativo Distancia Geodésica Coeficiente Jaccard Coeficiente Liben-Nowell y Kleinberg
0.928 1 0.929 0.963

4 Redes de Correlación

\(X\) es una variable aleatoria (continua) correspondiente a los vértices en \(V\). Así,una medida de similaridad entre pares de vértices es

\[ sim(i,j)=\rho_{ij}=corr({X_i,X_j})=\frac{\sigma_{ij}}{\sqrt{\sigma_{ii}\sigma_{jj}}} \]

Es decir, el coeficiente de correlación de Pearson entre \(X_i\) y \(X_j\) expresada en términos de las entradas de la matriz de covarianza \(\Sigma = \{\sigma_{ij}\}\) del vector aleatorio \((X_1,\dots ,X_{N_v})^T\) de los atributos nodales

El grafo \(G= (V,E)\) asociado se denomina grafo (red) de correlación (o covarianza) y se define de forma que \(\{i,j\}\in E\) si y solo si \(\rho_{ij}\) es diferente de cero.

Por lo tanto, para hacer inferencia sobre el grafo de asociación derivado de un conjunto de observaciones \(X_i\) se debe hacer inferencia sobre el conjunto de correlaciones distintas de cero:

\[ H_0 : \rho_{ij} = 0 {~frente~a~} H_1 : \rho_{ij} \neq 0 \]

Estadístico de prueba: transformación z de Fisher sobre el coeficiente de correlación

\[ z_{i,j}= arctanh(r_{ij})=\frac{1}{2} ~~log ~ \frac{1+r_{i,j}}{1-r_{i,j}} \]

donde \(r\) es la versión empírica de \(\rho\).

Distribución nula para la evaluación de la significancia: si el par de variables \((X_i,X_j)\) tienen distribución normal bivariada, la densidad de \(z_{ij}\) bajo la \(H_0 : \rho_{ij} = 0\) está bien aproximada por una variable aleatoria \(N \left(0,\frac{1}{(n-3)}\right)\)

Múltiples pruebas : en este caso, se deben realizar un gran número de pruebas de manera simultánea para \(N_v(N_v-1)/2\) (enlaces potenciales), lo cual implica que se deben ejecutar pruebas múltiples.

4.1 Ejemplo Red E.Coli (Libreria, SAND (R))

Niveles de expresión genética en la bacteria Escherichia coli (E. coli), medidos para 153 genes en cada una de las 40 condiciones experimentales diferentes.

Ecoli.expr es una matriz de 40 por 153 de niveles de expresión génica (logarítmica) en la bacteria Escherichia coli (E. coli), medidos para 153 factores de transcripción bajo cada una de las 40 condiciones experimentales diferentes, promediadas en tres réplicas de cada experimento. Los datos son un subconjunto de los publicados en RegulonDB (http://regulondb.ccg.unam.mx/). Se trata de experimentos de perturbación genética, en los que se “apaga” un gen determinado, para cada uno de los 40 genes diferentes

Red E.Coli

g.RegDB <- graph_from_adjacency_matrix(adjmatrix=regDB.adj, mode = "undirected")
plot(g.RegDB, vertex.label=NA, vertex.size =7.5, vertex.frame.color = "black")
title(main = "Red E.Coli")

vcount(g.RegDB)
## [1] 153
ecount(g.RegDB)
## [1] 209
library(sand)
data(Ecoli.data)
r <- cor(Ecoli.expr)
z <- 0.5 * log((1 + r) / (1 - r))
z.vec <- z[upper.tri(z)]
n <- dim(Ecoli.expr)[1]
corr.pvals <-2 * pnorm(abs(z.vec), 0, sqrt(1 / (n - 3)), lower.tail = FALSE)

Calculados estos \(valores~p\) y comparados con la distribución normal apropiada, verificamos la cantidad de pruebas simultaneas que se están realizando

length(corr.pvals)
## [1] 11628
  • Para las pruebas múltiples, si las decisiones sobres las hipótesis individuales se basan en valores p marginales no ajustados, suele haber una gran probabilidad de que se rechacen algunas de las hipótesis nulas verdaderas, es decir, hay un aumento sustancial del False Discovery Rate (FDR)

La función p.adjust proporciona herramientas para ajustar los \(valores ~p\) para las pruebas múltiples. En este caso, se aplica el ajuste de Benjamini-Hochberg para controlar el FDR

corr.pvals.adj <- p.adjust(corr.pvals, "BH")
table(corr.pvals.adj <0.05) #luego del ajuste
## 
## FALSE  TRUE 
##  6401  5227

Bajo este criterio (de correlación), estamos prediciendo 5227 enlaces que no se pudieron observar. Este numero es muy grande para ser realista y es un orden de magnitud mayor que en la red observada (el orden real es de 209)

Esta prueba tiene como supuesto la normalidad bivariada ¿la normalidad se satisface?

car::qqPlot(z.vec)

## [1]  903 6328
shapiro.test(x = sample(x = z.vec, size = 5000, replace = F))
## 
##  Shapiro-Wilk normality test
## 
## data:  sample(x = z.vec, size = 5000, replace = F)
## W = 0.99854, p-value = 0.0001497

Tanto el QQ-Plot, como la prueba de Shapiro- Wilk, rechazan la normalidad de los scores y podría ser el porqué de los resultados con los ajustes.

Ahora, si la distribución nula no es Normal, ¿entonces, cuál es? Metodología no parametrica, estudia la distribución nula en particular para coeficientes de correlación

library(fdrtool)

#calculo de los valores p, de nuevo
r.vec <- r[upper.tri(r)]
fdr <-fdrtool(x = r.vec,statistic = "correlation",plot = F,verbose = F)
length(fdr$qval)
## [1] 11628
head(fdr$qval) # arrojan no significativos
## [1] 1 1 1 1 1 1
table(fdr$qval <0.05)
## 
## FALSE 
## 11628

CONCLUSIÓN: Según este método, no se establece ningún enlace entre pares de genes

5 Redes de Correlación Parcial

Dos variables aleatorias pueden estar asociadas debido al efecto conjunto de otra (u otras) variables aleatorias de control eliminadas: La correlación parcial mide este tipo de asociación

Usar el coeficiente de correlación directamente sobre las dos variables de interés puede arrojar resultados engañosos si existe otra variable de confusión que esta relacionada con ambas variables.

Dos vértices \(i\) y \(j\) pueden tener atributos altamente correlacionados, \(X_i\) y \(X_j\) porque los vértices influyen entre sí de manera directa. Alternativamente, esta correlación puede ser alta, posiblemente, porque cada uno de ellos está fuertemente relacionado con un tercer vértice k. Por lo tanto \(X_i\) y \(X_j\) están altamente correlacionados con \(X_k\)

La correlación parcial de los atributos \(X_i\) y \(X_j\) de los vértices \(i,j \in V\) definida con respecto a los atributos \(X_{k_1},\dots, X_{k_m}\) de los vértices \(k_1,\dots;k_m \in V\), es la correlación entre \(X_i\) y \(X_j\) que queda después de descontar el efecto de \(X_{k_1},\dots, X_{k_m}\) común a \(X_i\) y \(X_j\)

\[ \rho_{i,j|k} = \frac{\rho_{i,j}-\rho_{i,k}\rho_{j,k}}{\sqrt{(1-\rho_{i,k}^2)(1-\rho_{j,k}^2)}} \]

Para definir los enlaces del grafo inducido por la correlación parcial (dado un valor de \(m\)), se procede de manera análoga que con el grafo de correlación de Pearson, pero reemplazando \(\rho_{ij}\) por \(\rho_{i,j|k}\) (cuando \(m=1\))

El grafo \(G=(V,E)\) correspondiente se denomina grafo (red) de correlación parcial y se define de forma que \(\{i,j\} \in E\) si y solo si \(\rho_{i,j|k}\) es diferente de cero para todo \(S_m \in V^{(m)}_{-(i,j)}\). Donde \(S_m = \{k_1,\dots,k_m\}\) y \(V^{(m)}_{-(i,j)}\) es la colección de todos los subconjuntos no ordenados de \(m\) vértices (distintos)

Por lo tanto, dado un conjunto de observaciones de cada \(X_i\), para hacer inferencia sobre el grafo de correlación parcial se debe hacer inferencia sobre el conjunto de correlaciones distintas de cero:

\[ H_0: \rho_{i,j|k} = 0 ~ para~ algún~ S_m \in V^{(m)}_{-(i,j)} ~~ frente~a ~~ H_0: \rho_{i,j|k} \neq 0 ~ para~ todo~ S_m \in V^{(m)}_{-(i,j)} \]

Estadístico de prueba: transformación \(z\) de Fisher sobre el coeficiente de correlación

\[ z_{i,j}= arctanh(r_{ij})=\frac{1}{2} ~~log ~ \frac{1+r_{i,j|S_m}}{1-r_{i,j|S_m}} \]

donde \(r{i,j|S_m}\) es la versión empírica de \(\rho_{i,j|S_m}\).

Distribución nula para la evaluación de la significancia:bajo el supuesto de normalidad multivariada de \((X_i,\dots,X_j)\) el estadistico de prueba sigue una distribución \(N \left(0,\frac{1}{(n-m-3)}\right)\)

Wille et al. (2004) propone construir una prueba a partir de pruebas de los sub-problemas a través de la agregación:

\[ p_{ij,max}= max\{p_{i,j|S_m} : S_m \in V^{(m)}_{-(i,j)}\} \]

#Tamaño de la muestra
n <- nrow(Ecoli.expr)
n
## [1] 40
#n de vertices
N <- ncol(Ecoli.expr)
N
## [1] 153

Agregación

#correlaciones agregadas y obteniendo el max
m <- 1
pcorr.pvals <- matrix(0, N, N)
for (i in 1:N) {
  for (j in 1:N) {
    rowi <- r[i,-c(i, j)]
    rowj <- r[j,-c(i, j)]
    tmp <- (r[i, j] - rowi * rowj) / sqrt((1 - rowi^2) * (1 - rowj^2)) #correlaciones parciales
    tmp.zvals <-(0.5) * log((1 + tmp) / (1 - tmp)) #estadístico de prueba
    tmp.pvals <-2 * pnorm(abs(tmp.zvals),mean = 0,sd = 1 / sqrt(n - m - 3),lower.tail = FALSE)
    pcorr.pvals[i, j] <-max(tmp.pvals)
  }
}

# valores p ajustados (pruebas multiples)

pcorr.pvals.vec <- pcorr.pvals[lower.tri(pcorr.pvals)]
pcorr.pvals.adj <- p.adjust(pcorr.pvals.vec, "BH")

# relaciones significativas

pcorr.edges <- (pcorr.pvals.adj < 0.05)
table(pcorr.edges)
## pcorr.edges
## FALSE  TRUE 
## 11603    25
fdr <-fdrtool(pcorr.pvals.vec,statistic = "pvalue",plot = FALSE,verbose = F)
pcorr.edges.2 <- (fdr$qval < 0.05)
table(pcorr.edges.2)
## pcorr.edges.2
## FALSE  TRUE 
## 11603    25

Verificamos que sean las mismas conexiones

identical(pcorr.edges,pcorr.edges.2)
## [1] TRUE

Comprobamos con el grafo observado

#grafo inducido
pcorr.A <- matrix(0,N,N)
pcorr.A[lower.tri(pcorr.A)] <- as.numeric(pcorr.pvals.adj < 0.05)
pcorr.A <- pcorr.A + t(pcorr.A)

#comparacion con grafo observado
table(regDB.adj[lower.tri(regDB.adj)],pcorr.A[lower.tri(pcorr.A)])
##    
##         0     1
##   0 11514    21
##   1    89     4

Hay 4 conexiones compartidas, 21 falsos positivos y 89 falsos negativos. Son pocas conexiones en común con el grafo original: esta topologia inducida está basada en las propiedades de los genes (atributos) y con correlacion parcial de orden \(m=1\)

6 Modelo gaussiano

Se denomina modelo gaussiano a un caso particular del método de correlación parcial en el que m = n-2, y los atributos siguen una distribución normal multivariada. En este modelo, la correlación parcial entre los atributos de dos vértices es condicional a la información de los atributos de los demás vértices . Por lo tanto, para \(G = (V, E)\) se tiene que:

En este modelo, los coeficientes de correlación parcial pueden expresarse como:

\[ \rho_{i j \mid V \backslash\{i, j\rangle}=\frac{-\omega_{i j}}{\sqrt{\omega_{i i} \omega_{j j}}}, \]

donde \(\omega_{i j}\) es la entrada \((i, j)\) de \(\Omega=\Sigma^{-1}\), la inversa de la matriz de covarianza \(\Sigma\) del vector \(\left(X_1, \ldots, X_{N_y}\right)^T\) también conocida como matriz de precisión.

El sistema de hipótesis a probar es el siguiente:

\[ H_0: \rho_{i j \mid V \backslash[i, j\}}=0 ~~\text { vs }~~ H_1: \rho_{i j \mid V \backslash\{i, j\}} \neq 0 \]

teniendo en cuenta las correlaciones parciales empíricas \(\hat{\rho}_{i j \mid V \backslash\{i, j\rangle}\). Sin embargo para redes extensas, suele considerarse la relación entre la correlación y la regresión lineal para probar la hipótesis anterior.

Consideremos el vector aleatorio \(\left(X_1, \ldots, X_{N_v}\right)^T\) de los atributos de los vértices que tiene una distribución normal multivariada con media \(\mathbf{0}\), y matriz de covarianzas \(\Sigma\), por lo tanto el vector de atributos \(X_{i}\) dado los valores de los atributos de los demás vértices tiene la siguiente esperanza condicional:

\[\mathbb{E}\left[X_i \mid \mathbf{X}^{(-i)}=\mathbf{x}^{(-i)}\right]=\left(\beta^{(-i)}\right)^T \mathbf{x}^{(-i)} \]

donde

De acuerdo a lo anterior, \((i, j) \in E\) si y solo si \(\beta_i^{(-i)}\) es diferente de 0. Estos coeficientes se estiman con métodos basados en regresión y selección de variables.

Sin embargo como hay \(N_v-1\) variables en la regresión por cada \(X_i\), y además \(n \ll N_v\), se ha propuesto la estimación de una regresión lasso de tal manera que :

\[ \hat{\beta}^{(-i)}=\arg \min _{\beta: \beta_i=0} \sum_{k=1}^n\left(x_{i k}-\left(\beta^{(-i)}\right)^T \mathbf{x}_k^{(-i)}\right)^2+\lambda \sum_{j \neq i}\left|\beta_j^{(-i)}\right| . \]

Esta estrategia no solo estima los coeficientes \(\beta^{(-i)}\), sino que asegura que \(\hat{\beta}_j^{(-i)}=0\) cuando la asociación de \(X_i\) y \(X_j\) es muy débil respecto a la elección del parámetro \(\lambda\).

Para que lo anterior pueda ponerse en práctica se debe:

El parquete huge permite abordar estos dos puntos.

# Librería huge para high-dimensional, undirected graph estimation
library(huge)
library(sand)
data(Ecoli.data)
set.seed(42)
huge.out = huge(Ecoli.expr)
## Conducting Meinshausen & Buhlmann graph estimation (mb)....done
plot(huge.out)

set.seed(42)

# rotational information criterion
huge.opt = huge.select(huge.out, criterion="ric")
## Conducting rotation information criterion (ric) selection....done
## Computing the optimal graph....done
g.huge = graph_from_adjacency_matrix(huge.opt$refit, "undirected")
huge.opt[["opt.lambda"]]
## [1] 155.8486
ecount(g.huge)
## [1] 0
set.seed(42)

# stability selection
huge.opt = huge.select(huge.out, criterion="stars")
g.huge = graph_from_adjacency_matrix(huge.opt$refit, "undirected")
huge.opt[["opt.lambda"]]
ecount(g.huge)
## [1] 759
# Intersección entre el grafo obtenido y el grafo real
g.regDB = graph_from_adjacency_matrix(regDB.adj,"undirected")
intersection(g.regDB, g.huge, byname=FALSE)
## IGRAPH 7c733b9 UN-- 153 22 -- 
## + attr: name (v/c)
## + edges from 7c733b9 (vertex names):
##  [1] yhiW_b3515_at--yhiX_b3516_at tdcA_b3118_at--tdcR_b3119_at
##  [3] rpoE_b2573_at--rpoH_b3461_at rpoD_b3067_at--tyrR_b1323_at
##  [5] rhaR_b3906_at--rhaS_b3905_at nac_b1988_at --putA_b1014_at
##  [7] marA_b1531_at--marR_b1530_at malT_b3418_at--rpoD_b3067_at
##  [9] hns_b1237_at --rcsB_b2217_at hipA_b1507_at--hipB_b1508_at
## [11] himA_b1712_at--himD_b0912_at gutM_b2706_at--srlR_b2707_at
## [13] fruR_b0080_at--mtlR_b3601_at flhD_b1892_at--lrhA_b2289_at
## [15] crp_b3357_at --srlR_b2707_at crp_b3357_at --pdhR_b0113_at
## + ... omitted several edges

7 Inferencia tomográfica de la topología de una red

La inferencia tomográfica se refiere a la inferencia de los componentes interiores de una red (tanto aristas como vértices), a partir de los datos obtenido de los vértices exteriores.

Este problema de inferencia no es fácil, pues sin el número de vértices y aristas internos, ni su forma de conexión, los vértices exteriores pueden generarse de infinitas topologías. Por este motivo, para inferir es necesario hacer suposiciones sobre la estructura interna de la red, la suposición más usual es restringir a la inferencia de redes en forma de árboles.

7.1 Topología de árboles

Se denomina un árbol (no dirigido) \(T=\left(V_T, E_T\right)\) a un grafo conectado sin ciclos

  • Un árbol rotado es un árbol con un vértice \(r \in V_T\) que funciona como raíz del árbol.
  • El conjunto de vértices \(R \subset V_T\) de 1 grado se conocen como hojas.
  • Los vértices tal que \(V \backslash\{\{r\} \cup R\}\) se denominan vértices internos.
  • Las aristas \(E_T\) se conocen como ramas.

7.2 Árbol binario

Un árbol binario es un árbol rotado que en dirección de la raíz a las hojas, cada vértice interno tiene máximo 2 hijos.

Para un grafo arbitrario \(G=(V, E)\), suponga un subconjunto de vértices \(N_l\) con \(n\) observaciones independientes e idénticamente distribuidas de variables aleatorias \(\left\{X_1, \ldots, X_{N_l}\right\}\). Bajo el supuesto de que estos vértices pueden ser identificados con las hojas \(R\) de un árbol \(T\), se quiere encontrar \(\hat{T}\) en el conjunto \(\mathscr{T}_{N_l}\) de todos los árboles binarios con \(N_l\) hojas que mejor explique los datos. Si se tiene conocimiento de la raíz del árbol \(r\), entonces las raíces de los árboles en \(\mathscr{T}_{N_l}\) pueden ser identificados con r.

7.3 Red informática

Consideremos la siguiente topología de una red de computadores (Experimento Castro, Coates, and Nowak (2004))

Coates et al. (2002) propone un esquema de medición denominado sándwich sondeo. Este consiste en enviar una secuencia de tres señales (dos pequeñas y una grande). Las dos señales pequeñas se envían a uno de los vértices de un par de hojas \(\{i, j\}\), mientras que la señal grande se envía al vértice restante. El primer envío es el de una señal pequeña, el segundo es de la grande que induce un retraso en la llegada de la tercera señal (pequeña) y la diferencia de llegada de las señales pequeñas reflejan cuantos caminos se comparten entre el vértice \(i\) y el \(j\).

Los datos de este experimento son los siguientes:

# DelayDiff: Diferencia en microsegundos entre los paquetes pequeños.
# SmallPktDest: Vértice de llegada del paquete pequeño
# BigPktDest: Vértice de llegada del paquete grande.
data(sandwichprobe)
delaydata[1:5, ]
##   DelayDiff SmallPktDest BigPktDest
## 1       757            3         10
## 2       608            6          2
## 3       242            8          9
## 4        84            1          8
## 5      1000            7          3
# Nombres de las máquinas
host.locs
##  [1] "IST"    "IT"     "UCBkly" "MSU1"   "MSU2"   "UIUC"   "UW1"    "UW2"   
##  [9] "Rice1"  "Rice2"
meanmat = with(delaydata, by(DelayDiff, delaydata[, 2:3], mean))

image(log(meanmat + t(meanmat)), xaxt="n", yaxt="n", col=heat.colors(16))

mtext(side=1, text=host.locs, at=seq(0.0,1.0,0.11), las=3)

mtext(side=2, text=host.locs, at=seq(0.0,1.0,0.11), las=1)

7.4 Métodos de inferencia en topologías de árboles

Los dos métodos más desarrollados son:

  1. Agrupamiento jerárquico.
  2. Métodos basados en verosimilitud.

En el agrupamiento jerárquico, se considera el conjunto de \(N_l\) hojas como los objetos a agrupar de acuerdo a una noción de similitud, el árbol resultante de este método es el árbol correspondiente al árbol inferido. En este contexto, se utilizan las \(n\) observaciones de las variables aleatorias \(\left\{X_1, \ldots, X_{N_l}\right\}\) en las hojas que permiten construir una matriz \(N_l \times N_l\) de similitudes.

# Distancia euclidiana entre los retrasos de tiempo
SSDelayDiff = with(delaydata, by(DelayDiff^2, delaydata[,2:3], sum))
x = as.dist(1 / sqrt(SSDelayDiff))
x
##               1            2            3            4            5
## 2  3.889975e-05                                                    
## 3  8.869402e-05 5.861594e-05                                       
## 4  9.163121e-05 5.249794e-05 8.225516e-05                          
## 5  6.988631e-05 5.601092e-05 7.504855e-05 6.184905e-05             
## 6  7.574652e-05 7.211238e-05 7.445629e-05 5.863732e-05 6.004042e-05
## 7  7.038040e-05 8.026426e-05 6.225794e-05 7.967455e-05 7.317794e-05
## 8  8.168043e-05 6.930725e-05 7.970223e-05 5.850447e-05 7.453598e-05
## 9  1.402884e-04 2.118930e-04 1.285978e-04 6.537296e-05 1.725446e-04
## 10 1.447694e-04 2.315999e-04 1.529694e-04 1.179254e-04 2.333809e-04
##               6            7            8            9
## 2                                                     
## 3                                                     
## 4                                                     
## 5                                                     
## 6                                                     
## 7  7.146507e-05                                       
## 8  7.937455e-05 5.075353e-05                          
## 9  1.973667e-04 1.053899e-04 1.582940e-04             
## 10 1.595379e-04 2.404842e-04 1.929773e-04 5.578967e-05
myclust = hclust(x, method="average")
plot(myclust, labels=host.locs, axes=FALSE, ylab=NULL, ann=FALSE)

Referencias

Castro, R. M., M. J. Coates, and R. D. Nowak. 2004. Likelihood Based Hierarchical Clustering. IEEE Transactions on Signal Processing. Vol. 52. 8. https://doi.org/10.1109/TSP.2004.831124.
Coates, Mark, Rui Castro, Robert Nowak, Manik Gadhiok, Ryan King, and Yolanda Tsang. 2002. Maximum Likelihood Network Topology Identification from Edge-Based Unicast Measurements. Sigmetrics Performance Evaluation Review - SIGMETRICS. Vol. 30. https://doi.org/10.1145/511399.511337.
Wille, Anja, Philip Zimmermann, Eva Vranová, Andreas Fürholz, Oliver Laule, Stefan Bleuler, Lars Hennig, et al. 2004. “Sparse Graphical Gaussian Modeling of the Isoprenoid Gene Network in Arabidopsis Thaliana.” Genome Biology 5 (October): R92. https://doi.org/10.1186/gb-2004-5-11-r92.