1 Introducción

En algunos casos, hay poca o ninguna incertidumbre al evaluar si existe o no un enlace entre dos vértices y podemos evaluar exhaustivamente la incidencia entre pares de vértices.

En otros casos, es posible que solo se tenga información sobre el estado de algunos de los posibles enlaces de la red, pero no de todos. Además, es posible que ni siquiera se tenga la capacidad de evaluar directamente si hay un enlace presente o no.

Puede ser que solo se puedan observar los atributos de los vértice o de las díadas que predicen el estado de los enlace. En tal caso el objetivo consiste en construir una red a partir de los datos disponibles como una tarea de inferencia estadística.

Inferencia sobre el estatus de aristas no observadas (presente o ausente): 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. 3. Predicción de enlaces dados algunos vértices y sus características.

2 Predicción de enlaces

\(\mathbf{Y}\): matriz de adyacencia de \(n\times n\) asociada con una red \(G=(V,E)\).

\(\mathbf{x}\): atributos de vértices.

Algunas de las entradas de \(\mathbf{Y}\) están observadas, \(\mathbf{Y}^{\text{obs}}\), mientras que las otras son “faltantes” o “perdidas”, \(\mathbf{Y}^{\text{miss}}\).

El objetivo es predecir (imputar) enlaces potenciales entre pares de vértices (presentes o ausentes) utilizando la información proporcionada por un subconjunto de enlaces observados y, si están disponibles, los atributos de los vértices.

El objetivo es predecir (imputar) los enlaces perdidos en \(\mathbf{Y}^{\text{miss}}\), dado que \(\mathbf{Y}^{\text{obs}} = \mathbf{y}^{\text{obs}}\) y posiblemente \(\mathbf{X} = \mathbf{x}\).

Un supuesto común es que la información faltante sobre la presencia/ausencia de las aristas es perdida al azar (missing at random, MAR). Esto significa que los datos perdidos se pueden explicar por medio de las variables donde si hay información completa.

Alternativas:

Modelos.
Puntajes (scores).

2.1 Puntajes

Para cada par de vértices \(i\) y \(j\) cuyo estado de enlace es desconocido, se calcula un puntaje \(s(i, j)\).

Se puede pronosticar por medio de un umbral o por medio de un ranking.

Se han propuesto muchos puntajes en la literatura:

Negativo de la distancia del camino más corto entre \(i\) y \(j\): \[ s(i,j) = -\textsf{dist}^{\text{obs}}(i,j) \] El signo negativo está presente para que los valores de puntuación más altos indiquen que es más “probable” que los pares de vértices compartan una arista.
Comparación de las vecindades \(\mathcal{N}_i\) y \(\mathcal{N}_j\) de \(i\) y \(j\): \[ s(i,j) = |\mathcal{N}^{\text{obs}}_i \cap \mathcal{N}^{\text{obs}}_j| \]
Coeficiente de Jaccard: \[ s(i,j) = \frac{|\mathcal{N}^{\text{obs}}_i \cap \mathcal{N}^{\text{obs}}_j|}{|\mathcal{N}^{\text{obs}}_i \cup \mathcal{N}^{\text{obs}}_j|} \]
Coeficiente de Liben-Nowell and Kleinberg: \[ s(i,j) = \sum_{k\in\mathcal{N}^{\text{obs}}_i \cap \mathcal{N}^{\text{obs}}_j} \frac{1}{\log|\mathcal{N}^{\text{obs}}_k|} \]

2.2 Ejemplo: Blogs de política

Red de blogs políticos franceses clasificados por el proyecto Observatoire Presidentielle en cuanto a afiliación política. Un enlace indica que al menos uno de los dos blogs hacía referencia al otro en su página web. Disponible en la libreria sand de R.

Fuente original: http://observatoire-presidentielle.fr/

suppressMessages(suppressWarnings(library(igraph)))
# datos
data(fblog)
# orden
vcount(fblog)

## [1] 192

# tamaño
ecount(fblog)

## [1] 1431

# dirigida?
is_directed(fblog)

## [1] FALSE

# data
suppressMessages(suppressWarnings(library(sand)))
n <- vcount(fblog)
A <- as.matrix(as_adjacency_matrix(fblog))
colnames(A) <- rownames(A) <- rep("", n)
# n diadas
n*(n-1/2)

## [1] 36768

ncn <- NULL
for (i in 1:(n-1)){
  ni <- ego(graph = fblog, order = 1, nodes = i)
  nj <- ego(graph = fblog, order = 1, nodes = (i+1):n)
  Nij <- mapply(FUN = intersect, ni, nj, SIMPLIFY = FALSE)
  temp <- unlist(lapply(Nij, length)) - 2*A[i, (i+1):n]
  ncn <- c(ncn, temp) 
}
length(ncn)

## [1] 18336

ncn[1:10]

##                     
## 1 2 1 2 2 2 0 0 0 0

# grafico 
suppressMessages(suppressWarnings(library(vioplot)))
Avec <- A[lower.tri(A)]
vioplot(ncn[Avec == 0], ncn[Avec == 1], names = c("Sin enlace", "Con enlace"), col = "mistyrose")
title(ylab="Number of Common Neighbors")

# prediccion: leave-one-out cross-validation
suppressMessages(suppressWarnings(library(ROCR)))
pred <- prediction(predictions = ncn, labels = Avec)
perf <- performance(prediction.obj = pred, measure = "auc")
slot(object = perf, name = "y.values")

## [[1]]
## [1] 0.9275179

3 Red de asociaciones

Las aristas surgen debido a un nivel suficiente de asociación entre ciertos atributos de los vértices incidentes (association networks).

Cada uno de los \(n\) vértices en \(V\) tiene asociado un vector \(\mathbf{x}\) de \(m\) atributos.

Existe una medida de similaridad \(\textbf{sim}(i,j)\) entre los vértices \(i\) y \(j\) en \(V\), junto con una idea clara de cuáles valores de \(\textbf{sim}\) constituyen un nivel de asociación no trivial entre \(i\) y \(j\).

\(\textbf{sim}\) no es observable directamente, pero los atributos \(\mathbf{x}\) contienen suficientemente información para hacer inferencias sobre \(\textbf{sim}\).

Dos medidas de asociación lineales comunes y populares son:

Correlación.
Correlación parcial.
Modelo de asociación Gaussiano.

3.1 Ejemplo: redes reguladoras de genes

Ecoli.expr: niveles de expresión génica (log nivel de actividad) en la bacteria Escherichia coli (E. coli), medidos para 153 genes en cada una de las 40 condiciones experimentales diferentes (promediado sobre tres repeticiones de cada experimento).

regDB.adj: matriz de adyacencia de relaciones regulatorias observadas entre genes de E. coli (http://regulondb.ccg.unam.mx/) al mismo tiempo que se recolectaron los datos experimentales.

Escherichia coli es una bacteria miembro de la familia de las enterobacterias y forma parte de la microbiota del tracto gastrointestinal de animales homeotermos, como por ejemplo el ser humano.

https://tecnosolucionescr.net/templates/yootheme/cache/31_Ecoli-18ff42b9.png

J. Faith, B. Hayete, J. Thaden, I. Mogno, J. Wierzbowski, G. Cottarel, S. Kasif, J. Collins, T. Gardner: Large-scale mapping and validation of Escherichia coli transcriptional regulation from a compendium of expression profiles. PLoS Biol. 5(1), e8 (2007).

# data
data(Ecoli.data)
dim(Ecoli.expr)

## [1]  40 153

dim(regDB.adj)

## [1] 153 153

# Rowv = NA : remover dendograma
# 40 experimentos: filas
# 153 genes: columnas
heatmap(x = scale(Ecoli.expr), Rowv = NA)

El dendograma evidencia que algunos genes muestran un comportamiento similar en ciertos subconjuntos de experimentos. Este hecho sugiere el valor de construir una red de asociación entre genes.

# red observada
suppressMessages(suppressWarnings(library(igraph)))
g.regDB <- graph_from_adjacency_matrix(adjmatrix = regDB.adj, mode = "undirected")
# orden
vcount(g.regDB)

## [1] 153

# tamaño
ecount(g.regDB)

## [1] 209

# dirigida?
is_directed(g.regDB)

## [1] FALSE

plot(g.regDB, vertex.size = 3, vertex.label = NA)

Los pares de genes reguladores observados están contenidos en gran medida dentro de un solo componente conectado.

3.2 Redes de correlación

Una medida estándar de similitud entre pares de vértices es \[ \textsf{sim}(i,j) = \textsf{corr}(X_i,X_j) = \rho_{i,j} \] donde \(X_i\) es una variable aleatoria continua asociada con el vértice \(i\) y \(\textsf{corr}\) es el coef. de correlación de Pearson.

El grafo \(G=(V,E)\) correspondiente se denomina grafo (red) de correlación (o covarianza) y se define de forma que \(\{i,j\}\in E\) si y sólo si \(\rho_{i,j}\neq 0\).

Por lo tanto, dado un conjunto de observaciones de cada \(X_i\), para hacer inferencia sobre el grafo de asociación se debe hacer inferencia sobre el conjunto de correlaciones distintas de cero: \[ H_0:\rho_{i,j} = 0\qquad\text{frente a}\qquad H_1:\rho_{i,j} \neq 0\,. \] Estadístico de prueba: estadistico de prueba: transformacion \(z\) de Fisher sobre coef. de correlación, \[ z_{i,j} = \textsf{arctanh}(r_{i,j}) = \frac12\,\log\frac{1+r_{i,j}}{1-r_{i,j}} \] donde \(r\) es la versión empírica de \(\rho\).

# coef. de correlacion muestral
r <- cor(Ecoli.expr)
# estadistico de prueba
z <- 0.5*log((1 + r)/(1 - r))
dim(r)

## [1] 153 153

dim(z)

## [1] 153 153

Distribución del estadístico bajo \(H_0\) (paradigma frecuentista): Si \((X_i,X_j)\sim\textsf{N}_2\) y \(n\) es lo suficientemente “grande”, entonces \(z_{i,j}\sim\textsf{N}_1(0,1/(n-3))\) bajo \(H_0\).

# tamaño de la muestra
n <- nrow(Ecoli.expr)
n

## [1] 40

# n de vertices
N <- ncol(Ecoli.expr)
N

## [1] 153

# estadisticos de prueba
z.vec <- z[upper.tri(z)]
length(z.vec)

## [1] 11628

N*(N-1)/2

## [1] 11628

# valores p
corr.pvals <- 2*pnorm(q = abs(z.vec), mean = 0, sd = sqrt(1/(n-3)), lower.tail = F)
# n de pruebas multiples
length(corr.pvals)

## [1] 11628

Pruebas múltiples (si las decisiones sobre las hipótesis individuales se basan en los valores p marginales no ajustados, entonces suele haber una gran probabilidad de que se rechacen algunas de las hipótesis nulas verdaderas, es decir, hay un aumento sustancial del FDR [False Descovery Rate]).

Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society Series B, 57, 289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x. https://www.jstor.org/stable/2346101.

# valores p ajustados: Benjamini-Hochberg
corr.pvals.adj <- p.adjust(p = corr.pvals, method = "BH")
length(corr.pvals.adj)

## [1] 11628

head(corr.pvals)

## [1] 2.428765e-03 1.576109e-02 1.075615e-05 5.599030e-03 9.408595e-06
## [6] 2.353950e-02

head(corr.pvals.adj)

## [1] 7.771515e-03 3.718383e-02 7.305638e-05 1.590267e-02 6.542958e-05
## [6] 5.189967e-02

# n de relaciones significativas
table(corr.pvals.adj < 0.05)

## 
## FALSE  TRUE 
##  6401  5227

Este número es demasiado grande para ser realista y un orden de magnitud mayor que en la red observada.

La normalidad se satisface?

# normalidad?
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.99741, p-value = 1.625e-07

Cuál es entonces la distribución nula?

suppressMessages(suppressWarnings(library(fdrtool)))
# https://cran.r-project.org/web/packages/fdrtool/index.html
# Estimates both tail area-based false discovery rates (Fdr) as well as local false discovery rates (fdr) for a variety of null models (p-values, z-scores, correlation coefficients, t-scores). 
# The proportion of null values and the parameters of the null distribution are adaptively estimated from the data. 
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)

## [1] 1 1 1 1 1 1

table(fdr$qval < 0.05)

## 
## FALSE 
## 11628

De acuerdo con este método no se estable ningún enlace entre pares de genes.

3.3 Redes de corrleación parcial

La correlación parcial mide la asociación entre dos variables aleatorias, con el efecto de un conjunto de variables aleatorias de control eliminadas.

Si se quiere encontrar hasta qué punto existe una relación entre dos variables de interés, usar el coeficiente de correlación directamente puede arrojar resultados engañosos si hay otra variable de confusión que está relacionada con ambas variables de interés.

Dos vértices \(i\) y \(j\) pueden tener atributos altamente correlacionados \(X_i\) y \(X_j\) porque los vértices de alguna manera influyen entre sí de manera directa. Alternativamente, su correlación puede ser alta principalmente porque por ejemplo cada uno de ellos está fuertemente influenciado por un tercer vértice \(k\), y 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},\ldots, X_{k_m}\) de los vértices \(k_1,\ldots, k_m \in V\{i, j\}\), es la correlación entre \(X_i\) y \(X_j\) que queda después de descontar el efecto de\(X_{k_1},\ldots, X_{k_m}\) común a \(X_i\) y \(X_j\).

Cuando \(m=0\), el coef. de correlación parcial se reduce al coef. de correlación de Pearson ordinario.

Cuando \(m=1\), el coef. de correlación parcial se puede calcular como \[ \rho_{i,j\mid k} = \frac{\rho_{i,j} - \rho_{i,k}\rho_{j,k}}{\sqrt{(1-\rho_{i,k}^2)(1-\rho_{j,k}^2)}} \]

El caso general se puede ver en: T. Anderson, An Introduction to Multivariate Statistical Analysis, Second Edition. New York: John Wiley & Sons, Inc., 1984.

Para definir los enlaces del grafo inducido por la correlación parcial (dado un valor de \(m\)) se procede de modo similar que con el grafo de correlación de Pearson, pero reemplazando \(\rho_{i,j}\) por \(\rho_{i,j\mid 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 sólo si \(\rho_{i,j\mid S_m}\neq 0\) para todo \(S_m\in V^{(m)}_{-(i,j)}\), donde \(S_m=\{k_1,\ldots, k_m \}\) y \(V^{(m)}_{-(i,j)}\) es la colección de todos los subconjuntos no ordenados de \(m\) vértices (distintos) de \(V-\{i,j\}\).

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\mid S_m} = 0\text{ para agún } S_m\in V^{(m)}_{-(i,j)}\qquad\text{frente a}\qquad H_1:\rho_{i,j\mid S_m} \neq 0\text{ para todo } S_m\in V^{(m)}_{-(i,j)}\,. \] Estadístico de prueba: estadistico de prueba: transformacion \(z\) de Fisher sobre coef. de correlación, \[ z_{i,j} = \textsf{arctanh}(r_{i,j}) = \frac12\,\log\frac{1+r_{i,j\mid S_m}}{1-r_{i,j\mid S_m}} \] donde \(r_{\mid S_m}\) es la versión empírica de \(\rho_{\mid S_m}\). Bajo el supuesto de normalidad multivariada de \(X_i,X_j,X_{k_1},\ldots,X_{k_m}\) el estadístico de prueba sigue una distribución normal con media 0 y varianza \(1/(1-n-m-3)\).

Wille et al. (2004) propone construir una prueba a partir de pruebas de los sub-problemas a través de la agregación: \[ p_{i,j,\max} = \max\left\{p_{i,j\mid S_m}: S_m\in V^{(m)}_{-(i,j)}\right\} \] A. Wille, P. Zimmermann, E. Vránova, A. Fürholz, O. Laule, S. Bleuler, L. Hennig, A. Prelic, P. Rohr, L. Thiele et al., “Sparse graphical Gaussian modeling of the isoprenoid gene network in arabidopsis thaliana,” Genome Biology, vol. 5, no. 11, p. R92, 2004.

# tamaño de la muestra
n <- nrow(Ecoli.expr)
n

## [1] 40

# n de vertices
N <- ncol(Ecoli.expr)
N

## [1] 153

# correlaciones agregadas con max
m <- 1
# r: correlaciones simples (de Pearson)
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)) # corr. parciales
     tmp.zvals <- 0.5*log((1 + tmp)/(1 - tmp)) # est. de prueba
     tmp.pvals <- 2 * pnorm(q = abs(tmp.zvals), mean = 0, sd = 1/ sqrt(n-m-3), lower.tail = F)
     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(p = pcorr.pvals.vec, method = "BH")

# relaciones significativas
pcorr.edges <- pcorr.pvals.adj < 0.05
table(pcorr.edges)

## pcorr.edges
## FALSE  TRUE 
## 11603    25

# usando fdr
fdr <- fdrtool(x = pcorr.pvals.vec, statistic = "pvalue", plot = F, verbose = F)
pcorr.edges.2 <- fdr$qval < 0.05
table(pcorr.edges)

## pcorr.edges
## FALSE  TRUE 
## 11603    25

# los mismas conecciones?
identical(pcorr.edges, pcorr.edges.2)

## [1] TRUE

# 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

g.regDB <- graph_from_adjacency_matrix(adjmatrix = regDB.adj, mode = "undirected")
g.pcorr <- graph_from_adjacency_matrix(adjmatrix = pcorr.A,   mode = "undirected")
intersection(g.regDB, g.pcorr, byname = FALSE)

## IGRAPH 4f782fc UN-- 153 4 -- 
## + attr: name (v/c)
## + edges from 4f782fc (vertex names):
## [1] yhiW_b3515_at--yhiX_b3516_at rhaR_b3906_at--rhaS_b3905_at
## [3] marA_b1531_at--marR_b1530_at gutM_b2706_at--srlR_b2707_at

Cuatro de estas 25 aristas están entre los establecidos en la red observada.

3.4 Modelo de asociación Gaussiano

Caso particular de correlación parcial cuando \(m=n-2\) y los atributos tiene distribución normal multivariada.

La correlación parcial \(\rho_{i,j\mid V_{-(i,j)}}\) entre los atributos de dos vértices es condicional a la información de los atributos de los demás vértices.

Bajo el supuesto de normalidad multivariada de los atributos, \(\rho_{i,j\mid V_{-(i,j)}} = 0\) si y sólo si \(X_i\) y \(X_j\) son condicionalmente independientes dado el esto de los atributos.

El grafo \(G=(V,E)\) correspondiente se denomina grafo de independencia condicional (o modelo gráfico gaussiano o grafo de concentración) y se define de forma que \(\{i,j\}\in E\) si y sólo si \(\rho_{i,j\mid V_{-(i,j)}}\neq 0\)

En este caso, los coeficientes de correlación parcial se pueden expresar como \[ \rho_{i,j\mid V_{-(i,j)}} = \frac{-\omega_{i,j}}{\sqrt{\omega_{i,i}}\sqrt{\omega_{j,j}}} \] donde \([\omega_{i,j}] = \Omega = \Sigma^{-1}\) es la matriz de precisión (conentración) de \((X_1,\ldots,X_n)\).

La construcción de \(G\) se denomina problema de selección de covarianza (covariance selection problem), donde se deben probar los sistemas \[ H_0:\rho_{i,j\mid V_{-(i,j)}} = 0\qquad\text{frente a}\qquad H_1:\rho_{i,j\mid V_{-(i,j)}} \neq 0\,. \]

Si \(X=(X_1,\ldots,X_n)\sim\textsf{N}_n(0,\Sigma)\), entonces \[ \mathsf{E}\left(X_i\mid X^{(-i)} = x^{(-i)}\right) = \beta^{(-i)\,\mathsf{T}} x^{(-i)} \] donde el superíndice \((-i)\) hace referencia a todas las entradas del vector menos la componente \(i\).

Las entradas de \(\beta^{(-i)}=\left[\beta^{(-i)}_j\right]\) se pueden expresar como \[ \beta^{(-i)}_j = \frac{-\omega_{i,j}}{\omega_{i,i}} \,. \]

Por lo tanto \(\{i,j\}\in E\) si y sólo si \(\beta^{(-i)}_j \neq 0\). Lo que quiere decir que para inferir \(G\) se debe hacer inferencia sobre coeficientes de la regresión condicional (métodos de estimación basados en regresión y selección de variables).

Como hay \(n-1\) variables regresoras para cada \(X_i\) con \(n=|V\), y además, el tamaño de la muestra puede ser menor que \(n\) se recomienda tener un enfoque de regresión penalizada (Lasso: realiza simultáneamente estimación y selección de variables, ): \[ \hat\beta^{(-i)} = \arg\min_{\beta:\beta_i=0}\sum_{k} \left( x_{i,k} - \beta^{(-i)\,\mathsf{T}} x^{(-i)} \right)^2 + \lambda \sum_{j=i} \beta^{(-i)}_j \] Se fuerza a \(\beta^{(-i)}_j\) donde la asociación entre \(X_i\) y \(X_j\) se considera demasiado débil.

Consideraciones iportantes en la práctica:

Normalidad multivariada.
Elección de la penalización \(\lambda\).

# estimacion
# huge: High-Dimensional Undirected Graph Estimation
# Data preprocessing, graph estimation, and model selection techniques 
# Nonparanormal(npn) transformation is applied to help relax the normality assumption. 
# Graph structure is estimated by Meinshausen-Buhlmann graph estimation or the graphical lasso
library(huge)

## Warning: package 'huge' was built under R version 4.1.2

set.seed(42)
huge.out <- huge(Ecoli.expr)

## Conducting Meinshausen & Buhlmann graph estimation (mb)....done

# rotational information criterion (prone to under-selection)
huge.opt <- huge.select(est = huge.out, criterion = "ric")

## Conducting rotation information criterion (ric) selection....done
## Computing the optimal graph....done

g.huge <- graph_from_adjacency_matrix(adjmatrix = huge.opt$refit, mode = "undirected")
ecount(g.huge)

## [1] 0

# stability selection (prone to over-selection)
huge.opt <- huge.select(est = huge.out, criterion = "stars")

## Conducting Subsampling....in progress:5% 
Conducting Subsampling....in progress:10% 
Conducting Subsampling....in progress:15% 
Conducting Subsampling....in progress:20% 
Conducting Subsampling....in progress:25% 
Conducting Subsampling....in progress:30% 
Conducting Subsampling....in progress:35% 
Conducting Subsampling....in progress:40% 
Conducting Subsampling....in progress:45% 
Conducting Subsampling....in progress:50% 
Conducting Subsampling....in progress:55% 
Conducting Subsampling....in progress:60% 
Conducting Subsampling....in progress:65% 
Conducting Subsampling....in progress:70% 
Conducting Subsampling....in progress:75% 
Conducting Subsampling....in progress:80% 
Conducting Subsampling....in progress:85% 
Conducting Subsampling....in progress:90% 
Conducting Subsampling....in progress:95% 
Conducting Subsampling....in progress:100% 
Conducting Subsampling....done.

g.huge <- graph_from_adjacency_matrix(huge.opt$refit, "undirected")
summary(g.huge)

## IGRAPH 5b09f21 U--- 153 759 --

# inter con correlacion parcial
intersection(g.pcorr, g.huge)

## IGRAPH 5b0b51e U--- 153 25 -- 
## + edges from 5b0b51e:
##  [1] 145--146 144--146 112--125 112--113 109--138 108--135  97--111  96--119
##  [9]  92--107  87--104  86-- 87  84--129  81--137  72--141  70-- 75  60--127
## [17]  46-- 77  42-- 43  38--153  37-- 98  27--123  21-- 33  12--135   9-- 90
## [25]   3-- 60

# inter con red observada
intersection(g.regDB, g.huge, byname=FALSE)

## IGRAPH 5b0c639 UN-- 153 22 -- 
## + attr: name (v/c)
## + edges from 5b0c639 (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

# ---

4 Inferencia tomográfica

Inferencia sobre los componentes “interiores” de la red. Ambos vértices y aristas, a partir de datos obtenidos en algún subconjunto de vértices “exteriores”.

Para un conjunto dado de observaciones, es posible que existan muchas topologías de red que podrían haberlas generado, dado que sin restricciones en aspectos como el número de vértices y aristas internas, y la manera en que se conectan entre sí, es difícil establecer soluciones específicas.

Una simplificación estructural clave ha sido la restricción a la inferencia de redes en forma de árboles.

Un árbol (no dirigido) \(T = (V_T, E_T)\) es un grafo conectado sin ciclos.

Un árbol enraizado (rooted tree) es un árbol en el que se destaca un vértice \(r ∈ V_T\):

El subconjunto de vértices \(R \subset V_T\) de grado uno se denomina hojas.
Los vértices en \(V_T - \{\{r\} \cup R\}\) se llaman vértices internos.
Las aristas en \(E_T\) se denominan ramas.

En las redes informáticas, los computadores de escritorio y portátiles son instancias de vértices “exteriores”, mientras que los routers de Internet a los que no se tiene acceso son vértices “interiores”.

R. Castro, M. Coates, and R. Nowak, “Likelihood-based hierarchical clustering,” IEEE Transactions on Signal Processing, vol. 52, no. 8, pp. 2308–2321, 2004.

Las máquinas de acceso por los usuarios son “exteriores” (conocidas). Los otros vértices (como los routers) y las ramas son todos “internos” dado que normalmente no los conocería un usuario estándar.

Un mapa de Internet completo no está disponible.

4.1 Objetivo

Un árbol binario es un árbol en el que al moverse desde la raíz hacia las hojas, cada vértice interno tiene como máximo dos hijos.

Para un conjunto de \(N_l\) vértices se tienen observaciones acerca de iid de variables aleatorias \(\{X_1,\ldots, X_{N_l}\}\). Bajo el supuesto de que estos vértices se pueden identificar con las hojas \(R\) de un árbol \(T\), el objetivo consiste en encontrar el árbol binario \(\hat{T}\) con \(N_l\) hojas que mejor explique los datos.

El objetivo radica en el aprendizaje de la estructura interna de árboles binarios.

Inferencia filogenética: describir relaciones evolutivas entre las especies biológicas.
Análisis de redes informáticas: caracterizar el árbol a lo largo del cual fluye el tráfico desde una dirección de Internet de origen hasta un conjunto de direcciones de destino.

4.2 Ejemplo: Redes informáticas

Se envía una secuencia de tres paquetes de información, dos pequeños y uno grande. El paquete grande se envía en segundo lugar, después del primer paquete pequeño y antes del segundo paquete pequeño.

Los dos paquetes pequeños se envían a uno de los vértices de un par de hojas \(\{i, j\}\), digamos \(i\), mientras que el paquete grande se envía al otro, es decir, \(j\).

El paquete grande induce una mayor demora en la llegada del segundo paquete pequeño a su destino, en comparación con la del primero, y que la diferencia en las demoras de los dos paquetes pequeños variará de manera que refleje cómo los caminos se comparten desde el origen hasta \(i\) y \(j\).

M. Coates, R. Castro, R. Nowak, M. Gadhiok, R. King, and Y. Tsang, “Maximum likelihood network topology identification from edge-based unicast measurements,” Proceedings of the 2002 ACM SIGMETRICS International Conference on Measurement and Modeling of Computer Systems, pp. 11–20, 2002.

# data
suppressMessages(suppressWarnings(library(sand)))
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

dim(delaydata)

## [1] 9567    3

# host machines
host.locs

##  [1] "IST"    "IT"     "UCBkly" "MSU1"   "MSU2"   "UIUC"   "UW1"    "UW2"   
##  [9] "Rice1"  "Rice2"

Diferencias medias de retrzo, simetrizadas para eliminar variaciones dentro de cada par debidas a la recepción de paquetes pequeños frente a paquetes grandes.

meanmat <- with(data = delaydata, expr = by(data = DelayDiff, INDICES = list(SmallPktDest, BigPktDest), FUN = mean))
dim(meanmat)

## [1] 10 10

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)

Los métodos de agrupamiento jerárquico requieren cierta noción de similitud que se construye a partir de los datos observados.

En el contexto del problema de inferencia tomográfica, es lógico utilizar las \(n\) observaciones de las variables aleatorias \(\{X_1,\ldots, X_{N_l}\}\) en las hojas para derivar una matriz \(N_l\times N_l\) de similitudes.

# distancias euclidianas
SSDelayDiff <- with(data = delaydata, expr = by(data = DelayDiff^2, INDICES = list(SmallPktDest, BigPktDest), FUN = sum))
dim(SSDelayDiff)

## [1] 10 10

# matrix de distancias (similitudes)
x <- as.dist(m = 1 / sqrt(SSDelayDiff)) # inversa de la distancia entre retrasos como matriz de similitud
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

Hay bastantes similaridades con la red observada orifinalmente:

# agrupamiento
myclust <- hclust(x, method = "average")
# grafico
plot(myclust, labels=host.locs, axes=FALSE, ylab=NULL, ann=FALSE)

Inferencia de topología de red

Juan Sosa PhD

jcsosam@unal.edu.co