1 Contexto y descripción de los datos

1.1 Fuente y contexto

Los datos utilizados en este análisis provienen de OpenFlights, una base de datos que recopila información geográfica y operacional sobre infraestructura y conexiones de transporte a nivel mundial. Según la descripción de la base, en enero de 2017 el conjunto de datos contenía más de 10.000 aeropuertos, estaciones de tren y terminales de ferry distribuidos alrededor del mundo. La base de rutas dejó de actualizarse en junio de 2014, por lo que los resultados describen la red aérea de ese momento y no la actual.

Para este trabajo se utilizan principalmente dos archivos de OpenFlights: la base de aeropuertos (airports.dat) y la base de rutas aéreas (routes.dat). La primera proporciona las características y la localización geográfica de los aeropuertos, mientras que la segunda permite identificar las conexiones existentes entre ellos.

Desde la perspectiva del análisis de redes, los datos pueden interpretarse de la siguiente manera:

  • Los aeropuertos representan los vértices (nodos) de la red.
  • Las rutas aéreas representan las aristas (conexiones) entre los vértices.
  • Las coordenadas geográficas permiten incorporar una dimensión espacial al análisis de la red.
  • La información sobre país y continente permite estudiar la distribución geográfica de los vértices y sus patrones de conectividad.

Aunque la base original contempla diferentes tipos de infraestructura de transporte, para la construcción de la red analizada se consideran únicamente los registros correspondientes a aeropuertos. De esta manera, la red representa la estructura de conectividad del transporte aéreo a escala mundial.

1.2 Información de los aeropuertos

Cada registro de la base de aeropuertos corresponde a una infraestructura aeroportuaria e incluye información sobre su identificación, localización geográfica y características básicas. Entre las variables utilizadas en este trabajo se encuentran el identificador único del aeropuerto, su nombre, ciudad, país, códigos IATA e ICAO y sus coordenadas geográficas.

Variable Descripción
Airport ID Identificador único asignado por OpenFlights al aeropuerto.
Nombre Nombre del aeropuerto.
City Ciudad principal atendida por el aeropuerto.
País País o territorio donde se encuentra ubicado el aeropuerto.
IATA Código de tres letras asignado por la IATA, cuando está disponible.
ICAO Código de cuatro letras asignado por la ICAO, cuando está disponible.
Latitud Coordenada geográfica de latitud expresada en grados decimales.
Longitud Coordenada geográfica de longitud expresada en grados decimales.
Altitud Elevación del aeropuerto, expresada en pies.
Timezone Diferencia horaria respecto al UTC.
DST Información relacionada con el horario de verano.
Tz database timezone Zona horaria en formato de la base de datos tz (Olson).
Tipo Tipo de infraestructura registrada. En este análisis se consideran únicamente los registros clasificados como airport.
Source Fuente de la información utilizada para el registro.

1.3 Interpretación de los datos como una red

La estructura de los datos permite transformar la información aeroportuaria en una representación de red. En esta representación, cada aeropuerto constituye un vértice, mientras que una ruta entre un aeropuerto de origen y uno de destino constituye una arista.

Inicialmente, las rutas pueden representarse mediante una red dirigida, dado que cada registro distingue explícitamente entre un aeropuerto de origen y uno de destino.

2 Construcción del grafo y control de calidad

# ==============================================================================
# AUDITORÍA Y CONTROL DE CALIDAD
# ==============================================================================
invalid_ids <- airports[is.na(airport_id) | airport_id %in% c("NA", "NULL", "", "\\N")]
invalid_coords <- airports[
  is.na(latitude) | is.na(longitude) |
    latitude < -90 | latitude > 90 |
    longitude < -180 | longitude > 180
]
invalid_routes_na <- routes[
  is.na(source_id) | is.na(destination_id) |
    source_id %in% c("NA", "NULL", "", "\\N") |
    destination_id %in% c("NA", "NULL", "", "\\N")
]
self_loops <- routes[source_id == destination_id & !is.na(source_id)]

# Limpieza
airports[, airport_id := trimws(as.character(airport_id))]
routes[, source_id := trimws(as.character(source_id))]
routes[, destination_id := trimws(as.character(destination_id))]

airports_clean <- airports[
  !is.na(airport_id) & !(airport_id %in% c("NA", "NULL", "", "\\N")) &
    !is.na(latitude) & !is.na(longitude) &
    latitude >= -90 & latitude <= 90 &
    longitude >= -180 & longitude <= 180
]
airports_clean <- unique(airports_clean, by = "airport_id")
setcolorder(airports_clean, "airport_id")

valid_src <- routes$source_id %in% airports_clean$airport_id
valid_dst <- routes$destination_id %in% airports_clean$airport_id
routes_clean <- routes[valid_src & valid_dst]

# Distancias Haversine
coords_src <- airports_clean[match(routes_clean$source_id, airport_id), .(longitude, latitude)]
coords_dst <- airports_clean[match(routes_clean$destination_id, airport_id), .(longitude, latitude)]
routes_clean[, distance_km := distHaversine(coords_src, coords_dst) / 1000]

# Data frames finales para igraph
edges_df <- data.frame(
  from = routes_clean$source_id,
  to   = routes_clean$destination_id,
  routes_clean[, !(names(routes_clean) %in% c("source_id", "destination_id")), with = FALSE],
  stringsAsFactors = FALSE
)
vertices_df <- as.data.frame(airports_clean, stringsAsFactors = FALSE)

Antes de construir la red se realizó un proceso de auditoría y depuración de los datos, con el objetivo de verificar la validez de los aeropuertos y las rutas utilizadas en el análisis. En particular, se revisaron los identificadores de los aeropuertos, sus coordenadas geográficas y la correspondencia entre los aeropuertos y las rutas registradas.

# 1. Red dirigida inicial
g_directed <- graph_from_data_frame(d = edges_df, vertices = vertices_df, directed = TRUE)

# 2. Eliminar lazos (manteniendo por ahora las aristas múltiples)
g_directed <- igraph::simplify(g_directed, remove.loops = TRUE, remove.multiple = FALSE)

# 3. Atributos para la agregación
E(g_directed)$registros       <- 1
E(g_directed)$aerolineas_list <- E(g_directed)$airline

# 4. Colapsar a red no dirigida
g_undirected <- as.undirected(
  g_directed,
  mode = "collapse",
  edge.attr.comb = list(
    registros       = "sum",
    aerolineas_list = function(x) length(unique(x)),
    distance_km     = "mean",
    "ignore"
  )
)
g_undirected <- set_edge_attr(g_undirected, "num_registros", value = E(g_undirected)$registros)
g_undirected <- set_edge_attr(g_undirected, "num_aerolineas", value = E(g_undirected)$aerolineas_list)
g_undirected <- delete_edge_attr(g_undirected, "registros")
g_undirected <- delete_edge_attr(g_undirected, "aerolineas_list")

# Peso
E(g_undirected)$weight <- E(g_undirected)$num_registros

La auditoría muestra que los 7.698 aeropuertos considerados como vértices presentan identificadores y coordenadas válidas, por lo que no fue necesario excluir ningún aeropuerto por estas razones. En cuanto a las rutas, se excluyeron 892 registros debido a la presencia de aeropuertos sin correspondencia en el conjunto de vértices o a valores faltantes.

Durante este proceso también se identificó 1 lazo (self-loops), correspondientes a rutas cuyo aeropuerto de origen y destino coinciden. Una vez realizada la depuración, se obtuvieron 66.771 aristas dirigidas.

2.1 Grafo no dirigido

Posteriormente se construyó la red principal no dirigida, en la que dos aeropuertos están conectados si existe al menos una ruta entre ellos, colapsando los duplicados.

Esta representación permite estudiar propiedades estructurales de la red mundial de transporte aéreo, como su conectividad, distribución geográfica, centralidad, componentes y vulnerabilidad frente a la pérdida de determinados vértices o conexiones.

La red final no dirigida mantiene su orden de 7.698 aeropuertos, mientras que su tamaño es de 18.858 aristas.

Como peso de las aristas (\(w_{ij}\)) se utilizó el número de registros de vuelo, que actúa como una aproximación a la densidad de tráfico, el volumen de operaciones y la carga logística sostenida por la infraestructura aeroportuaria en el par origen-destino. En el cálculo de la fuerza del nodo (\(s_i = \sum_{j} w_{ij}\)), este peso opera como una medida de intensidad operativa, reflejando la oferta total de vuelos acumulada por cada aeropuerto dentro de la red.

Códigos de resultados de la auditoría y control de calidad
cat("1. Aeropuertos excluidos por ID inválido:", nrow(invalid_ids), "\n",
    "2. Aeropuertos excluidos por coordenadas:", nrow(invalid_coords), "\n",
    "3. Vértices válidos finales (orden):", nrow(vertices_df), "\n",
    "4. Rutas excluidas (huérfanas/NAs):", nrow(routes) - nrow(routes_clean), "\n",
    "5. Lazos detectados:", nrow(self_loops), "\n",
    "6. Aristas dirigidas:", nrow(edges_df), "\n",
    "7. Número de aristas (tamaño) red no dirigida:", ecount(g_undirected))
## 1. Aeropuertos excluidos por ID inválido: 0 
##  2. Aeropuertos excluidos por coordenadas: 0 
##  3. Vértices válidos finales (orden): 7698 
##  4. Rutas excluidas (huérfanas/NAs): 892 
##  5. Lazos detectados: 1 
##  6. Aristas dirigidas: 66771 
##  7. Número de aristas (tamaño) red no dirigida: 18858

Los principales resultados del proceso de control de calidad se resumen en la siguiente tabla:

Control realizado Resultado
Aeropuertos excluidos por ID inválido 0
Aeropuertos excluidos por coordenadas inválidas 0
Vértices válidos finales 7.698
Rutas excluidas por aeropuertos huérfanos o valores faltantes 892
Lazos detectados 1
Aristas dirigidas 66.771
Aristas no dirigidas finales 18.858

3 Descripción y visualización

Se calculan algunas medidas útiles para caracterizar la red, como el orden, el tamaño, la densidad y el grado promedio. Dado que la red es muy grande, los cálculos se realizan también para su componente gigante (la componente conexa de mayor tamaño), como se muestra en la siguiente tabla.

cl <- components(g_undirected)

num_componentes <- cl$no
ordenes         <- cl$csize
orden_gigante   <- max(ordenes)
total_nodos     <- vcount(g_undirected)

prop_gigante <- orden_gigante / total_nodos
pct_gigante  <- prop_gigante * 100

largest_cluster_id <- which.max(cl$csize)
giant_component    <- induced_subgraph(g_undirected, which(cl$membership == largest_cluster_id))

total_relaciones <- ecount(g_undirected)
pct_tamgigante   <- 100 * ecount(giant_component) / total_relaciones
Métricas descriptivas comparativas entre la red completa y la componente gigante.
Métrica Red completa Componente gigante
Orden (número de vértices) 7.698 3.188
Tamaño (número de aristas) 18.858 18.833
Densidad 0,001 0,004
Grado promedio 4,899 11,815
Número total de componentes 4.491 1 (gigante)
Orden relativo (% de la red) 100,00 % 41,41 %
Tamaño relativo (% de la red) 100,00 % 99,87 %

La componente gigante está conformada por 3.188 aeropuertos, que representan aproximadamente el 41,4 % de los 7.698 aeropuertos de la red completa. Sin embargo, esta componente contiene 18.833 de las 18.858 aristas de la red, equivalentes a aproximadamente el 99,9 % de las conexiones.

Este resultado evidencia una marcada concentración de las conexiones en una única estructura conectada. Aunque una proporción importante de los aeropuertos no pertenece a la componente gigante, estos concentran muy pocas conexiones en comparación con los que forman parte de ella.

La densidad de 0,004, equivalente al 0,4 %, muestra que la componente gigante es una red dispersa: solo una pequeña proporción de las conexiones potenciales entre sus aeropuertos está presente como conexión directa. Esto no implica una baja conectividad en términos de alcanzabilidad, ya que los aeropuertos de la componente gigante están conectados entre sí mediante caminos, aunque generalmente requieren pasar por otros aeropuertos para alcanzar destinos sin conexión directa.

Por su parte, el grado promedio de 11,815 indica que cada aeropuerto de la componente gigante tiene, en promedio, aproximadamente 12 conexiones directas. Este valor debe interpretarse conjuntamente con la baja densidad: la red posee pocas conexiones directas frente al total posible, pero estas son suficientes para mantener una única estructura conectada que concentra prácticamente todas las aristas.

3.1 Distribución del grado

# Ajuste de una ley de potencias discreta (poweRlaw)
m_pl <- displ$new(deg_giant)
m_pl$setXmin(estimate_xmin(m_pl))
alpha  <- m_pl$getPars()
d_min  <- m_pl$getXmin()

A continuación se analiza la distribución del grado de los aeropuertos pertenecientes a la componente gigante. Como se observa en los gráficos, la distribución presenta una fuerte asimetría hacia la derecha: la mayoría de los aeropuertos tiene pocas conexiones directas, mientras que un grupo reducido concentra un número elevado de conexiones.

Para caracterizar el comportamiento de la cola se ajustó una ley de potencias discreta mediante el paquete poweRlaw. El valor de \(d_{\min}\) se seleccionó minimizando la distancia de Kolmogorov-Smirnov entre la distribución empírica y la ajustada, obteniendo

\[ \hat{d}_{\min}=2 \qquad\text{y}\qquad \hat{\alpha}=1,71. \]

Por tanto, el ajuste considera los aeropuertos con grado \(k\geq 2\), que corresponden al 77,9 % de la distribución analizada. El valor de \(\hat{\alpha}\) determina la forma de decrecimiento de la frecuencia de los grados a medida que estos aumentan. La línea azul representa el ajuste teórico y la línea vertical punteada indica el valor estimado de \(d_{\min}\).

En el panel log-log se aprecia que, a mayor grado, mayor es la dispersión relativa de las frecuencias: el grado presenta un comportamiento heterocedástico, cuya variabilidad aumenta con el grado, lo que indica una alta dispersión entre los aeropuertos de más conexiones.

df_grados <- data.table(grado = deg_giant)

d_tail  <- d_min:max(deg_giant)
p_tail  <- sum(deg_giant >= d_min) / length(deg_giant)
n_total <- length(deg_giant)
df_tail_fit <- data.table(
  grado    = d_tail,
  N_fitted = p_tail * dist_pdf(m_pl, d_tail) * n_total
)

df_deg_counts <- df_grados[, .(N = .N), by = grado][grado > 0]

p_deg_lin <- ggplot(df_grados, aes(x = grado)) +
  geom_histogram(binwidth = 2, fill = "#2b5c8f", color = "white", alpha = 0.85) +
  geom_line(data = df_tail_fit, aes(x = grado, y = N_fitted), color = "royalblue", linewidth = 1) +
  geom_vline(xintercept = d_min, linetype = "dashed", color = "firebrick", linewidth = 0.8) +
  coord_cartesian(xlim = c(0, 50)) +
  labs(title = "A. Escala lineal",
       x = "Grado k (N° de destinos directos)", y = "Número de aeropuertos") +
  theme_minimal(base_size = 11) +
  theme(panel.grid.minor = element_blank(), plot.title = element_text(face = "bold"))

p_deg_log <- ggplot(df_deg_counts, aes(x = grado, y = N)) +
  geom_point(color = "#d9381e", alpha = 0.7, size = 2.5) +
  geom_line(data = df_tail_fit, aes(x = grado, y = N_fitted), color = "royalblue", linewidth = 1.1) +
  geom_vline(xintercept = d_min, linetype = "dashed", color = "firebrick", linewidth = 0.8) +
  scale_x_log10(labels = trans_format("log10", math_format(10^.x))) +
  scale_y_log10(labels = trans_format("log10", math_format(10^.x))) +
  labs(title = "B. Escala log-log", x = "Grado k (log10)", y = "Frecuencia N(k) (log10)") +
  theme_minimal(base_size = 11) +
  theme(panel.grid.minor = element_blank(), plot.title = element_text(face = "bold"))

p_deg_lin + p_deg_log +
  plot_annotation(
    title = "Distribución del grado (componente gigante)",
    subtitle = paste0("Línea azul: ajuste ley de potencias | Línea roja punteada: k_min = ", d_min),
    theme = theme(
      plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
      plot.subtitle = element_text(size = 10, hjust = 0.5, color = "gray30")
    )
  )

3.2 Distribución de la fuerza

st_giant <- strength(giant_component, weights = E(giant_component)$num_registros)

df_fuerza    <- data.table(fuerza = st_giant)
df_st_counts <- df_fuerza[, .(N = .N), by = fuerza][fuerza > 0]

p_st_lin <- ggplot(df_fuerza, aes(x = fuerza)) +
  geom_histogram(bins = 30, fill = "#1b9e77", color = "white", alpha = 0.85) +
  coord_cartesian(xlim = c(0, quantile(st_giant, 0.95))) +
  labs(title = "A. Escala lineal",
       x = "Fuerza s (suma de registros de vuelo)", y = "Número de aeropuertos") +
  theme_minimal(base_size = 11) +
  theme(panel.grid.minor = element_blank(), plot.title = element_text(face = "bold"))

p_st_log <- ggplot(df_st_counts, aes(x = fuerza, y = N)) +
  geom_point(color = "#d95f02", alpha = 0.7, size = 2.5) +
  geom_smooth(method = "lm", color = "black", linetype = "dashed", se = FALSE, linewidth = 0.7) +
  scale_x_log10(labels = trans_format("log10", math_format(10^.x))) +
  scale_y_log10(labels = trans_format("log10", math_format(10^.x))) +
  labs(title = "B. Escala log-log", x = "Fuerza s (log10)", y = "Frecuencia N(s) (log10)") +
  theme_minimal(base_size = 11) +
  theme(panel.grid.minor = element_blank(), plot.title = element_text(face = "bold"))

p_st_lin + p_st_log +
  plot_annotation(
    title = "Distribución de la fuerza (tráfico) en la componente gigante",
    theme = theme(plot.title = element_text(size = 14, face = "bold", hjust = 0.5))
  )

En conjunto, los resultados muestran una red globalmente conectada pero localmente dispersa, caracterizada por una fuerte concentración de las conexiones en una única componente gigante. Esta estructura resulta particularmente relevante para el análisis posterior de la vulnerabilidad de la red, ya que cambios en determinados aeropuertos pueden afectar la comunicación entre diferentes partes de la red aeroportuaria internacional.

3.3 Longitudes de las rutas

Dado que la distribución de la longitud de las rutas presenta una fuerte asimetría hacia la derecha, con una alta densidad en distancias cortas, se aplicó una transformación logarítmica (scale_x_log10) para expandir la escala visual y estabilizar la varianza.

df_edges_undir <- igraph::as_data_frame(g_undirected, what = "edges")

cont_v1 <- V(g_undirected)$continent[match(df_edges_undir$from, V(g_undirected)$name)]
cont_v2 <- V(g_undirected)$continent[match(df_edges_undir$to,   V(g_undirected)$name)]

df_edges_undir_long <- rbind(
  data.table(distance_km = df_edges_undir$distance_km, continent = cont_v1),
  data.table(distance_km = df_edges_undir$distance_km, continent = cont_v2)
)
df_edges_undir_long <- df_edges_undir_long[!is.na(continent) & continent != "Otros"]

ggplot(df_edges_undir_long, aes(x = distance_km, fill = continent)) +
  geom_density(alpha = 0.5) +
  scale_x_log10(labels = comma) +
  labs(
    title = "Distribución de longitud de rutas (km) por continente",
    subtitle = "Red no dirigida: cada ruta se asigna a ambos continentes",
    x = "Distancia física en km (escala log10)",
    y = "Densidad de la longitud",
    fill = "Continente"
  ) +
  theme_minimal()

Esta transformación resultó fundamental al analizar regiones como Oceanía: en la escala real, la concentración de rutas insulares muy cortas genera un pico de densidad extremadamente alto que puede dar la falsa impresión de que la región domina el volumen global de vuelos. La escala logarítmica permite revelar la verdadera estructura del perfil continental, separando la concentración local en trayectos cortos del volumen total de operaciones.

Como cada curva está normalizada (su área es igual a 1), el gráfico compara la forma de las distribuciones de cada continente y no su volumen de rutas: Oceanía se distingue por una mayor proporción de rutas muy cortas, asociadas a conexiones entre islas, mientras que en los demás continentes la mayor parte de las rutas se concentra entre unos 500 y 3.000 km.

df_edges_dir <- igraph::as_data_frame(g_directed, what = "edges")

df_edges_dir$cont_out <- V(g_directed)$continent[match(df_edges_dir$from, V(g_directed)$name)]
df_edges_dir$cont_in  <- V(g_directed)$continent[match(df_edges_dir$to,   V(g_directed)$name)]

df_dir_long <- rbind(
  data.table(distance_km = df_edges_dir$distance_km, continent = df_edges_dir$cont_out,
             tipo_flujo = "Vuelos de salida (origen)"),
  data.table(distance_km = df_edges_dir$distance_km, continent = df_edges_dir$cont_in,
             tipo_flujo = "Vuelos de entrada (destino)")
)
df_dir_long <- df_dir_long[!is.na(continent) & continent != "Otros"]

ggplot(df_dir_long, aes(x = distance_km, fill = continent)) +
  geom_density(alpha = 0.5) +
  scale_x_log10(labels = comma) +
  facet_wrap(~ tipo_flujo, ncol = 2) +
  labs(
    title = "Densidad de longitud de rutas: vuelos de salida vs. entrada",
    subtitle = "Red dirigida: separada por origen y destino de cada ruta",
    x = "Distancia física en km (escala log10)",
    y = "Densidad de la longitud",
    fill = "Continente"
  ) +
  theme_minimal() +
  theme(strip.text = element_text(size = 11, face = "bold"),
        plot.title = element_text(size = 14, face = "bold"))

Al comparar las distribuciones de longitud de ruta entre los vuelos de salida (origen) y de entrada (destino) en la red dirigida, se observa una simetría prácticamente perfecta. Este comportamiento se explica por la alta reciprocidad del transporte aéreo comercial: la oferta de rutas en un sentido (\(A \to B\)) se corresponde de forma bilateral con el trayecto de retorno (\(B \to A\)), manteniendo una distribución de distancias equivalente entre los flujos de entrada y salida de cada continente.

3.4 Visualizaciones

Primero se presenta una visualización geográfica que permite ubicar en el espacio los aeropuertos de mayor grado.

world_map <- map_data("world") %>% filter(region != "Antarctica")

giant_tidy <- as_tbl_graph(giant_component) %>%
  mutate(
    deg      = centrality_degree(),
    deg_norm = (deg - min(deg)) / (max(deg) - min(deg))
  )

# Nota: se usa geom_edge_link en lugar de geom_edge_arc porque dibujar ~19.000
# arcos es mucho más lento. Para la versión final puede volver a geom_edge_arc.
ggraph(giant_tidy, layout = "manual",
       x = V(giant_component)$longitude, y = V(giant_component)$latitude) +
  geom_polygon(data = world_map, aes(x = long, y = lat, group = group),
               fill = "#f2f2f2", color = "#cccccc", linewidth = 0.2) +
  geom_edge_link(color = "#2b5c8f", alpha = 0.06, width = 0.2) +
  geom_node_point(aes(size = deg, color = continent, alpha = deg), shape = 16) +
  scale_size_continuous(range = c(0.4, 2.8), name = "Grado (conexiones)") +
  scale_alpha_continuous(range = c(0.35, 0.95), guide = "none") +
  scale_color_brewer(palette = "Set1", name = "Continente") +
  coord_fixed(1.3) +
  theme_void() +
  labs(
    title = "Red geográfica global — componente gigante",
    subtitle = "Ubicación espacial (latitud/longitud) con degradado por intensidad de conectividad"
  ) +
  theme(legend.position = "bottom", legend.box = "horizontal",
        plot.title = element_text(face = "bold", size = 14),
        plot.subtitle = element_text(color = "grey30"))

Para evitar el colapso visual (hairball effect) de la visualización topológica de los 3.188 nodos, se aplicó un filtro basado en el criterio de atipicidad de Tukey (\(k \ge Q_3 + 1{,}5 \times \text{IQR}\)) sobre la distribución del grado.

q3         <- quantile(deg_giant, 0.75)
iqr        <- IQR(deg_giant)
umbral_iqr <- as.numeric(q3 + 1.5 * iqr)

giant_tidy_sub <- giant_tidy %>%
  mutate(
    deg          = centrality_degree(),
    strength_val = centrality_degree(weights = weight)
  ) %>%
  filter(deg >= umbral_iqr) %>%
  mutate(label_hub = if_else(rank(-strength_val) <= 10, iata, NA_character_))

set.seed(123)
ggraph(giant_tidy_sub, layout = "fr", weights = weight) +
  geom_edge_link(aes(alpha = weight, width = weight), color = "#999999") +
  scale_edge_width(range = c(0.15, 0.9), name = "Peso ruta") +
  scale_edge_alpha(range = c(0.08, 0.35), guide = "none") +
  geom_node_point(aes(size = strength_val, color = continent), alpha = 0.88) +
  scale_size_continuous(range = c(1.8, 6.5), name = "Fuerza (strength)") +
  scale_color_brewer(palette = "Set1", name = "Continente") +
  geom_node_text(aes(label = label_hub), repel = TRUE, size = 3, fontface = "bold",
                 color = "#222222", max.overlaps = 15) +
  theme_void() +
  labs(
    title = "Red topológica global — componente gigante",
    subtitle = paste0("Subgrafo filtrado por criterio de Tukey (k ≥ ", round(umbral_iqr, 1),
                      "). Layout de Fruchterman-Reingold con jerarquía por fuerza.")
  ) +
  theme(legend.position = "bottom", legend.box = "horizontal",
        plot.title = element_text(face = "bold", size = 13),
        plot.subtitle = element_text(color = "grey30", size = 10),
        legend.title = element_text(face = "bold", size = 9),
        legend.text = element_text(size = 8)) +
  guides(color = guide_legend(override.aes = list(size = 4)),
         size = guide_legend(title.position = "top"))

Al contrastar la representación geográfica con la topológica (layout Fruchterman-Reingold ponderado por la fuerza del nodo), se evidencia un fenómeno clave: para los nodos de grado alto, la disposición libre por fuerzas de atracción agrupa a América, Europa y Asia en clústeres bien definidos y casi autónomos, imitando una lógica de proximidad espacial a pesar de no usar coordenadas. En contraste, Oceanía se muestra dispersa entre América y Asia, mientras que África aparece fragmentada y conectada directamente hacia Europa y Asia. Este patrón refleja cómo la eficiencia en tiempos de vuelo y consumo de combustible fomenta la integración regional directa entre los bloques principales, dejando a regiones como Oceanía y África en una posición de dependencia funcional supeditada a la cercanía geográfica con los hubs continentales dominantes.

nombres    <- V(giant_component)$name
continente <- V(giant_component)$continent
grados     <- deg_giant

indices_validos <- which(!is.na(continente) & continente != "Otros")
orden_indices <- indices_validos[order(continente[indices_validos], -grados[indices_validos])]

adj_mat      <- as.matrix(as_adjacency_matrix(giant_component, sparse = FALSE))
adj_ordenada <- adj_mat[orden_indices, orden_indices]
rm(adj_mat); invisible(gc())

conts_ordenados <- continente[orden_indices]
conts_unicos    <- unique(conts_ordenados)
cont_code <- as.numeric(factor(conts_ordenados, levels = conts_unicos))

misma_comunidad <- outer(conts_ordenados, conts_ordenados, FUN = "==")

adj_color <- matrix(0, nrow = nrow(adj_ordenada), ncol = ncol(adj_ordenada))
adj_color[adj_ordenada == 1 & misma_comunidad] <-
  cont_code[row(adj_color)[adj_ordenada == 1 & misma_comunidad]]
codigo_inter <- length(conts_unicos) + 1
adj_color[adj_ordenada == 1 & !misma_comunidad] <- codigo_inter

paleta_set1    <- brewer.pal(n = max(length(conts_unicos), 3), name = "Set1")[1:length(conts_unicos)]
color_inter    <- "grey20"
colores_matriz <- c("white", paleta_set1, color_inter)

image(1:nrow(adj_color), 1:ncol(adj_color), adj_color,
      col = colores_matriz, useRaster = TRUE,
      xlab = "Aeropuertos (ordenados por continente y grado ↓)",
      ylab = "Aeropuertos (ordenados por continente y grado ↓)",
      main = "Matriz de adyacencia — relaciones intra e intercontinentales")
legend("topleft", legend = c(conts_unicos, "Intercontinental"),
       fill = c(paleta_set1, color_inter), bg = "white", cex = 0.8,
       title = "Tipo de relación")

A partir del ordenamiento por continente y grado descendente, la matriz de adyacencia (renderizada como mapa de bits con useRaster por sus dimensiones) muestra un claro patrón de bloques densos a lo largo de la diagonal principal. Esta concentración confirma una marcada segregación espacial y una alta densidad de vuelos intracontinentales. Los puntos fuera de la diagonal son notablemente más dispersos y corresponden a las rutas intercontinentales directas, sostenidas casi exclusivamente por los aeropuertos de mayor grado, ubicados en la cúspide de cada bloque.

3.5 Redundancia operacional

Para evaluar qué tanto depende la red de rutas operadas por una sola aerolínea, se compara la componente gigante con la subred que conserva únicamente las rutas operadas por dos o más aerolíneas.

g_redundante <- subgraph.edges(
  giant_component,
  eids = E(giant_component)[E(giant_component)$num_aerolineas >= 2],
  delete.vertices = FALSE
)

comp_full <- components(giant_component)
comp_red  <- components(g_redundante)

m_full <- c(nodos = vcount(giant_component), aristas = ecount(giant_component),
            dens = edge_density(giant_component), comps = comp_full$no,
            gc = max(comp_full$csize))
m_red  <- c(nodos = vcount(g_redundante), aristas = ecount(g_redundante),
            dens = edge_density(g_redundante), comps = comp_red$no,
            gc = max(comp_red$csize))

fmt_int <- function(x) format(x, big.mark = ".", decimal.mark = ",", scientific = FALSE, trim = TRUE)
fmt_dec <- function(x, d = 4) format(round(x, d), nsmall = d, decimal.mark = ",", scientific = FALSE, trim = TRUE)
fmt_var <- function(new, old) {
  v <- (new - old) / old * 100
  paste0(ifelse(v > 0, "+", ""), format(round(v, 1), decimal.mark = ",", nsmall = 1), " %")
}

tabla_comparativa <- data.frame(
  Métrica = c("Total de nodos (aeropuertos)", "Total de aristas (rutas activas)",
              "Densidad de la red", "Número de componentes",
              "Tamaño de la componente gigante (nodos)"),
  `Red completa (≥ 1 aerolínea)` = c(fmt_int(m_full["nodos"]), fmt_int(m_full["aristas"]),
                                     fmt_dec(m_full["dens"]), fmt_int(m_full["comps"]),
                                     fmt_int(m_full["gc"])),
  `Red redundante (≥ 2 aerolíneas)` = c(fmt_int(m_red["nodos"]), fmt_int(m_red["aristas"]),
                                        fmt_dec(m_red["dens"]), fmt_int(m_red["comps"]),
                                        fmt_int(m_red["gc"])),
  Variación = c(fmt_var(m_red["nodos"], m_full["nodos"]),
                fmt_var(m_red["aristas"], m_full["aristas"]),
                fmt_var(m_red["dens"], m_full["dens"]),
                paste0("+", fmt_int(m_red["comps"] - m_full["comps"])),
                fmt_var(m_red["gc"], m_full["gc"])),
  check.names = FALSE
)

kable(tabla_comparativa, align = c("l", "r", "r", "r"), row.names = FALSE,
      caption = "Evaluación de redundancia operacional: red global vs. red sin monopolios de ruta (≥ 2 aerolíneas)")
Evaluación de redundancia operacional: red global vs. red sin monopolios de ruta (≥ 2 aerolíneas)
Métrica Red completa (≥ 1 aerolínea) Red redundante (≥ 2 aerolíneas) Variación
Total de nodos (aeropuertos) 3.188 3.188 0,0 %
Total de aristas (rutas activas) 18.833 8.362 -55,6 %
Densidad de la red 0,0037 0,0016 -55,6 %
Número de componentes 1 1.353 +1.352
Tamaño de la componente gigante (nodos) 3.188 1.834 -42,5 %
set.seed(123)
layout_fijo <- create_layout(as_tbl_graph(giant_component), layout = "fr")

p_orig <- ggraph(layout_fijo) +
  geom_edge_link(alpha = 0.1, color = "#2b5c8f") +
  geom_node_point(aes(color = continent), size = 1.5, alpha = 0.7) +
  scale_color_brewer(palette = "Set1", name = "Continente") +
  theme_void() +
  labs(title = "Panel A: red global completa (≥ 1 aerolínea)",
       subtitle = paste0("Rutas: ", fmt_int(ecount(giant_component)),
                         " | Incluye monopolios de ruta y conexiones únicas"))

p_filt <- ggraph(as_tbl_graph(g_redundante), layout = "manual",
                 x = layout_fijo$x, y = layout_fijo$y) +
  geom_edge_link(alpha = 0.2, color = "#e41a1c") +
  geom_node_point(aes(color = continent), size = 1.5, alpha = 0.7) +
  scale_color_brewer(palette = "Set1", name = "Continente") +
  theme_void() +
  labs(title = "Panel B: red con redundancia operacional (≥ 2 aerolíneas)",
       subtitle = paste0("Rutas: ", fmt_int(ecount(g_redundante))))

(p_orig + p_filt) + plot_layout(guides = "collect") & theme(legend.position = "bottom")

Al conservar solo las rutas operadas por al menos dos aerolíneas, el número de rutas pasa de 18.833 a 8.362 (-55,6 %), y la red se divide en 1.353 componentes, cuya componente mayor reúne 1.834 de los 3.188 aeropuertos (57,5 %). Es decir, una parte considerable de la conectividad de la red depende de rutas operadas por una única aerolínea, lo que constituye una fuente de fragilidad operacional.

4 Distancia geográfica y distancia en la red

4.1 Marco teórico y métricas

Para analizar la eficiencia del transporte aéreo en la componente gigante, se comparan tres métricas para cada par de aeropuertos \((i,j)\):

Medida Definición Naturaleza
\(d_G(i,j)\) Número mínimo de aristas (vuelos o tramos). Topológica (sin peso)
\(d_E(i,j)\) Distancia geográfica en línea recta sobre la superficie terrestre (km). Geodésica (Haversine)
\(d_{\text{red}}(i,j)\) Distancia mínima ponderada recorrida sobre la red (km). Ruta óptima en red

4.1.1 Distancia de gran círculo (\(d_E\)) e índice de rodeo (\(S_{ij}\))

La distancia de gran círculo \(d_E(i,j)\) es el trayecto directo más corto sobre la esfera terrestre entre las coordenadas de latitud (\(\varphi\)) y longitud (\(\lambda\)) de dos aeropuertos, calculada con la fórmula de Haversine:

\[d_E(i,j) = 2R \arcsin \sqrt{\sin^2\!\left(\frac{\varphi_j-\varphi_i}{2}\right) + \cos\varphi_i\,\cos\varphi_j\,\sin^2\!\left(\frac{\lambda_j-\lambda_i}{2}\right)}\]

donde \(R = 6\,378{,}137\text{ km}\) es el radio terrestre usado por geosphere. La distancia ponderada sobre la red se define como

\[d_{\text{red}}(i,j) = \min_{\gamma:\, i \to j} \sum_{(u,v)\in\gamma} \ell_{uv},\]

donde \(\ell_{uv}\) es la longitud en kilómetros de la arista \((u,v)\). Por la desigualdad triangular esférica, \(d_{\text{red}}(i,j) \ge d_E(i,j)\).

Para evaluar la ineficiencia espacial se define el índice de rodeo (\(S_{ij}\)) y su sobrecosto porcentual:

\[S_{ij} = \frac{d_{\text{red}}(i,j)}{d_E(i,j)} \ge 1 \quad \implies \quad \text{Sobrecosto (\%)} = (S_{ij} - 1) \times 100\]


4.2 Muestreo probabilístico

Dado el costo computacional de evaluar todos los pares, la eficiencia espacial se evaluó sobre un muestreo aleatorio simple (MAS) de pares de aeropuertos, con tamaño

\[n = \frac{Z^2 \sigma^2}{e^2},\]

donde \(Z = 1{,}96\) (confianza del 95 %), \(e = 100\text{ km}\) es el error absoluto máximo tolerable y \(\sigma^2\) la varianza empírica de la longitud de las rutas (distance_km). Como esta varianza se calcula sobre las aristas y no sobre los pares de aeropuertos, que en general están más alejados entre sí, se fija además un mínimo de 5.000 pares.

km <- E(giant_component)$distance_km
sigma_dist <- sd(km, na.rm = TRUE)
media_dist <- mean(km, na.rm = TRUE)

Z_alpha2  <- 1.96
error_abs <- 100

n_calculado <- ceiling((Z_alpha2^2 * sigma_dist^2) / error_abs^2)
N_MUESTRA   <- max(n_calculado, 5000)

kable(
  data.frame(
    Parámetro = c("Desviación estándar de distance_km (km)", "Media de distance_km (km)",
                  "Error tolerable (km)", "Error relativo (% de la media)",
                  "Nivel de confianza", "n calculado (pares)", "n final ajustado (pares)"),
    Valor = c(fmt_num(sigma_dist, 2), fmt_num(media_dist, 2), error_abs,
              fmt_num(100 * error_abs / media_dist, 1), "95 % (Z = 1,96)",
              fmt_num(n_calculado, 0), fmt_num(N_MUESTRA, 0))
  ),
  caption = "Justificación y diseño del tamaño muestral"
)
Justificación y diseño del tamaño muestral
Parámetro Valor
Desviación estándar de distance_km (km) 1.962,84
Media de distance_km (km) 1.771,77
Error tolerable (km) 100
Error relativo (% de la media) 5,6
Nivel de confianza 95 % (Z = 1,96)
n calculado (pares) 1.481
n final ajustado (pares) 5.000

4.3 Propiedades topológicas globales (\(d_G\))

# distance_table() devuelve cuántos pares hay a cada distancia (1, 2, 3, ...)
# sin construir la matriz completa de distancias de n x n.
dt_G <- distance_table(giant_component, directed = FALSE)$res

L_prom   <- sum(seq_along(dt_G) * dt_G) / sum(dt_G)
diametro <- length(dt_G)
p90      <- which(cumsum(dt_G) / sum(dt_G) >= 0.90)[1]

kable(
  data.frame(
    Métrica = c("Distancia promedio (L)", "Diámetro topológico (D)",
                "Diámetro efectivo (percentil 90)"),
    Valor   = c(fmt_num(L_prom, 3), diametro, p90)
  ),
  caption = "Métricas de distancia geodésica en la componente gigante"
)
Métricas de distancia geodésica en la componente gigante
Métrica Valor
Distancia promedio (L) 3,959
Diámetro topológico (D) 12
Diámetro efectivo (percentil 90) 5

La distancia geodésica promedio (\(L \approx 3,96\)) y el diámetro efectivo (5) reflejan una estructura de mundo pequeño (small-world): el 90 % de los pares de aeropuertos se conecta con 5 vuelos o menos. La presencia de hubs centralizadores permite conectar la gran mayoría de los pares con muy pocas escalas.

ggplot(data.frame(d_G = factor(seq_along(dt_G)), n = dt_G), aes(x = d_G, y = n)) +
  geom_col(fill = "#2b5c8f", color = "white", alpha = 0.85) +
  scale_y_continuous(labels = label_number(big.mark = ".", decimal.mark = ",")) +
  theme_minimal() +
  theme(panel.grid.minor = element_blank()) +
  labs(title = "Distribución de distancias geodésicas (d_G)",
       x = "Número de escalas (aristas)", y = "Frecuencia de pares")


4.4 Ineficiencia espacial e índice de rodeo (\(S_{ij}\))

set.seed(42)

n_nodos <- vcount(giant_component)
coords  <- cbind(V(giant_component)$longitude, V(giant_component)$latitude)

# Muestra de pares sin pares repetidos ni pares (i, i)
N_SOBRE <- ceiling(N_MUESTRA * 1.3)
pares <- data.frame(
  i = sample(n_nodos, N_SOBRE, replace = TRUE),
  j = sample(n_nodos, N_SOBRE, replace = TRUE)
) %>%
  filter(i != j) %>%
  mutate(a = pmin(i, j), b = pmax(i, j)) %>%
  distinct(a, b) %>%
  slice_head(n = N_MUESTRA) %>%
  transmute(i = a, j = b)

# Distancias calculadas solo desde los orígenes de la muestra (no la matriz completa)
src  <- sort(unique(pares$i))
pares$fila <- match(pares$i, src)

dist_G_s  <- distances(giant_component, v = src, mode = "all", weights = NA)
dist_km_s <- distances(giant_component, v = src, mode = "all", weights = km)

df_muestra <- pares %>%
  mutate(
    d_G   = dist_G_s[cbind(fila, j)],
    d_red = dist_km_s[cbind(fila, j)],
    d_E   = distHaversine(coords[i, ], coords[j, ]) / 1000
  ) %>%
  filter(is.finite(d_G), d_E > 10) %>%      # se excluyen pares a menos de 10 km
  mutate(S_ij = d_red / d_E)

4.4.1 Asociación topológica vs. geográfica

test_spearman <- cor.test(df_muestra$d_G, df_muestra$d_E, method = "spearman", exact = FALSE)

kable(
  data.frame(
    `Variable 1` = "Distancia geodésica (d_G)",
    `Variable 2` = "Distancia gran círculo (d_E)",
    `Coeficiente (rho)` = fmt_num(test_spearman$estimate, 3),
    `p-valor` = ifelse(test_spearman$p.value < 0.005, "< 0,005", fmt_num(test_spearman$p.value, 4)),
    `Muestra (n)` = fmt_num(nrow(df_muestra), 0),
    check.names = FALSE
  ),
  align = "c",
  row.names = FALSE,
  caption = "Correlación de rangos de Spearman entre escalas topológicas y distancia física"
)
Correlación de rangos de Spearman entre escalas topológicas y distancia física
Variable 1 Variable 2 Coeficiente (rho) p-valor Muestra (n)
Distancia geodésica (d_G) Distancia gran círculo (d_E) 0,373 < 0,005 5.000
ggplot(df_muestra, aes(x = factor(d_G), y = d_E)) +
  geom_boxplot(fill = "#2b5c8f", alpha = 0.6, outlier.alpha = 0.3) +
  theme_minimal() +
  labs(title = "Distancia de gran círculo según número de escalas",
       x = "Escalas topológicas (d_G)", y = "Distancia geográfica d_E (km)")

La correlación de Spearman (\(\rho = 0,373\)) refleja una relación débil entre el número de escalas y la distancia física. Esta desconexión es característica de las redes aeroportuarias organizadas en hubs: dos destinos distantes pueden requerir pocas escalas si ambos conectan con centros neurálgicos, mientras que trayectos regionales cortos pueden requerir múltiples conexiones.

4.4.2 Distribución del índice \(S_{ij}\)

med_S <- median(df_muestra$S_ij)

resumen_S <- data.frame(
  Estadístico = c("Mínimo", "Mediana", "Media", "Percentil 90", "Percentil 99", "Máximo"),
  valor = c(min(df_muestra$S_ij), med_S, mean(df_muestra$S_ij),
            quantile(df_muestra$S_ij, c(0.90, 0.99)), max(df_muestra$S_ij))
) %>%
  transmute(Estadístico,
            `S_ij` = fmt_num(valor, 3),
            `Sobrecosto (%)` = paste0(fmt_num((valor - 1) * 100, 1), " %"))

kable(resumen_S, caption = "Distribución del índice de rodeo y sobrecosto espacial")
Distribución del índice de rodeo y sobrecosto espacial
Estadístico S_ij Sobrecosto (%)
Mínimo 1,000 0,0 %
Mediana 1,105 10,5 %
Media 1,216 21,6 %
Percentil 90 1,486 48,6 %
Percentil 99 2,762 176,2 %
Máximo 9,392 839,2 %
ggplot(df_muestra, aes(x = d_E, y = S_ij)) +
  geom_point(alpha = 0.3, color = "#e41a1c", size = 1) +
  geom_hline(yintercept = 1, linetype = "dashed", linewidth = 0.8) +
  geom_smooth(method = "gam", color = "#377eb8", se = TRUE) +
  scale_y_log10() +
  theme_minimal() +
  theme(panel.grid.minor = element_blank()) +
  labs(title = paste0("Índice de rodeo vs. distancia de gran círculo [n = ",
                      fmt_num(nrow(df_muestra), 0), "]"),
       x = "Distancia directa d_E (km)", y = "S_ij = d_red / d_E (escala log)")

En la mediana, volar a través de la red implica recorrer apenas un 10,5 % de kilómetros adicionales frente a la trayectoria en línea recta, lo que indica una elevada eficiencia espacial global. No obstante, el sesgo a la derecha revela pares periféricos con sobrecostos severos: el 1 % de los pares con mayor rodeo supera un sobrecosto del 176 %.


4.5 Conflicto de optimización: escalas vs. distancia

Para evaluar el costo de oportunidad entre reducir el número de conexiones o la distancia recorrida, se comparan dos políticas de enrutamiento:

  1. Ruta óptima en escalas: minimiza \(d_G\); se evalúa la distancia física total recorrida (km_minsaltos).
  2. Ruta óptima en kilómetros: minimiza \(d_{\text{red}}\); se evalúa la cantidad de escalas requeridas (saltos_minkm).

Ambas se obtienen con pesos lexicográficos, sin recorrer los caminos uno por uno:

  • Con pesos \(M + \ell_{uv}\) (\(M\) muy grande), el camino de menor peso es el de menos escalas y, entre los empatados, el de menos kilómetros. Los kilómetros se recuperan como \(\text{peso} - M\,d_G\).
  • Con pesos \(\ell_{uv} + \varepsilon\) (\(\varepsilon\) muy pequeño), el camino de menor peso es el de menos kilómetros, y su número de escalas se recupera como \((\text{peso} - d_{\text{red}})/\varepsilon\).
M   <- 1e6     # mayor que cualquier diferencia posible de km entre caminos
eps <- 1e-4    # km por escala: no altera en la práctica la elección del camino

dist_lex_s <- distances(giant_component, v = src, mode = "all", weights = km + M)
dist_eps_s <- distances(giant_component, v = src, mode = "all", weights = km + eps)

idx_m <- cbind(df_muestra$fila, df_muestra$j)

# Verificación: el número de saltos del camino lexicográfico coincide con d_G
stopifnot(all(floor(dist_lex_s[idx_m] / M) == df_muestra$d_G))

df_muestra <- df_muestra %>%
  mutate(
    km_minsaltos = dist_lex_s[idx_m] - M * d_G,
    saltos_minkm = round((dist_eps_s[idx_m] - d_red) / eps),
    penal_km     = km_minsaltos / d_red,
    saltos_extra = saltos_minkm - d_G
  )

stopifnot(all(df_muestra$km_minsaltos >= df_muestra$d_red - 1e-6),
          all(df_muestra$saltos_extra >= 0))

# Resúmenes que se usan en el texto
pct_coinc  <- mean(abs(df_muestra$penal_km - 1) < 1e-6)
med_penal  <- median(df_muestra$penal_km)
med_extra  <- median(df_muestra$saltos_extra)
max_penal  <- max(df_muestra$penal_km)

kable(
  data.frame(
    Métrica = c("Pares donde ambos caminos coinciden en kilómetros",
                "Pares con penalización en distancia > 50 % (penal_km > 1,5)",
                "Pares con al menos 2 saltos adicionales (saltos_extra ≥ 2)",
                "Mediana de penalización por minimizar saltos (penal_km)",
                "Mediana de saltos extra por minimizar kilómetros"),
    Valor = c(fmt_por(pct_coinc, 1),
              fmt_por(mean(df_muestra$penal_km > 1.5), 1),
              fmt_por(mean(df_muestra$saltos_extra >= 2), 1),
              fmt_num(med_penal, 3),
              med_extra)
  ),
  align = c("l", "r"),
  caption = "Trade-off entre minimización de saltos y minimización de distancia"
)
Trade-off entre minimización de saltos y minimización de distancia
Métrica Valor
Pares donde ambos caminos coinciden en kilómetros 43,3 %
Pares con penalización en distancia > 50 % (penal_km > 1,5) 1,8 %
Pares con al menos 2 saltos adicionales (saltos_extra ≥ 2) 28,1 %
Mediana de penalización por minimizar saltos (penal_km) 1,004
Mediana de saltos extra por minimizar kilómetros 1
p_a <- ggplot(df_muestra, aes(x = d_red, y = km_minsaltos)) +
  geom_point(alpha = 0.3, color = "#2b5c8f", size = 1) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "red") +
  theme_minimal() +
  labs(title = "Costo en distancia por minimizar escalas",
       x = "Distancia mínima sobre la red d_red (km)",
       y = "Distancia en ruta de mínimas escalas (km)")

p_b <- ggplot(df_muestra, aes(x = factor(saltos_extra))) +
  geom_bar(fill = "#e41a1c", alpha = 0.8) +
  theme_minimal() +
  labs(title = "Escalas adicionales por minimizar distancia",
       x = "Escalas extra (camino de km mínimos - d_G)", y = "Frecuencia de pares")

p_a | p_b

En el 43 % de los pares, el camino de menos escalas coincide exactamente con el de menos kilómetros, y cuando difieren la diferencia suele ser pequeña: la penalización mediana por minimizar escalas es de apenas 0,4 % en distancia, y solo el 1,8 % de los pares recorre más de un 50 % adicional. En sentido inverso, seguir la ruta más corta en kilómetros exige típicamente una escala adicional. El conflicto entre ambos criterios es, por tanto, poco frecuente, pero puede ser severo en pares periféricos.

4.5.1 Pares con altas discrepancias

limpiar_nombre <- function(x) {
  x %>%
    gsub("(?i)\\s+(International|Domestic|Regional|Municipal|County|National|Executive|Memorial|Intercontinental)?\\s*(Airport|Airfield|Airstrip|Air Base|Aerodrome|Heliport)$", "", ., perl = TRUE) %>%
    gsub("(?i)\\s+Airport", "", ., perl = TRUE) %>%
    gsub("\\s*\\([^)]+\\)", "", .) %>%
    trimws()
}

casos_dif_5 <- df_muestra %>%
  filter(penal_km > 1.5 | saltos_extra >= 2) %>%
  arrange(desc(penal_km)) %>%
  slice_head(n = 5) %>%
  mutate(
    Origen  = limpiar_nombre(V(giant_component)$name[i]),
    Destino = limpiar_nombre(V(giant_component)$name[j]),
    Par     = paste(Origen, "→", Destino)
  ) %>%
  select(
    `Origen → Destino`       = Par,
    `Escalas mín.`           = d_G,
    `Escalas en km mín.`     = saltos_minkm,
    `Dist. óptima red (km)`  = d_red,
    `Dist. mín. escalas (km)` = km_minsaltos,
    `Sobrecosto (x)`         = penal_km,
    `Escalas extra`          = saltos_extra
  )

kable(casos_dif_5, digits = 2,
      align = c("l", "c", "c", "r", "r", "c", "c"),
      format.args = list(decimal.mark = ",", big.mark = "."),
      caption = "Cinco pares con mayor conflicto entre minimizar escalas y distancia recorrida")
Cinco pares con mayor conflicto entre minimizar escalas y distancia recorrida
Origen → Destino Escalas mín. Escalas en km mín. Dist. óptima red (km) Dist. mín. escalas (km) Sobrecosto (x) Escalas extra
Cayenne-Rochambeau → General José Antonio Anzoategui 3 4 1.509,69 6.205,52 4,11 1
Sangster → Ignacio Agramonte 2 3 1.300,42 5.290,03 4,07 1
Ottawa Macdonald-Cartier → San Tomé 3 4 4.315,33 14.439,25 3,35 1
Stewart → Bishop 2 3 1.022,27 3.228,29 3,16 1
Juan Gualberto Gomez → Arenal 3 5 2.057,64 6.150,85 2,99 2

La inspección de los pares con mayor conflicto revela que priorizar la ruta con menor número de saltos genera desviaciones severas respecto al óptimo en kilómetros. En los cinco casos, el camino más corto en distancia requiere más tramos que el de mínimas escalas, por lo que el conflicto entre ambos criterios es directo:

  • El caso más pronunciado corresponde a Cayenne-Rochambeau \(\to\) General José Antonio Anzoátegui. El camino óptimo en red suma \(1.509{,}69\text{ km}\) con \(4\) tramos, mientras que la ruta de mínimas escalas (\(3\) tramos) alcanza los \(6.205{,}52\text{ km}\), un sobrecosto de \(4{,}11\times\) (\(311\ \%\) de desviación).
  • Un patrón similar aparece en Sangster \(\to\) Ignacio Agramonte y Stewart \(\to\) Bishop, pares relativamente cercanos en los que ahorrar un tramo (\(2\) en lugar de \(3\)) multiplica la distancia recorrida por \(4{,}07\times\) (\(5.290{,}03\text{ km}\) frente a \(1.300{,}42\text{ km}\)) y \(3{,}16\times\) (\(3.228{,}29\text{ km}\) frente a \(1.022{,}27\text{ km}\)), respectivamente.
  • En Ottawa Macdonald-Cartier \(\to\) San Tomé, aceptar \(1\) escala extra (\(4\) tramos en lugar de \(3\)) reduce el trayecto de \(14.439{,}25\text{ km}\) a \(4.315{,}33\text{ km}\) (\(3{,}35\times\)).
  • Finalmente, Juan Gualberto Gómez \(\to\) Arenal es el único caso que requiere \(2\) escalas extra (\(5\) tramos en lugar de \(3\)) para recorrer la ruta más corta. A cambio, la distancia baja de \(6.150{,}85\text{ km}\) a \(2.057{,}64\text{ km}\) (\(2{,}99\times\)).

Todos los pares extremos involucran aeropuertos regionales de América, en particular del Caribe, Centroamérica y el norte de Sudamérica. Allí, la ruta de menos escalas se desvía a través de hubs alejados del corredor geográfico directo, mientras que la ruta más corta encadena varios vuelos regionales.

5 Restricción espacial

Una red de transporte aéreo está sujeta a un costo físico: conectar aeropuertos lejanos es más caro que conectar aeropuertos cercanos. En esta sección se estudia (i) cómo decae la probabilidad de conexión con la distancia geográfica, (ii) qué tan distinta es la red observada de redes nulas que conservan los grados pero ignoran la geografía, y (iii) si los hubs mantienen rutas sistemáticamente más largas una vez se descuenta el efecto trivial de tener más aristas. Todo el análisis se realiza sobre la componente gigante.

5.1 Probabilidad de conexión según la distancia

Las distancias de gran círculo \(d_E\) entre todos los pares de aeropuertos de la componente gigante se agrupan en \(K = 25\) intervalos de ancho logarítmico (de 0 a unos 20.137 km, ligeramente por encima de la semicircunferencia terrestre). En cada intervalo \(b\) se estima

\[ \hat P(b) = \frac{E_b}{N_b}, \]

donde \(E_b\) es el número de aristas cuya longitud cae en \(b\) y \(N_b\) es el número de pares de aeropuertos (conectados o no) con \(d_E\) en \(b\). Los intervalos de confianza son de Clopper-Pearson. Para resumir la tasa de decaimiento se ajusta una regresión logística cuasi-binomial sobre los conteos,

\[ \operatorname{logit} P(d) = \alpha + \beta \log_{10} d, \]

de modo que \(\beta<0\) indica que la probabilidad de conexión cae, aproximadamente, como una potencia de la distancia. Se usa la familia cuasi-binomial porque las aristas que comparten un hub no son independientes, y la binomial subestimaría los errores estándar.

g_obs <- giant_component
stopifnot(is_simple(g_obs))

n_v    <- vcount(g_obs)
m_e    <- ecount(g_obs)
k_obs  <- degree(g_obs)
coords <- cbind(V(g_obs)$longitude, V(g_obs)$latitude)
cont   <- V(g_obs)$continent
cont[is.na(cont)] <- "Otros"

el_obs <- as_edgelist(g_obs, names = FALSE)

largo_km <- function(el) {
  distHaversine(coords[el[, 1], ], coords[el[, 2], ]) / 1000
}

prop_inter <- function(el) {
  ci <- cont[el[, 1]]; cj <- cont[el[, 2]]
  ok <- ci != "Otros" & cj != "Otros"
  mean(ci[ok] != cj[ok])
}

d_obs <- largo_km(el_obs)

# Intervalos logarítmicos de distancia
R_MAX  <- pi * 6378137 / 1000
brk    <- c(0, 10^seq(1, log10(R_MAX * 1.005), length.out = 25))
n_bins <- length(brk) - 1
ctr    <- sqrt(pmax(brk[-length(brk)], 1) * brk[-1])

bin_count <- function(d) tabulate(findInterval(d, brk), nbins = n_bins)

# Pares posibles (N_b) y aristas observadas (E_b) por intervalo
d_mat <- distm(coords, fun = distHaversine) / 1000
d_all <- d_mat[upper.tri(d_mat)]
rm(d_mat); invisible(gc())

pares_bin   <- bin_count(d_all)
aristas_bin <- bin_count(d_obs)
rm(d_all); invisible(gc())

stopifnot(sum(pares_bin) == choose(n_v, 2), sum(aristas_bin) == m_e)

5.2 Redes nulas por rewiring

Para decidir si el efecto de la distancia es notable, se compara la red observada con redes nulas: redes en las que cada aeropuerto conserva exactamente su número de conexiones, pero elige sus destinos sin considerar la geografía. Se generan mediante rewiring con preservación de grados:

  1. Se eligen al azar dos aristas, \((a,b)\) y \((c,d)\).
  2. Se reemplazan por \((a,d)\) y \((c,b)\), es decir, los extremos se “cruzan”.
  3. Si el intercambio crea un lazo (\(a=d\) o \(c=b\)) o una arista que ya existe, se descarta y se prueba con otro par.

Como cada aeropuerto pierde una conexión y gana otra, su grado no cambia, de modo que la secuencia de grados se preserva de forma exacta (y con ella el número de aristas y la densidad). Las coordenadas de los aeropuertos tampoco cambian, pero las aristas sí; por eso las longitudes se recalculan con Haversine en cada réplica. Dado que la mayoría de los pares de aeropuertos del mundo están muy lejos entre sí, los intercambios aleatorios tienden a alargar las rutas, y ese es el efecto que se mide.

Tras \(10\,m\) intentos de intercambio (\(m\) = número de aristas), la red es aproximadamente una muestra aleatoria de las redes simples con esa secuencia de grados. Se repite el proceso \(B = 500\) veces, partiendo siempre de la red observada. Si la red observada se parece a las nulas en una medida, esa medida puede explicarse solo por los grados; si difiere de forma marcada, la diferencia no se debe al grado y es coherente con una restricción espacial.

Para cada réplica \(r\) se calculan la longitud total \(L^{(r)}=\sum_{(i,j)\in E^{(r)}} d_E(i,j)\), la longitud promedio \(\bar L^{(r)}=L^{(r)}/m\) y la proporción de aristas intercontinentales. También se guardan el histograma de longitudes y la longitud media incidente de cada nodo, que se usan más adelante. Para un estadístico \(T\) con valor observado \(T_{obs}\) y réplicas \(T^{(1)},\dots,T^{(B)}\), se reporta \(z=(T_{obs}-\bar T)/s_T\) y el valor \(p\) empírico bilateral

\[ p = \frac{1 + \#\{r : |T^{(r)}-\bar T| \ge |T_{obs}-\bar T|\}}{B+1}, \]

cuyo mínimo posible es \(1/(B+1)\).

B <- 500

# Las réplicas se guardan en disco: solo se calculan la primera vez.
# Si cambia la red o el diseño, borre el archivo para recalcular.
archivo_nulos <- file.path("cache_resultados", paste0("nulos_rewiring_B", B, ".rds"))

# TRUE: usa varios núcleos (requiere 'future.apply'). Los resultados son
# reproducibles pero no idénticos a los de la versión secuencial.
USAR_PARALELO <- FALSE

g_lean <- make_graph(as.vector(t(el_obs)), n = n_v, directed = FALSE)

una_replica <- function(b) {
  g_rw <- rewire(g_lean, keeping_degseq(loops = FALSE, niter = 10 * m_e))
  el_b <- as_edgelist(g_rw, names = FALSE)
  d_b  <- largo_km(el_b)
  list(
    grados_ok = all(degree(g_rw) == k_obs),
    total     = sum(d_b),
    media     = mean(d_b),
    inter     = prop_inter(el_b),
    hist      = bin_count(d_b),
    mean_nodo = strength(g_rw, weights = d_b) / k_obs
  )
}

if (file.exists(archivo_nulos)) {
  nulos <- readRDS(archivo_nulos)
} else {
  if (USAR_PARALELO && requireNamespace("future.apply", quietly = TRUE)) {
    future::plan(future::multisession, workers = max(1, parallel::detectCores() - 1))
    nulos <- future.apply::future_lapply(seq_len(B), una_replica,
                                         future.seed = 2026,
                                         future.packages = c("igraph", "geosphere"))
    future::plan(future::sequential)
  } else {
    set.seed(2026)
    nulos <- lapply(seq_len(B), una_replica)
  }
  saveRDS(nulos, archivo_nulos)
}

stopifnot(all(vapply(nulos, `[[`, logical(1), "grados_ok")))   # grados intactos

res <- data.frame(
  total = vapply(nulos, `[[`, numeric(1), "total"),
  media = vapply(nulos, `[[`, numeric(1), "media"),
  inter = vapply(nulos, `[[`, numeric(1), "inter")
)
hist_null      <- do.call(rbind, lapply(nulos, `[[`, "hist"))
mean_nodo_null <- do.call(rbind, lapply(nulos, `[[`, "mean_nodo"))

5.2.1 Probabilidad de conexión y distribución de longitudes: observada vs. nula

df_prob <- data.frame(ctr, pares = pares_bin, aristas = aristas_bin,
                      nulo = colMeans(hist_null)) %>%
  filter(pares > 0) %>%
  mutate(
    p_obs  = aristas / pares,
    p_nulo = nulo / pares,
    ic_lo  = ifelse(aristas == 0, 0, qbeta(0.025, aristas, pares - aristas + 1)),
    ic_hi  = qbeta(0.975, aristas + 1, pares - aristas)
  )

p_prob <- ggplot(filter(df_prob, aristas > 0), aes(x = ctr)) +
  geom_errorbar(aes(ymin = ic_lo, ymax = ic_hi), width = 0.04, color = "grey50") +
  geom_point(aes(y = p_obs, color = "Observada"), size = 2.5) +
  geom_line(data = filter(df_prob, nulo > 0),
            aes(y = p_nulo, color = "Nula (rewiring)"), linewidth = 0.9) +
  scale_x_log10(labels = comma) +
  scale_y_log10() +
  scale_color_manual(values = c("Observada" = "#e41a1c", "Nula (rewiring)" = "grey30"), name = NULL) +
  labs(title = "A. Probabilidad de conexión vs. distancia",
       x = "Distancia de gran círculo d_E (km, escala log)",
       y = "P(conexión) = aristas / pares posibles (log)") +
  theme_minimal() +
  theme(panel.grid.minor = element_blank(), legend.position = "bottom",
        plot.title = element_text(face = "bold"))

prop_null <- hist_null / m_e
q_null    <- apply(prop_null, 2, quantile, probs = c(0.025, 0.975))

df_hist <- data.frame(ctr, obs = aristas_bin / m_e, nulo = colMeans(prop_null),
                      lo = q_null[1, ], hi = q_null[2, ])

p_hist <- ggplot(df_hist, aes(x = ctr)) +
  geom_ribbon(aes(ymin = lo, ymax = hi, fill = "Nula: IC 95 %"), alpha = 0.6) +
  geom_line(aes(y = nulo, color = "Nula: media"), linewidth = 0.8) +
  geom_line(aes(y = obs, color = "Observada"), linewidth = 1) +
  geom_point(aes(y = obs, color = "Observada"), size = 2) +
  scale_x_log10(labels = comma) +
  scale_fill_manual(values = c("Nula: IC 95 %" = "grey75"), name = NULL) +
  scale_color_manual(values = c("Nula: media" = "grey30", "Observada" = "#e41a1c"), name = NULL) +
  labs(title = "B. Distribución de longitudes de las rutas",
       x = "Longitud de la ruta (km, escala log)", y = "Proporción de aristas en el intervalo") +
  theme_minimal() +
  theme(panel.grid.minor = element_blank(), legend.position = "bottom",
        plot.title = element_text(face = "bold"))

p_prob | p_hist

La figura A muestra la probabilidad de conexión en función de la distancia de gran círculo. En la red observada, esta probabilidad es alta a distancias cortas (del orden de 0,5 entre aeropuertos separados por unas decenas de kilómetros) y decrece de forma sostenida hasta valores cercanos a \(10^{-5}\) en las distancias más largas. En las redes nulas, en cambio, la probabilidad es mucho menor a distancias cortas y permanece aproximadamente constante a partir de unos 100 km. Las dos curvas se cruzan en torno a los 3.500 km: por debajo de esa distancia la red observada presenta más conexiones que la referencia nula, y por encima presenta menos. Los intervalos de confianza son amplios en los primeros intervalos de distancia, donde hay muy pocos pares de aeropuertos.

La figura B compara la distribución de longitudes de las rutas. La distribución observada es unimodal, con moda cercana a 1.400 km, y cae con rapidez a partir de unos 3.000 km, aunque conserva una pequeña masa de rutas entre 6.000 y 9.000 km. La distribución bajo rewiring está desplazada hacia distancias mucho mayores, con moda cercana a 9.000 km. La banda del 95 % de las réplicas nulas es muy estrecha, por lo que la diferencia entre ambas distribuciones es mucho mayor que la variabilidad entre réplicas. Dado que los intervalos son logarítmicos, la comparación debe hacerse por la forma y la posición de las curvas, no por la altura absoluta de cada punto.

ajusta_beta <- function(e) {
  d <- data.frame(e = e, ne = pares_bin - e, x = log10(ctr))[pares_bin > 0, ]
  coef(suppressWarnings(glm(cbind(e, ne) ~ x, family = quasibinomial, data = d)))[["x"]]
}

d_fit    <- data.frame(e = aristas_bin, ne = pares_bin - aristas_bin, x = log10(ctr))[pares_bin > 0, ]
fit_obs  <- glm(cbind(e, ne) ~ x, family = quasibinomial, data = d_fit)
beta_obs <- coef(fit_obs)[["x"]]
ci_beta  <- confint.default(fit_obs)["x", ]
beta_nul <- apply(hist_null, 1, ajusta_beta)

kable(
  data.frame(
    Red = c("Observada", "Nula (rewiring, media de las réplicas)"),
    `Pendiente β` = c(fmt_num(beta_obs, 3), fmt_num(mean(beta_nul), 3)),
    `IC 95 %` = c(paste0("[", fmt_num(ci_beta[1], 3), "; ", fmt_num(ci_beta[2], 3), "]"),
                  paste0("[", fmt_num(quantile(beta_nul, .025), 3), "; ",
                         fmt_num(quantile(beta_nul, .975), 3), "]")),
    check.names = FALSE
  ),
  align = c("l", "r", "r"),
  caption = "Decaimiento de la probabilidad de conexión con log10(d): regresión logística cuasi-binomial."
)
Decaimiento de la probabilidad de conexión con log10(d): regresión logística cuasi-binomial.
Red Pendiente β IC 95 %
Observada -3,478 [-3,946; -3,010]
Nula (rewiring, media de las réplicas) -0,771 [-0,802; -0,739]

5.2.2 Longitud total, longitud promedio y proporción intercontinental

Con la clasificación de countrycode, “Americas” es un único continente, por lo que las rutas entre Norteamérica y Sudamérica no cuentan como intercontinentales.

obs_stats <- c(total = sum(d_obs), media = mean(d_obs), inter = prop_inter(el_obs))

comp_tab <- do.call(rbind, lapply(names(obs_stats), function(s) {
  x0 <- res[[s]]; xo <- obs_stats[[s]]
  data.frame(
    estadistico = s, obs = xo,
    media_nula  = mean(x0), sd_nula = sd(x0),
    q025 = unname(quantile(x0, 0.025)), q975 = unname(quantile(x0, 0.975)),
    z = (xo - mean(x0)) / sd(x0),
    p = (1 + sum(abs(x0 - mean(x0)) >= abs(xo - mean(x0)))) / (B + 1)
  )
}))
rownames(comp_tab) <- comp_tab$estadistico

fmt_stat <- list(
  total = function(x) fmt_num(x, 0),
  media = function(x) fmt_num(x, 1),
  inter = function(x) fmt_por(x, 1)
)
f <- function(col) {
  vapply(seq_len(nrow(comp_tab)),
         function(i) fmt_stat[[comp_tab$estadistico[i]]](comp_tab[[col]][i]),
         character(1))
}
p_min <- 1 / (B + 1)
fmt_p <- function(p) ifelse(p <= p_min + 1e-12, paste0("≤ ", fmt_num(p_min, 3)), fmt_num(p, 3))

kable(
  data.frame(
    Estadístico = c("Longitud total (km)", "Longitud promedio (km)", "Proporción intercontinental"),
    Observado   = f("obs"),
    `Media nula` = f("media_nula"),
    `IC 95 % nulo` = paste0("[", f("q025"), "; ", f("q975"), "]"),
    z = fmt_num(comp_tab$z, 1),
    `Valor p` = fmt_p(comp_tab$p),
    check.names = FALSE
  ),
  align = c("l", "r", "r", "r", "r", "r"),
  caption = paste0("Red observada vs. ", B, " réplicas con secuencia de grados preservada.")
)
Red observada vs. 500 réplicas con secuencia de grados preservada.
Estadístico Observado Media nula IC 95 % nulo z Valor p
Longitud total (km) 33.367.658 141.321.478 [140.318.584; 142.196.241] -225,5 ≤ 0,002
Longitud promedio (km) 1.771,8 7.503,9 [7.450,7; 7.550,4] -225,5 ≤ 0,002
Proporción intercontinental 14,3 % 72,0 % [71,4 %; 72,6 %] -190,5 ≤ 0,002

5.3 Grado y longitud de las rutas incidentes

Para el nodo \(i\) de grado \(k_i\) se define la longitud media incidente \(\bar\ell_i = \frac{1}{k_i}\sum_{j\sim i} d_E(i,j)\). Al relacionar el grado con la longitud de las rutas hay dos efectos mecánicos que no dicen nada sobre los hubs:

  • La longitud total incidente \(\sum_j d_E(i,j)\) crece linealmente con \(k_i\) por simple suma.
  • La longitud máxima incidente crece con \(k_i\) porque el máximo de más observaciones es mayor.

La media \(\bar\ell_i\) no depende mecánicamente de \(k_i\) en esperanza, así que es el estadístico adecuado. Aun así, un aeropuerto grande puede tener rutas largas simplemente por estar lejos de otros (por ejemplo, Honolulu), no por ser hub. Por eso se controla con la red nula: como el rewiring conserva el grado y la ubicación de cada nodo, la longitud media esperada \(E_0[\bar\ell_i]\) bajo conexión sin preferencia geográfica ya incorpora tanto \(k_i\) como la posición del aeropuerto. El exceso de longitud se mide con el cociente

\[ \rho_i = \frac{\bar\ell_i^{\,obs}}{E_0[\bar\ell_i]}, \]

donde \(\rho_i<1\) indica rutas más cortas de lo que daría el azar. Se evalúa si \(\rho_i\) depende del grado mediante la correlación de Spearman y una regresión \(\log\rho_i = a + c\,\log k_i\). Además, para cada clase de grado se compara el promedio observado de \(\bar\ell_i\) con su distribución bajo las \(B\) réplicas, y se agrega una clase de hubs definida con el criterio de Tukey (\(k \ge Q_3 + 1{,}5\,\mathrm{IQR}\)) usado en las visualizaciones.

inc <- data.frame(nodo = c(el_obs[, 1], el_obs[, 2]), d = rep(d_obs, 2)) %>%
  group_by(nodo) %>%
  summarise(k = n(), total = sum(d), maxi = max(d), media = mean(d), .groups = "drop") %>%
  arrange(nodo)
stopifnot(nrow(inc) == n_v, all(inc$k == k_obs))

mean_obs <- inc$media
mean_exp <- colMeans(mean_nodo_null)
razon    <- mean_obs / mean_exp

# 1. Efectos mecánicos vs. media
rho_s <- function(y) suppressWarnings(cor(k_obs, y, method = "spearman"))
kable(
  data.frame(
    Medida = c("Longitud total incidente (suma)", "Longitud máxima incidente",
               "Longitud media incidente", "Cociente ρ = media observada / esperada (nula)"),
    `Spearman con el grado` = fmt_num(c(rho_s(inc$total), rho_s(inc$maxi),
                                        rho_s(mean_obs), rho_s(razon)), 3),
    check.names = FALSE
  ),
  align = c("l", "r"),
  caption = "Relación con el grado: efectos mecánicos frente a medidas controladas."
)
Relación con el grado: efectos mecánicos frente a medidas controladas.
Medida Spearman con el grado
Longitud total incidente (suma) 0,892
Longitud máxima incidente 0,770
Longitud media incidente 0,635
Cociente ρ = media observada / esperada (nula) 0,632
# 2. Regresión log-log del cociente
fit_rho <- lm(log(razon) ~ log(k_obs))
ci_rho  <- confint(fit_rho)["log(k_obs)", ]


# 3. Contraste por clase de grado
clase <- cut(k_obs, breaks = c(0, 1, 2, 5, 10, 20, 50, Inf),
             labels = c("1", "2", "3-5", "6-10", "11-20", "21-50", ">50"))
umbral_hub <- quantile(k_obs, 0.75) + 1.5 * IQR(k_obs)

test_clase <- function(idx, etiqueta) {
  T_obs <- mean(mean_obs[idx])
  T_nul <- rowMeans(mean_nodo_null[, idx, drop = FALSE])
  data.frame(
    Clase = etiqueta, n = length(idx), grado_med = mean(k_obs[idx]),
    obs = T_obs, nulo = mean(T_nul),
    lo = unname(quantile(T_nul, 0.025)), hi = unname(quantile(T_nul, 0.975)),
    razon = T_obs / mean(T_nul),
    p = (1 + sum(abs(T_nul - mean(T_nul)) >= abs(T_obs - mean(T_nul)))) / (B + 1)
  )
}

idx_clase <- Filter(length, split(seq_len(n_v), clase))
tab_clase <- do.call(rbind, Map(test_clase, idx_clase, names(idx_clase)))
tab_hub   <- test_clase(which(k_obs >= umbral_hub),
                        paste0("Hubs (k ≥ ", fmt_num(umbral_hub, 1), ")"))

tab_all <- rbind(tab_clase, tab_hub)
kable(
  data.frame(
    `Clase de grado` = tab_all$Clase,
    Nodos = fmt_num(tab_all$n, 0),
    `Grado medio` = fmt_num(tab_all$grado_med, 1),
    `Longitud media observada (km)` = fmt_num(tab_all$obs, 0),
    `Esperada bajo nula (km)` = fmt_num(tab_all$nulo, 0),
    Cociente = fmt_num(tab_all$razon, 3),
    `Valor p` = fmt_p(tab_all$p),
    check.names = FALSE
  ),
  align = c("l", rep("r", 6)),
  caption = "Longitud media de las rutas incidentes por clase de grado: observada vs. rewiring."
)
Longitud media de las rutas incidentes por clase de grado: observada vs. rewiring.
Clase de grado Nodos Grado medio Longitud media observada (km) Esperada bajo nula (km) Cociente Valor p
1 705 1,0 485 8.474 0,057 ≤ 0,002
2 638 2,0 521 8.483 0,061 ≤ 0,002
3-5 760 3,7 669 8.369 0,080 ≤ 0,002
6-10 389 7,7 1.069 8.005 0,134 ≤ 0,002
11-20 263 15,0 1.343 7.707 0,174 ≤ 0,002
21-50 249 32,8 1.632 7.642 0,214 ≤ 0,002
>50 184 96,7 2.192 7.161 0,306 ≤ 0,002
Hubs (k ≥ 19,5) 455 58,0 1.845 7.450 0,248 ≤ 0,002
# Gráficos
df_nodo <- data.frame(k = k_obs, media = mean_obs, razon = razon)

p_a <- ggplot(df_nodo, aes(x = k, y = media)) +
  geom_point(alpha = 0.25, color = "#2b5c8f", size = 1) +
  geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"), color = "#e41a1c", se = TRUE) +
  scale_x_log10() + scale_y_log10(labels = comma) +
  labs(title = "A. Longitud media incidente vs. grado (observada)",
       x = "Grado k (escala log)", y = "Longitud media (km, escala log)") +
  theme_minimal() + theme(panel.grid.minor = element_blank())

p_b <- ggplot(df_nodo, aes(x = k, y = razon)) +
  geom_point(alpha = 0.25, color = "#2b5c8f", size = 1) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"), color = "#e41a1c", se = TRUE) +
  scale_x_log10() + scale_y_log10() +
  labs(title = "B. Cociente observada / esperada bajo rewiring",
       x = "Grado k (escala log)", y = "ρ (escala log); 1 = igual que el azar") +
  theme_minimal() + theme(panel.grid.minor = element_blank())

tab_clase_plot <- tab_clase
tab_clase_plot$Clase <- factor(tab_clase_plot$Clase, levels = levels(clase))
p_c <- ggplot(tab_clase_plot, aes(x = Clase)) +
  geom_errorbar(aes(ymin = lo, ymax = hi), width = 0.25, color = "grey40") +
  geom_point(aes(y = nulo, shape = "Nula: media e IC 95 %"), size = 3, color = "grey30") +
  geom_point(aes(y = obs, shape = "Observada"), size = 3, color = "#e41a1c") +
  scale_shape_manual(values = c("Nula: media e IC 95 %" = 1, "Observada" = 16), name = NULL) +
  scale_y_continuous(labels = comma) +
  labs(title = "C. Longitud media incidente por clase de grado",
       x = "Clase de grado", y = "Longitud media (km)") +
  theme_minimal() + theme(panel.grid.minor = element_blank(), legend.position = "bottom")

(p_a | p_b) / p_c

La regresión de \(\log\rho_i\) sobre \(\log k_i\) da una pendiente de 0,50 (IC 95 %: [0,48; 0,52]): al duplicar el grado, el cociente entre la longitud observada y la esperada aumenta en promedio un 41 %.

La figura A muestra que la longitud media de las rutas incidentes crece con el grado: la curva suavizada pasa de unos 300 km en los aeropuertos de grado 1 a cerca de 3.000 km en los de mayor grado. La dispersión es muy amplia entre los aeropuertos pequeños (de menos de 100 km a más de 10.000 km) y se reduce con el grado, algo esperable porque un promedio sobre más rutas varía menos. La figura B muestra el cociente entre la longitud observada y la esperada bajo rewiring: casi todos los aeropuertos quedan por debajo de 1 y el cociente aumenta con el grado, de aproximadamente 0,04 en grado 1 a cerca de 0,45 en los hubs.

La figura C resume lo anterior por clase de grado. La longitud media esperada bajo rewiring disminuye con el grado (de unos 8.474 km en grado 1 a unos 7.161 km en la clase de más de 50 conexiones), porque los hubs están ubicados en regiones centrales, más cercanas al resto del mundo. La longitud observada, en cambio, aumenta (de unos 485 km a unos 2.192 km). Las redes nulas quedan muy por encima de lo observado en todas las clases, con intervalos del 95 % muy estrechos.

6 Caracterización de los aeropuertos

A continuación se mide la importancia de cada aeropuerto de la componente gigante con siete medidas de centralidad y se comparan los rankings que producen. Después se calcula el coeficiente de participación continental para distinguir los hubs que operan dentro de una región de los que conectan continentes.

Las medidas utilizadas son:

Medida Definición Qué captura
Grado \(d_i=\sum_{j\neq i}y_{ij}\) Conectividad directa (número de destinos)
Fuerza \(s_i=\sum_{j}w_{ij}\) Suma de los pesos de las rutas incidentes
Centralidad armónica \(C_{H}(i)=\frac{1}{n-1}\sum_{j\neq i}\frac{1}{d_G(i,j)}\) Cercanía (en saltos) al resto de la red
Intermediación \(C_{B}(i)=\frac{2}{(n-1)(n-2)}\sum_{s\neq i\neq t}\frac{\sigma_{st}(i)}{\sigma_{st}}\) Papel de puente en los caminos más cortos
Centralidad propia \(C_{E}(i)=\frac{1}{\lambda_1}\sum_{j}y_{ij}\,C_{E}(j)\) Conexión con aeropuertos que a su vez son centrales
PageRank \(CP(i)=\frac{1-\delta}{n}+\delta\sum_{j\sim i}\frac{CP(j)}{d_j},\ \delta=0{,}85\) Prestigio según un paseo aleatorio
Coreness mayor \(k\) tal que \(i\) pertenece al \(k\)-core Pertenencia al núcleo más denso

Donde \(\mathbf{Y}=[y_{ij}]\) es la matriz de adyacencia, \(\lambda_1\) su mayor valor propio, \(\sigma_{st}\) el número de caminos más cortos entre \(s\) y \(t\), y \(\sigma_{st}(i)\) el número de esos caminos que pasan por \(i\).

Todas las medidas se calculan sobre la red no dirigida, sin pesos, salvo la fuerza. El peso de las aristas (número de registros ruta–aerolínea) mide intensidad, no costo, así que usarlo como distancia en la intermediación o la armónica invertiría su significado: un peso alto haría “más lejos” a un par muy conectado.

g      <- giant_component
n_v    <- vcount(g)
el     <- as_edgelist(g, names = FALSE)
coords <- cbind(V(g)$longitude, V(g)$latitude)
cont   <- V(g)$continent
cont[is.na(cont)] <- "Otros"
etq    <- ifelse(is.na(V(g)$iata) | V(g)$iata == "", V(g)$icao, V(g)$iata)

meds_id  <- c("grado", "fuerza", "armonica", "intermediacion", "propia", "pagerank", "coreness")
meds_lab <- c("Grado", "Fuerza", "Armónica", "Intermediación", "Propia", "PageRank", "Coreness")

fmt_valor <- list(
  grado          = function(x) fmt_num(x, 0),
  fuerza         = function(x) fmt_num(x, 0),
  armonica       = function(x) fmt_num(x, 2),
  intermediacion = function(x) fmt_por(x, 1),
  propia         = function(x) fmt_num(x, 2),
  pagerank       = function(x) fmt_por(x, 2),
  coreness       = function(x) fmt_num(x, 0)
)

df_med <- data.frame(
  iata           = etq,
  ciudad         = V(g)$city,
  pais           = V(g)$country,
  continente     = cont,
  grado          = degree(g),
  fuerza         = strength(g, weights = E(g)$num_registros),
  armonica       = harmonic_centrality(g, mode = "all", weights = NA, normalized = TRUE),
  intermediacion = betweenness(g, directed = FALSE, weights = NA, normalized = TRUE),
  propia         = eigen_centrality(g, directed = FALSE, scale = TRUE, weights = NA)$vector,
  pagerank       = page_rank(g, directed = FALSE, damping = 0.85, weights = NA)$vector,
  coreness       = coreness(g, mode = "all"),
  stringsAsFactors = FALSE
)

top5 <- sapply(meds_id, function(m) {
  o <- order(-df_med[[m]], -df_med$grado)[1:5]
  paste0(df_med$iata[o], " (", fmt_valor[[m]](df_med[[m]][o]), ")")
})
colnames(top5) <- meds_lab
rownames(top5) <- 1:5

kable(top5, align = "l",
      caption = "Cinco primeros aeropuertos (IATA) y valor de cada medida de centralidad. Intermediación y PageRank en %.")
Cinco primeros aeropuertos (IATA) y valor de cada medida de centralidad. Intermediación y PageRank en %.
Grado Fuerza Armónica Intermediación Propia PageRank Coreness
AMS (248) ATL (1.826) FRA (0,47) CDG (6,5 %) AMS (1,00) ATL (0,49 %) AMS (31)
FRA (244) ORD (1.108) CDG (0,46) LAX (6,1 %) FRA (0,99) ORD (0,45 %) FRA (31)
CDG (240) PEK (1.061) AMS (0,46) ANC (5,9 %) CDG (0,95) ISL (0,45 %) CDG (31)
ISL (233) LHR (1.047) LHR (0,45) DXB (5,6 %) MUC (0,88) DFW (0,45 %) ISL (31)
ATL (217) CDG (1.041) DXB (0,45) FRA (5,3 %) FCO (0,81) DEN (0,44 %) MUC (31)

Las medidas no coinciden en quién encabeza la red, y cada una destaca un tipo de importancia distinto. Las que miden conectividad directa y pertenencia al núcleo denso (grado, centralidad propia y coreness), junto con la cercanía al resto de la red (armónica), están lideradas por hubs europeos: Ámsterdam (AMS), Fráncfort (FRA) y París (CDG) aparecen en las cuatro. La fuerza, que suma los registros de ruta, la encabezan Atlanta (ATL), Chicago (ORD) y Pekín (PEK), mientras que el PageRank, que premia estar conectado a aeropuertos importantes, lo encabezan hubs estadounidenses (ATL, ORD, DFW, DEN) y Estambul (ISL). La intermediación, que mide el papel de puente, muestra otro perfil: Los Ángeles (LAX) y Anchorage (ANC) aparecen solo en esta medida, lo que sugiere que su importancia está en enlazar regiones y no en el número de destinos.

6.1 Comparación de rankings

Los rankings se comparan de dos maneras complementarias. La correlación de Spearman compara el orden de todos los aeropuertos, pero está dominada por la gran masa de aeropuertos periféricos. La superposición de los primeros \(k\) mide cuántos aeropuertos aparecen en el top-\(k\) de dos medidas, que es lo que importa para identificar hubs: dos medidas pueden estar muy correlacionadas globalmente y discrepar en quiénes encabezan la lista.

Elección de \(k\). Inicialmente se consideró usar el número de aeropuertos que superan el umbral de Tukey sobre el grado (\(d_v\ge Q_3+1{,}5\,\mathrm{IQR}\)), como en las visualizaciones, pero ese número resulta demasiado grande para representar solo a los hubs. Por ello se toma \(k=\lceil 0{,}01\,n\rceil\) (el 1 % de los aeropuertos, 32 en este caso) como valor principal. Es lo bastante pequeño para que la lista corresponda a los hubs reales, pues con \(k\) grande la superposición crece solo porque los conjuntos cubren buena parte de la red (por azar se esperaría una superposición de \(k/n\)), y lo bastante grande para no depender de uno o dos intercambios de posición. Se evalúa además la sensibilidad con \(k\in\{10, 25, 32, 50, 100, 5\,\%\ \text{de}\ n\}\). Los empates se rompen por grado y fuerza; esto afecta sobre todo a la coreness, que es una medida gruesa con muchos aeropuertos en cada clase.

rho <- cor(df_med[, meds_id], method = "spearman")
dimnames(rho) <- list(meds_lab, meds_lab)

top_set <- function(x, k) order(-x, -df_med$grado, -df_med$fuerza)[seq_len(k)]

solap <- function(k) {
  S <- lapply(meds_id, function(m) top_set(df_med[[m]], k))
  M <- outer(seq_along(S), seq_along(S),
             Vectorize(function(a, b) length(intersect(S[[a]], S[[b]])) / k))
  dimnames(M) <- list(meds_lab, meds_lab)
  M
}

k_main <- ceiling(0.01 * n_v)
k_grid <- sort(unique(c(10, 25, k_main, 50, 100, ceiling(0.05 * n_v))))
M_main <- solap(k_main)


sens <- do.call(rbind, lapply(k_grid, function(k) {
  M <- solap(k); v <- M[upper.tri(M)]
  data.frame(k = k, media = mean(v), minimo = min(v), maximo = max(v), azar = k / n_v)
}))

# Resúmenes que se usan en el texto
rho_rango <- range(rho[upper.tri(rho)])
bloque1   <- rho[c("Grado", "Fuerza", "Coreness"), c("Grado", "Fuerza", "Coreness")]
rho_b1    <- range(bloque1[upper.tri(bloque1)])

heat <- function(M, titulo, lim) {
  d <- as.data.frame(as.table(M), stringsAsFactors = FALSE)
  names(d) <- c("x", "y", "v")
  d$x <- factor(d$x, levels = colnames(M))
  d$y <- factor(d$y, levels = rev(rownames(M)))
  ggplot(d, aes(x, y, fill = v)) +
    geom_tile(color = "white") +
    geom_text(aes(label = formatC(v, format = "f", digits = 2, decimal.mark = ",")), size = 3) +
    scale_fill_gradient(low = "white", high = "#2b5c8f", limits = lim, name = NULL) +
    labs(title = titulo, x = NULL, y = NULL) +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 40, hjust = 1), panel.grid = element_blank(),
          plot.title = element_text(face = "bold", size = 11))
}

pA <- heat(rho, "A. Correlación de Spearman (todos los aeropuertos)", c(min(0, min(rho)), 1))
pB <- heat(M_main, paste0("B. Superposición del top-", k_main), c(0, 1))

pC <- ggplot(sens, aes(x = k)) +
  geom_ribbon(aes(ymin = minimo, ymax = maximo), fill = "grey80", alpha = 0.7) +
  geom_line(aes(y = media), color = "#2b5c8f", linewidth = 1) +
  geom_point(aes(y = media), color = "#2b5c8f", size = 2) +
  geom_line(aes(y = azar), linetype = "dashed", color = "grey40") +
  geom_vline(xintercept = k_main, linetype = "dotted") +
  scale_x_log10() +
  labs(title = "C. Superposición según k (media y rango entre pares de medidas)",
       subtitle = "Línea punteada horizontal: superposición esperada por azar (k/n)",
       x = "k (escala log)", y = "Proporción de aeropuertos en común") +
  theme_minimal() + theme(panel.grid.minor = element_blank(),
                          plot.title = element_text(face = "bold", size = 11))

(pA | pB) / pC

Las medidas de centralidad están correlacionadas en todo el ranking (Spearman de 0,42 a 0,98) y forman dos bloques: grado, fuerza y coreness (0,94–0,98), y armónica y centralidad propia (0,97). La intermediación es la más distinta de las demás. La coincidencia en el extremo superior es mucho menor que la correlación global: grado y coreness (0,98) comparten solo el 62 % del top-32, y las menores coincidencias se dan entre la intermediación y la coreness (0,28) y entre la intermediación y la propia (0,31). El patrón se mantiene al variar \(k\): la superposición media se ubica entre 0,57 y 0,69 desde \(k = 10\) hasta el 5 % de los aeropuertos, siempre muy por encima de la esperada por azar.

Las comparaciones que involucran a la coreness deben leerse con cautela: 93 aeropuertos comparten su valor máximo (31), más que el propio \(k = 32\), por lo que su top-32 queda definido por el criterio de desempate (grado y fuerza). Precisamente por ese desempate, la baja coincidencia con el grado es informativa: varios de los aeropuertos con más destinos no pertenecen al núcleo de máxima coreness, lo que indica que parte de su grado proviene de conexiones con aeropuertos poco interconectados entre sí.

6.2 Aeropuertos con posiciones muy distintas entre medidas

Para cada aeropuerto se calcula su posición en cada ranking (1 = el más central; los empates reciben la posición promedio) y se mide la discrepancia como la diferencia entre su mejor y su peor posición. Para no destacar ruido de la periferia, solo se consideran los aeropuertos que están entre los 100 primeros en al menos una medida. La coreness, al ser una medida gruesa, tiende a ampliar la discrepancia.

rk <- function(x) rank(-x, ties.method = "average")
R  <- sapply(meds_id, function(m) rk(df_med[[m]]))
colnames(R) <- meds_lab

# Posición de un aeropuerto (por código IATA) en una medida; se usa en el texto
pos <- function(codigo, medida) {
  i <- match(codigo, df_med$iata)
  if (is.na(i)) NA else round(R[i, medida])
}

k_cand <- 100
cand   <- which(apply(R, 1, min) <= k_cand)
spread <- apply(R, 1, max) - apply(R, 1, min)
mejor  <- meds_lab[apply(R, 1, which.min)]
peor   <- meds_lab[apply(R, 1, which.max)]
n_tab <- 8
idx    <- cand[order(-spread[cand])][1:n_tab]

kable(
  data.frame(Aeropuerto = df_med$iata[idx], Ciudad = df_med$ciudad[idx], País = df_med$pais[idx],
             round(R[idx, ]), `Mejor en` = mejor[idx], `Peor en` = peor[idx],
             check.names = FALSE),
  row.names = FALSE,
  caption = paste0("Los ", n_tab, " aeropuertos con mayor discrepancia de posición entre medidas (posición 1 = más central).")
)
Los 8 aeropuertos con mayor discrepancia de posición entre medidas (posición 1 = más central).
Aeropuerto Ciudad País Grado Fuerza Armónica Intermediación Propia PageRank Coreness Mejor en Peor en
BET Bethel United States 408 572 2331 59 2526 78 1506 Intermediación Propia
HIR Honiara Solomon Islands 468 556 1739 85 2217 129 1195 Intermediación Propia
VLI Port-vila Vanuatu 424 492 1358 76 1771 118 1020 Intermediación Propia
FAI Fairbanks United States 496 572 1192 67 1682 125 1506 Intermediación Propia
YZF Yellowknife Canada 496 476 1680 68 1568 163 1195 Intermediación Armónica
YQT Thunder Bay Canada 600 770 917 63 1284 333 1195 Intermediación Propia
AEP Buenos Aires Argentina 218 256 1098 53 1223 46 723 PageRank Propia
PPT Papeete French Polynesia 328 396 524 47 1037 95 1020 Intermediación Propia
dif_gb  <- abs(log(R[, "Grado"] / R[, "Intermediación"]))
cand_gb <- which(pmin(R[, "Grado"], R[, "Intermediación"]) <= k_cand)
lab_gb  <- cand_gb[order(-dif_gb[cand_gb])][1:10]

df_rr <- data.frame(rk_grado = R[, "Grado"], rk_inter = R[, "Intermediación"], etq = "")
df_rr$etq[lab_gb] <- df_med$iata[lab_gb]

ggplot(df_rr, aes(rk_grado, rk_inter)) +
  geom_point(alpha = 0.25, color = "#2b5c8f", size = 1) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
  geom_point(data = df_rr[lab_gb, ], color = "#e41a1c", size = 2) +
  geom_text_repel(aes(label = etq), size = 3, max.overlaps = Inf, min.segment.length = 0) +
  scale_x_log10() + scale_y_log10() +
  labs(title = "Posición por grado vs. posición por intermediación",
       subtitle = "Los puntos por debajo de la diagonal son más centrales como puentes que por su número de destinos",
       x = "Posición por grado (1 = mayor; escala log)",
       y = "Posición por intermediación (1 = mayor; escala log)") +
  theme_minimal() + theme(panel.grid.minor = element_blank())

En la tabla, 7 de los 8 aeropuertos con mayor discrepancia alcanzan su mejor posición en intermediación (entre los puestos 47 y 85), mientras que en grado ocupan posiciones intermedias (entre 218 y 600) y en centralidad propia y armónica quedan muy atrás, en muchos casos por debajo del puesto 1.000. Todos son aeropuertos periféricos que funcionan como única puerta de entrada a una subred aislada: Bethel y Fairbanks en Alaska, Yellowknife y Thunder Bay en el norte de Canadá, Honiara y Port Vila en las islas del Pacífico y Papeete en la Polinesia Francesa. Tienen pocos destinos y no están conectados con el núcleo de la red, por lo que la centralidad propia y la armónica los ubican en la periferia, pero una gran cantidad de caminos más cortos pasa necesariamente por ellos, lo que eleva su intermediación. La única excepción es Buenos Aires-Aeroparque (AEP), cuya mejor posición es en PageRank (puesto 46): concentra la red doméstica argentina, de modo que muchos aeropuertos pequeños dependen de él, aunque esté poco conectado con los grandes hubs. El gráfico de posiciones por grado frente a intermediación confirma este patrón: las mayores discrepancias entre medidas no corresponden a los grandes hubs, sino a aeropuertos de tamaño medio cuya importancia está en su posición como puente.

La nube de puntos se distribuye en torno a la diagonal, lo que indica una asociación positiva entre ambas medidas (Spearman de 0,82): los aeropuertos con más destinos tienden a ocupar también posiciones altas en intermediación. Sin embargo, la dispersión es considerable y varios nodos se alejan de forma marcada de esta tendencia. Los casos más notables se ubican por debajo de la diagonal. Anchorage (ANC) ocupa el puesto 288 por grado, pero el 3 por intermediación. Los Ángeles (LAX) pasa del puesto 24 en grado al 2 en intermediación. Un patrón similar, aunque menos extremo, se observa en São Paulo (GRU), Brisbane (BNE), Papeete (PPT) y en un grupo de aeropuertos regionales de Alaska (BET, FAI) y del norte de Canadá (YWG, YZF, YQT).

6.3 Participación continental

Para distinguir hubs regionales de conectores intercontinentales se calcula el coeficiente de participación continental

\[ P_i = 1 - \sum_{c=1}^{C}\left(\frac{k_{i,c}}{k_i}\right)^2, \]

donde \(k_{i,c}\) es el número de vecinos de \(i\) en el continente \(c\) y \(k_i=\sum_c k_{i,c}\). Vale 0 cuando todos los vecinos están en un solo continente y crece cuando los vecinos se reparten entre varios; su máximo es \(1-1/C\) (0,8 con cinco continentes). Los vecinos sin continente asignado (“Otros”) se excluyen del conteo, y “Americas” es un único continente en la clasificación usada.

\(P_i\) tiene una limitación mecánica: con \(k_i\) vecinos, \(P_i\le 1-1/k_i\), de modo que un aeropuerto con pocas rutas no puede alcanzar valores altos. Por eso las comparaciones se hacen solo con aeropuertos de grado \(k_i\ge C\), y además se reporta la correlación parcial de Spearman controlando por el grado (correlación entre los residuos de los rangos de cada variable sobre el rango del grado).

Se compara \(P_i\) con la intermediación, la longitud media de las rutas incidentes \(\bar\ell_i\), el número de continentes alcanzados (continentes con al menos un vecino) y la coreness.

lev <- sort(setdiff(unique(cont), "Otros"))
C   <- length(lev)
ci  <- match(cont, lev)

a  <- c(el[, 1], el[, 2])
b  <- c(el[, 2], el[, 1])
ok <- !is.na(ci[b])
Kmat <- unclass(table(factor(a[ok], levels = seq_len(n_v)),
                      factor(ci[b[ok]], levels = seq_len(C))))
Kmat <- matrix(as.numeric(Kmat), nrow = n_v, ncol = C, dimnames = list(NULL, lev))

k_val <- rowSums(Kmat)
P     <- ifelse(k_val > 0, 1 - rowSums((Kmat / pmax(k_val, 1))^2), NA)

d_km       <- distHaversine(coords[el[, 1], ], coords[el[, 2], ]) / 1000
long_media <- as.numeric(strength(g, weights = d_km)) / df_med$grado

df_med$P          <- P
df_med$n_cont     <- rowSums(Kmat > 0)
df_med$long_media <- long_media

base <- df_med %>% filter(grado >= C, !is.na(P))

parcial <- function(x, y, z) {
  cor(resid(lm(rank(x) ~ rank(z))), resid(lm(rank(y) ~ rank(z))))
}

vars_cmp <- c("Grado (referencia)" = "grado", "Intermediación" = "intermediacion",
              "Longitud media de las rutas" = "long_media",
              "Continentes alcanzados" = "n_cont", "Coreness" = "coreness")

tab_cmp <- data.frame(
  Variable = names(vars_cmp),
  Spearman = sapply(vars_cmp, function(v) cor(base$P, base[[v]], method = "spearman")),
  `Spearman parcial (controlando grado)` = sapply(vars_cmp, function(v) {
    if (v == "grado") NA_real_ else parcial(base$P, base[[v]], base$grado)
  }),
  check.names = FALSE, row.names = NULL
)
kable(tab_cmp, digits = 3, format.args = list(decimal.mark = ","),
      caption = paste0("Asociación de P con otras medidas (aeropuertos con grado ≥ ", C,
                       "; n = ", nrow(base), ")."))
Asociación de P con otras medidas (aeropuertos con grado ≥ 5; n = 1239).
Variable Spearman Spearman parcial (controlando grado)
Grado (referencia) 0,475 —
Intermediación 0,412 0,132
Longitud media de las rutas 0,678 0,564
Continentes alcanzados 0,950 0,941
Coreness 0,455 0,062

6.4 Hubs regionales y conectores intercontinentales

Se consideran hubs los aeropuertos que superan el umbral de Tukey sobre el grado (\(k\ge Q_3+1{,}5\,\mathrm{IQR}\)), el mismo usado en las visualizaciones. Entre ellos, se clasifican como conectores intercontinentales los que tienen \(P_i\) por encima de la mediana de \(P\) entre los hubs, y como hubs regionales los demás. Es una clasificación relativa: el corte en la mediana es una convención, y el umbral P_cut puede modificarse para evaluar la sensibilidad. Lo que se examina es si los dos grupos difieren también en intermediación, longitud de las rutas, continentes alcanzados y coreness.

grado_v    <- df_med$grado
umbral_hub <- as.numeric(quantile(grado_v, 0.75) + 1.5 * IQR(grado_v))
P_cut      <- median(base$P[base$grado >= umbral_hub])

base <- base %>%
  mutate(tipo = factor(
    ifelse(grado >= umbral_hub,
           ifelse(P >= P_cut, "Conector intercontinental", "Hub regional"), "No hub"),
    levels = c("No hub", "Hub regional", "Conector intercontinental")))



res_tipo <- base %>%
  filter(tipo != "No hub") %>%
  group_by(tipo) %>%
  summarise(Aeropuertos = n(),
            `Grado (mediana)` = median(grado),
            `Intermediación (mediana)` = median(intermediacion),
            `Longitud media, km (mediana)` = median(long_media),
            `Continentes alcanzados (mediana)` = median(n_cont),
            `Coreness (mediana)` = median(coreness),
            `P (mediana)` = median(P),
            .groups = "drop")

kable(res_tipo %>% rename(Tipo = tipo), digits = c(0, 0, 1, 4, 0, 1, 1, 3),
      format.args = list(decimal.mark = ",", big.mark = "."),
      caption = "Hubs regionales vs. conectores intercontinentales.")
Hubs regionales vs. conectores intercontinentales.
Tipo Aeropuertos Grado (mediana) Intermediación (mediana) Longitud media, km (mediana) Continentes alcanzados (mediana) Coreness (mediana) P (mediana)
Hub regional 227 36 0,0005 1.266 1 24 0,000
Conector intercontinental 228 56 0,0029 2.137 3 25 0,322
hubs <- filter(base, tipo != "No hub")
kable(table(Tipo = droplevels(hubs$tipo), `Continente del aeropuerto` = hubs$continente),
      caption = "Distribución de los hubs por continente.")
Distribución de los hubs por continente.
Africa Americas Asia Europe Oceania
Hub regional 2 84 67 71 3
Conector intercontinental 27 33 63 94 11
top_tipo <- hubs %>%
  group_by(tipo) %>% arrange(desc(grado), .by_group = TRUE) %>% slice_head(n = 8) %>% ungroup() %>%
  select(tipo, iata, ciudad, pais, grado, P, n_cont, long_media, intermediacion, coreness)

kable(top_tipo, digits = c(0, 0, 0, 0, 0, 3, 0, 0, 4, 0),
      format.args = list(decimal.mark = ",", big.mark = "."),
      col.names = c("Tipo", "IATA", "Ciudad", "País", "Grado", "P", "Continentes",
                    "Longitud media (km)", "Intermediación", "Coreness"),
      caption = "Ocho hubs de mayor grado de cada tipo.")
Ocho hubs de mayor grado de cada tipo.
Tipo IATA Ciudad País Grado P Continentes Longitud media (km) Intermediación Coreness
Hub regional DFW Dallas-Fort Worth United States 188 0,102 4 2.038 0,0304 30
Hub regional DEN Denver United States 169 0,046 3 1.603 0,0264 30
Hub regional STN London United Kingdom 153 0,112 3 1.333 0,0050 31
Hub regional CLT Charlotte United States 141 0,082 2 1.489 0,0084 30
Hub regional DTW Detroit United States 135 0,127 3 1.612 0,0092 30
Hub regional LAS Las Vegas United States 133 0,087 3 2.254 0,0124 30
Hub regional MSP Minneapolis United States 131 0,060 3 1.459 0,0139 30
Hub regional PMI Palma de Mallorca Spain 127 0,031 2 1.304 0,0041 31
Conector intercontinental AMS Amsterdam Netherlands 248 0,593 4 3.361 0,0519 31
Conector intercontinental FRA Frankfurt Germany 244 0,648 4 3.681 0,0533 31
Conector intercontinental CDG Paris France 240 0,687 4 3.718 0,0646 31
Conector intercontinental ISL Istanbul Turkey 233 0,637 4 2.616 0,0431 31
Conector intercontinental ATL Atlanta United States 217 0,156 4 2.053 0,0304 30
Conector intercontinental ORD Chicago United States 206 0,237 3 2.295 0,0473 30
Conector intercontinental PEK Beijing China 205 0,403 5 3.178 0,0506 30
Conector intercontinental MUC Munich Germany 192 0,454 4 2.326 0,0165 31
col_tipo <- c("No hub" = "grey75", "Hub regional" = "#2b5c8f", "Conector intercontinental" = "#e41a1c")
lab_top  <- hubs %>% slice_max(intermediacion, n = 10)

p1 <- ggplot(filter(base, intermediacion > 0), aes(P, intermediacion, color = tipo)) +
  geom_point(alpha = 0.6, size = 1.4) +
  geom_vline(xintercept = P_cut, linetype = "dashed", color = "grey40") +
  geom_text_repel(data = lab_top, aes(label = iata), size = 3, color = "black", max.overlaps = Inf) +
  scale_y_log10() +
  scale_color_manual(values = col_tipo, name = NULL) +
  labs(title = "A. Participación continental vs. intermediación",
       x = "Participación continental P", y = "Intermediación normalizada (escala log)") +
  theme_minimal() + theme(panel.grid.minor = element_blank(),
                          plot.title = element_text(face = "bold", size = 11))

p2 <- ggplot(base, aes(P, long_media, color = tipo)) +
  geom_point(alpha = 0.6, size = 1.4) +
  geom_vline(xintercept = P_cut, linetype = "dashed", color = "grey40") +
  scale_y_log10(labels = comma) +
  scale_color_manual(values = col_tipo, name = NULL) +
  labs(title = "B. Participación continental vs. longitud media de las rutas",
       x = "Participación continental P", y = "Longitud media incidente (km, escala log)") +
  theme_minimal() + theme(panel.grid.minor = element_blank(),
                          plot.title = element_text(face = "bold", size = 11))

(p1 | p2) + plot_layout(guides = "collect") & theme(legend.position = "bottom")

Con el criterio de Tukey se consideran hubs los 455 aeropuertos con grado de al menos 20, y la mediana de \(P\) entre ellos, usada como corte, es 0,137 (línea punteada). Este corte separa a los hubs en dos grupos: hubs regionales, con conexiones concentradas en su propio continente, y conectores intercontinentales, con conexiones distribuidas entre varios continentes. Se trata de una clasificación relativa: hubs con valores de \(P\) cercanos al corte, como Dallas-Fort Worth (\(P = 0{,}10\)) y Atlanta (\(P = 0{,}16\)), quedan en grupos distintos a pesar de tener perfiles similares.

En el panel A, los aeropuertos con mayor intermediación corresponden casi en su totalidad a conectores intercontinentales, entre ellos Chicago (ORD), Los Ángeles (LAX), Pekín (PEK), Toronto (YYZ), Ámsterdam (AMS), Fráncfort (FRA), Estambul (ISL), Dubái (DXB) y París (CDG). Los hubs regionales presentan valores de intermediación intermedios, mayoritariamente entre \(10^{-4}\) y \(10^{-2}\). Anchorage (ANC) constituye la principal excepción: clasificado como hub regional, con \(P \approx 0\), alcanza una de las intermediaciones más altas de la red,

7 Cohesión y conectividad

Esta sección estudia qué tan “cerrada” es la red en el nivel local (tríadas, triángulos y agrupamiento) y qué tan frágil es su conectividad en el nivel global (puentes, puntos de articulación y conectividad). El censo triádico se calcula sobre la componente gigante dirigida (rutas origen → destino, sin rutas múltiples ni lazos), y el resto del análisis sobre la componente gigante no dirigida y sin pesos. El número de aerolíneas y la distancia en kilómetros se usan solo para caracterizar las aristas críticas. Todas las distancias en la red se miden en aristas (vuelos directos).

# ---- Red no dirigida (componente gigante) ------------------------------------
g36    <- giant_component
n36    <- vcount(g36)
m36    <- ecount(g36)
k36    <- degree(g36)
el36   <- as_edgelist(g36, names = FALSE)
iata36 <- ifelse(is.na(V(g36)$iata) | V(g36)$iata == "", V(g36)$icao, V(g36)$iata)
cont36 <- V(g36)$continent
cont36[is.na(cont36)] <- "Otros"

# Carpeta donde se guardan las simulaciones (solo se calculan la primera vez)
dir.create("cache_resultados", showWarnings = FALSE)

# ---- Red dirigida binaria y simple (mismos aeropuertos) ----------------------
g_dir_gc <- induced_subgraph(g_directed, which(cl$membership == largest_cluster_id))
g_dir_gc <- igraph::simplify(g_dir_gc, remove.multiple = TRUE, remove.loops = TRUE,
                             edge.attr.comb = "ignore")

# ---- Función común: comparación con un modelo nulo por simulación -------------
# obs: vector observado; sims: matriz B x length(obs) con las réplicas nulas.
# Devuelve media y desviación nulas, puntaje Z, perfil normalizado y valor p de
# Monte Carlo bilateral: p = min(1, 2 * min(p_superior, p_inferior)).
resumen_nulo <- function(obs, sims) {
  B_s  <- nrow(sims)
  mu   <- colMeans(sims)
  sdv  <- apply(sims, 2, sd)
  z    <- ifelse(sdv > 0, (obs - mu) / sdv, NA_real_)
  nz   <- sqrt(sum(z^2, na.rm = TRUE))
  zn   <- if (nz > 0) z / nz else z
  obs_mat <- matrix(obs, nrow = B_s, ncol = length(obs), byrow = TRUE)
  p_up <- (1 + colSums(sims >= obs_mat)) / (B_s + 1)
  p_lo <- (1 + colSums(sims <= obs_mat)) / (B_s + 1)
  data.frame(medida = names(obs), obs = obs, media_nula = mu, sd_nula = sdv,
             z = z, z_norm = zn, p_mc = pmin(1, 2 * pmin(p_up, p_lo)),
             row.names = NULL)
}

# Formato del valor p de Monte Carlo (su mínimo es 2 / (B + 1))
fmt_pmc <- function(p, B_s) {
  ifelse(p <= 2 / (B_s + 1) + 1e-12, paste0("≤ ", fmt_num(2 / (B_s + 1), 4)), fmt_num(p, 4))
}

7.1 Censo triádico y comparación con un modelo nulo

Una tríada es cualquier conjunto de tres aeropuertos. En la red dirigida, cada tríada pertenece a una de 16 clases según la configuración de sus arcos (clases MAN de Holland y Leinhardt, donde los dígitos cuentan díadas Mutuas, Asimétricas y Nulas).

Para saber qué configuraciones son más o menos frecuentes de lo esperado, el censo observado se compara con \(B\) redes simuladas que conservan exactamente el grado de salida y de entrada de cada aeropuerto (sample_degseq con el método edge.switching.simple). Para cada clase \(c\) se reportan:

  • el puntaje \(Z_c = (N_c^{obs} - \bar N_c)/s_c\), donde \(\bar N_c\) y \(s_c\) son la media y la desviación estándar bajo el modelo nulo;
  • el perfil normalizado \(\tilde Z_c = Z_c / \sqrt{\sum_{c'} Z_{c'}^2}\), que permite comparar el patrón de sobre y subrepresentación independientemente del tamaño de la red;
  • el valor \(p\) de Monte Carlo bilateral, \(p_c = \min\{1,\ 2\min(p_c^{+}, p_c^{-})\}\), con \(p_c^{+} = \frac{1 + \#\{b: N_c^{(b)} \ge N_c^{obs}\}}{B+1}\) y \(p_c^{-}\) análogo con \(\le\).

La regla \(|Z| \gtrsim 2\) se usa solo como referencia descriptiva, no como prueba formal; no se aplica al perfil normalizado, que está acotado entre \(-1\) y \(1\).

triad_labels <- c("003", "012", "102", "021D", "021U", "021C", "111D", "111U",
                  "030T", "030C", "201", "120D", "120U", "120C", "210", "300")

desc_man <- c(
  "003" = "Vacía", "012" = "Un arco asimétrico", "102" = "Un par mutuo",
  "021D" = "A ← B → C", "021U" = "A → B ← C", "021C" = "A → B → C",
  "111D" = "A ↔ B ← C", "111U" = "A ↔ B → C",
  "030T" = "Triángulo transitivo sin pares mutuos", "030C" = "Ciclo A → B → C → A",
  "201" = "Dos pares mutuos (tríada abierta)", "120D" = "A ← B → C, A ↔ C",
  "120U" = "A → B ← C, A ↔ C", "120C" = "A → B → C, A ↔ C",
  "210" = "A → B ↔ C, A ↔ C", "300" = "Tres pares mutuos"
)

tc_obs <- triad_census(g_dir_gc)
names(tc_obs) <- triad_labels

# ---- Modelo nulo: secuencia de grados de salida y entrada ---------------------
# 200 réplicas bastan para Z y valores p estables; súbalo solo si lo necesita.
B_tri       <- 200
archivo_tri <- file.path("cache_resultados", paste0("nulos_triadas_dirigidas_B", B_tri, ".rds"))

if (file.exists(archivo_tri)) {
  null_counts <- readRDS(archivo_tri)
} else {
  out_deg <- degree(g_dir_gc, mode = "out")
  in_deg  <- degree(g_dir_gc, mode = "in")
  null_counts <- matrix(NA_real_, nrow = B_tri, ncol = 16, dimnames = list(NULL, triad_labels))
  set.seed(123)
  for (b in seq_len(B_tri)) {
    g_b <- sample_degseq(out_deg, in_deg, method = "edge.switching.simple")
    null_counts[b, ] <- triad_census(g_b)
  }
  saveRDS(null_counts, archivo_tri)
}

res_tri <- resumen_nulo(tc_obs, null_counts)

kable(
  res_tri %>%
    transmute(Clase = medida, Configuración = desc_man[medida],
              Observado = fmt_num(obs, 0), `Media nula` = fmt_num(media_nula, 1),
              `DE nula` = fmt_num(sd_nula, 1), Z = fmt_num(z, 1),
              `Z normalizado` = fmt_num(z_norm, 3), `p (Monte Carlo)` = fmt_pmc(p_mc, B_tri)),
  row.names = FALSE, align = c("l", "l", rep("r", 6)),
  caption = paste0("Censo triádico dirigido frente a ", B_tri,
                   " redes con grados de salida y entrada preservados.")
)
Censo triádico dirigido frente a 200 redes con grados de salida y entrada preservados.
Clase Configuración Observado Media nula DE nula Z Z normalizado p (Monte Carlo)
003 Vacía 5.336.149.445 5.287.186.508,5 102.059,6 479,7 0,302 ≤ 0,0100
012 Un arco asimétrico 2.505.815 99.091.934,7 194.564,8 -496,4 -0,312 ≤ 0,0100
102 Un par mutuo 55.379.681 5.143.499,2 103.856,0 483,7 0,304 ≤ 0,0100
021D A ← B → C 523 723.028,3 5.084,5 -142,1 -0,089 ≤ 0,0100
021U A → B ← C 509 718.509,6 5.055,9 -142,0 -0,089 ≤ 0,0100
021C A → B → C 858 1.468.827,1 10.202,6 -143,9 -0,090 ≤ 0,0100
111D A ↔︎ B ← C 23.777 251.876,6 3.860,6 -59,1 -0,037 ≤ 0,0100
111U A ↔︎ B → C 29.484 253.097,5 3.873,6 -57,7 -0,036 ≤ 0,0100
030T Triángulo transitivo sin pares mutuos 13 84.667,8 1.405,1 -60,3 -0,038 ≤ 0,0100
030C Ciclo A → B → C → A 52 29.651,4 487,1 -60,8 -0,038 ≤ 0,0100
201 Dos pares mutuos (tríada abierta) 852.092 23.463,7 993,0 834,5 0,525 ≤ 0,0100
120D A ← B → C, A ↔︎ C 96 13.617,9 208,8 -64,8 -0,041 ≤ 0,0100
120U A → B ← C, A ↔︎ C 167 13.670,2 202,0 -66,8 -0,042 ≤ 0,0100
120C A → B → C, A ↔︎ C 251 28.097,3 365,6 -76,2 -0,048 ≤ 0,0100
210 A → B ↔︎ C, A ↔︎ C 3.837 11.296,7 409,9 -18,2 -0,011 ≤ 0,0100
300 Tres pares mutuos 96.236 1.089,6 93,6 1.016,0 0,639 ≤ 0,0100
par(mfrow = c(1, 2), mar = c(8, 4, 3, 1))

z_plot <- setNames(res_tri$z, res_tri$medida)
barplot(z_plot, las = 2, col = ifelse(z_plot > 0, "#2b5c8f", "#e41a1c"),
        main = "Puntajes Z por clase de tríada", ylab = "Z")
abline(h = c(-2, 0, 2), lty = c(2, 1, 2), col = 2)

zn_plot <- setNames(res_tri$z_norm, res_tri$medida)
barplot(zn_plot, las = 2, col = ifelse(zn_plot > 0, "#2b5c8f", "#e41a1c"),
        main = "Perfil normalizado", ylab = "Z normalizado", ylim = c(-1, 1))
abline(h = 0)

par(mfrow = c(1, 1))

# Resúmenes para el texto
top_sobre <- res_tri %>% filter(!is.na(z)) %>% arrange(desc(z)) %>% head(3)
top_bajo  <- res_tri %>% filter(!is.na(z)) %>% arrange(z) %>% head(3)
clases_asim <- c("012", "021D", "021U", "021C", "111D", "111U", "030T", "030C",
                 "120D", "120U", "120C", "210")
prop_asim <- sum(tc_obs[clases_asim]) / sum(tc_obs[names(tc_obs) != "003"])

El modelo nulo conserva cuántas rutas salen y llegan a cada aeropuerto, pero no si esas rutas son recíprocas ni con quién se conectan. Frente a esa referencia, las clases sobrerrepresentadas son 003, 102, 201 y 300, mientras que las que contienen al menos un arco asimétrico quedan por debajo de lo esperado.

Este patrón refleja, en primer lugar, la reciprocidad de la red. La clase 102 (un par mutuo y un tercer aeropuerto desconectado) muestra que las rutas se operan casi siempre en ambos sentidos: donde existe un vuelo \(A \to B\) suele existir también \(B \to A\). El modelo nulo, en cambio, asigna la ruta de ida y la de vuelta de forma independiente, por lo que genera muchos más arcos asimétricos. La sobrerrepresentación de 003 (tríada vacía) es una consecuencia de lo anterior: como los arcos se concentran en pares de ida y vuelta, se ocupan menos pares de aeropuertos de los que se ocuparían al azar, y quedan más tríadas sin ninguna conexión.

La clase 201 (dos pares mutuos sin cierre, \(A \leftrightarrow B \leftrightarrow C\)) corresponde a la estructura de hub y radios: un aeropuerto central conectado en ambos sentidos con dos destinos que no tienen ruta directa entre sí. Es la forma en que los hubs permiten viajar entre aeropuertos no conectados mediante una escala.

En segundo lugar, el perfil normalizado muestra un cierre triádico superior al esperado. El triángulo completamente recíproco (300) presenta el mayor exceso de todas las clases (Z = 1.016), lo que indica que los aeropuertos que comparten un vecino tienden a estar conectados entre sí más de lo que explican sus grados. Como en la red dirigida la reciprocidad y el cierre triádico se confunden (las clases 102, 201 y 300 se benefician de ambos), la siguiente subsección separa el segundo efecto usando la red no dirigida.

7.2 Transitividad

Sobre la red no dirigida se calculan tres medidas de cierre triádico con transitivity():

  • Transitividad global: \(T = 3\,t / \tau\), donde \(t\) es el número de triángulos y \(\tau = \sum_i \binom{k_i}{2}\) el número de tríadas conectadas (caminos de longitud dos).
  • Agrupamiento local: \(C_i = 2\,t_i / [k_i(k_i - 1)]\), donde \(t_i\) es el número de triángulos que contienen a \(i\); está definido solo para \(k_i \ge 2\).
  • Agrupamiento local promedio: \(\bar C = \frac{1}{|V_{k \ge 2}|}\sum_{i: k_i \ge 2} C_i\).

Para la comparación nula se repite el procedimiento anterior con \(B\) redes no dirigidas que conservan el grado de cada aeropuerto. A diferencia de la red dirigida, en este caso la reciprocidad no interviene, de modo que la comparación aísla el cierre triádico. Además, con los grados fijos, \(\tau\) es constante, así que la transitividad global depende únicamente del número de triángulos.

# ---- Medidas observadas (igraph) ---------------------------------------------
n_triang <- sum(count_triangles(g36)) / 3
T_glob   <- transitivity(g36, type = "global")
C_loc    <- transitivity(g36, type = "local", isolates = "NaN")
Cbar     <- transitivity(g36, type = "average")

p36 <- edge_density(g36)

# ---- Modelo nulo: secuencia de grados (red no dirigida) ----------------------
B_trans       <- 200
archivo_trans <- file.path("cache_resultados", paste0("nulos_transitividad_B", B_trans, ".rds"))

if (file.exists(archivo_trans)) {
  null_trans <- readRDS(archivo_trans)
} else {
  null_trans <- matrix(NA_real_, nrow = B_trans, ncol = 3,
                       dimnames = list(NULL, c("triangulos", "trans_global", "trans_promedio")))
  set.seed(123)
  for (b in seq_len(B_trans)) {
    g_b <- sample_degseq(k36, method = "edge.switching.simple")
    null_trans[b, ] <- c(sum(count_triangles(g_b)) / 3,
                         transitivity(g_b, type = "global"),
                         transitivity(g_b, type = "average"))
  }
  saveRDS(null_trans, archivo_trans)
}

obs_trans <- c(triangulos = n_triang, trans_global = T_glob, trans_promedio = Cbar)
res_trans <- resumen_nulo(obs_trans, null_trans)

kable(
  data.frame(
    Medida = c("Número de triángulos", "Transitividad global", "Agrupamiento local promedio (k ≥ 2)"),
    Observado = c(fmt_num(res_trans$obs[1], 0), fmt_num(res_trans$obs[2:3], 3)),
    `Media nula` = c(fmt_num(res_trans$media_nula[1], 1), fmt_num(res_trans$media_nula[2:3], 4)),
    `DE nula` = c(fmt_num(res_trans$sd_nula[1], 1), fmt_num(res_trans$sd_nula[2:3], 4)),
    Z = fmt_num(res_trans$z, 1),
    `p (Monte Carlo)` = fmt_pmc(res_trans$p_mc, B_trans),
    check.names = FALSE
  ),
  row.names = FALSE, align = c("l", rep("r", 5)),
  caption = paste0("Triángulos y transitividad de la red no dirigida frente a ", B_trans,
                   " redes con secuencia de grados preservada. Densidad de referencia: ",
                   fmt_num(p36, 4), ".")
)
Triángulos y transitividad de la red no dirigida frente a 200 redes con secuencia de grados preservada. Densidad de referencia: 0,0037.
Medida Observado Media nula DE nula Z p (Monte Carlo)
Número de triángulos 100.652 34.484,3 481,8 137,3 ≤ 0,0100
Transitividad global 0,250 0,0856 0,0012 137,3 ≤ 0,0100
Agrupamiento local promedio (k ≥ 2) 0,633 0,1166 0,0036 142,8 ≤ 0,0100

La red tiene 100.652 triángulos y una transitividad global de 0,250: el 25,0 % de las tríadas conectadas (caminos de dos aristas) está cerrado, a pesar de que la densidad es de apenas 0,37 %. El agrupamiento local promedio (0,633) es mayor que la transitividad global, porque ambas medidas ponderan de forma distinta. El promedio da el mismo peso a cada aeropuerto, y los aeropuertos pequeños, con pocos vecinos, suelen tenerlos conectados entre sí. La transitividad global, en cambio, cuenta tríadas, y la mayoría de las tríadas abiertas se concentran en los hubs, que conectan destinos sin relación directa entre sí.

Frente a las 200 redes simuladas con los mismos grados, la red observada tiene 2,9 veces más triángulos de lo esperado (Z = 137,3), una transitividad global 2,9 veces mayor (Z = 137,3) y un agrupamiento promedio 5,4 veces mayor (Z = 142,8). En los tres casos el valor observado supera a todas las réplicas (p ≤ 0,0100, el mínimo posible con 200 simulaciones). Por lo tanto, el exceso de cierre triádico no se explica por la distribución del grado. Es coherente con la restricción espacial: los aeropuertos cercanos tienden a compartir vecinos cercanos y, con ello, a cerrar triángulos, mientras que el modelo nulo, al ignorar la geografía, conecta aeropuertos distantes que rara vez tienen destinos en común.

7.3 Agrupamiento, grado y continentes

df_cl <- data.frame(k = k36, C = C_loc, cont = cont36) %>% filter(k >= 2)
      
# C(k): promedio del agrupamiento local por grado
ck <- df_cl %>% group_by(k) %>% summarise(C = mean(C), n = n(), .groups = "drop")
fit_ck  <- lm(log(C) ~ log(k), data = filter(ck, C > 0), weights = n)
beta_ck <- coef(fit_ck)[["log(k)"]]
ci_ck   <- confint(fit_ck)["log(k)", ]
rho_ck  <- cor(df_cl$k, df_cl$C, method = "spearman")

p_ck <- ggplot(df_cl, aes(k, C)) +
  geom_point(alpha = 0.15, color = "#2b5c8f", size = 1) +
  geom_point(data = filter(ck, C > 0), aes(k, C), color = "#e41a1c", size = 1.8) +
  geom_abline(intercept = coef(fit_ck)[1] / log(10), slope = beta_ck,
              linetype = "dashed", color = "black") +
  scale_x_log10() + scale_y_log10() +
  labs(title = "A. Agrupamiento local vs. grado",
       subtitle = "Rojo: promedio C(k) por grado | Línea: ajuste log-log ponderado",
       x = "Grado k (escala log)", y = "Coeficiente de agrupamiento local C (escala log)") +
  theme_minimal() + theme(panel.grid.minor = element_blank(),
                          plot.title = element_text(face = "bold", size = 11))

p_ck_cont <- ggplot(filter(df_cl, cont != "Otros"), aes(k, C, color = cont)) +
  geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"), se = FALSE, linewidth = 1) +
  scale_x_log10() +
  scale_color_brewer(palette = "Set1", name = "Continente") +
  labs(title = "B. C(k) suavizado por continente",
       x = "Grado k (escala log)", y = "Coeficiente de agrupamiento local C") +
  theme_minimal() + theme(panel.grid.minor = element_blank(),
                          plot.title = element_text(face = "bold", size = 11))

p_ck | p_ck_cont

# Triángulos según cuántos continentes tocan
tri <- matrix(as.integer(triangles(g36)), ncol = 3, byrow = TRUE)
c1 <- cont36[tri[, 1]]; c2 <- cont36[tri[, 2]]; c3 <- cont36[tri[, 3]]
n_cont_tri    <- 1 + (c2 != c1) + (c3 != c1 & c3 != c2)
prop_tri_mono <- mean(n_cont_tri == 1)

conts36 <- sort(setdiff(unique(cont36), "Otros"))
tab_cont36 <- do.call(rbind, lapply(conts36, function(cc) {
  idx  <- which(cont36 == cc)
  sub  <- induced_subgraph(g36, idx)
  toca <- c1 == cc | c2 == cc | c3 == cc
  data.frame(
    Continente = cc,
    Aeropuertos = length(idx),
    `Grado medio` = mean(k36[idx]),
    `C local promedio (k ≥ 2)` = mean(C_loc[idx], na.rm = TRUE),
    `Transitividad de la subred intracontinental` = transitivity(sub, type = "global"),
    `Triángulos que lo tocan` = sum(toca),
    `% de esos triángulos dentro del continente` = 100 * mean((c1 == cc & c2 == cc & c3 == cc)[toca]),
    check.names = FALSE
  )
}))

kable(tab_cont36, digits = c(0, 0, 1, 3, 3, 0, 1), row.names = FALSE,
      format.args = list(decimal.mark = ",", big.mark = "."),
      caption = "Agrupamiento y triángulos por continente.")
Agrupamiento y triángulos por continente.
Continente Aeropuertos Grado medio C local promedio (k ≥ 2) Transitividad de la subred intracontinental Triángulos que lo tocan % de esos triángulos dentro del continente
Africa 288 7,7 0,682 0,224 6.487 12,2
Americas 1.228 9,1 0,634 0,261 31.033 75,3
Asia 796 12,7 0,683 0,298 34.119 53,2
Europe 603 21,4 0,520 0,317 54.526 61,0
Oceania 266 4,7 0,697 0,185 1.333 35,3

El panel A muestra que el agrupamiento local disminuye con el grado: el ajuste \(C(k) \propto k^{\beta}\) da \(\hat\beta = -0,26\) (IC 95 %: [-0,28; -0,25]), y la correlación de Spearman entre \(C_i\) y \(k_i\) es de -0,50. Los aeropuertos con pocos destinos presentan agrupamientos cercanos a 0,7–0,8, es decir, sus vecinos suelen estar conectados entre sí, mientras que en los hubs el agrupamiento cae a valores de 0,1–0,2, porque conectan muchos destinos sin ruta directa entre ellos y la mayoría de sus tríadas quedan abiertas. Este decaimiento es característico de una organización jerárquica, con grupos regionales muy cohesionados unidos por hubs que actúan como intermediarios, y es el mismo mecanismo que hace que la transitividad global sea menor que el agrupamiento promedio.

Los triángulos son, en su mayoría, regionales: el 75,6 % de ellos tiene sus tres aeropuertos en el mismo continente. El panel B y la tabla por continente muestran que el patrón de decaimiento se repite en todas las regiones, aunque con intensidades distintas. Como \(C\) depende del grado, el agrupamiento promedio de cada continente refleja en buena medida su distribución de grados: Oceanía y África tienen los promedios más altos (0,697 y 0,682), pero también los grados medios más bajos (4,7 y 7,7), mientras que Europa tiene el menor promedio (0,520) y el mayor grado medio (21,4). Al comparar a igual grado, Europa presenta la curva más plana y la mayor transitividad intracontinental (0,317), en coherencia con el núcleo denso de hubs europeos identificado en la caracterización de los aeropuertos. Oceanía y África, en cambio, pierden agrupamiento con rapidez y tienen la menor transitividad intracontinental (0,185 y 0,224), lo que refleja estructuras en estrella. Además, solo el 12,2 % de los triángulos que involucran un aeropuerto africano y el 35,3 % de los de Oceanía se cierran dentro del propio continente, frente al 75,3 % de América: los aeropuertos africanos y de Oceanía forman triángulos principalmente con hubs de otros continentes, lo que confirma su dependencia funcional de Europa y Asia. Las curvas de los extremos deben leerse con cautela, ya que en esos rangos de grado hay pocos aeropuertos por continente.

7.4 Clanes

Un clan (clique) es un conjunto de aeropuertos en el que todos tienen ruta directa con todos. Un clan es maximal si no puede ampliarse con otro aeropuerto, y máximo si no existe uno de mayor tamaño; el tamaño de los clanes máximos es el número clan \(\omega(G)\). No se enumeran todos los clanes maximales, porque en una red de este tamaño son demasiados y su cálculo es muy costoso; se reportan solo el número clan y los clanes máximos.

w_clique <- clique_num(g36)
lc       <- largest_cliques(g36)                      # clanes máximos

comp_lc <- vapply(lc, function(v) length(unique(cont36[v])), numeric(1))
ej_lc   <- as.integer(lc[[1]])

kable(
  data.frame(Aeropuerto = iata36[ej_lc], Ciudad = V(g36)$city[ej_lc],
             País = V(g36)$country[ej_lc], Continente = cont36[ej_lc],
             Grado = k36[ej_lc]) %>% arrange(desc(Grado)),
  row.names = FALSE,
  caption = paste0("Uno de los ", length(lc), " clanes máximos (", w_clique, " aeropuertos todos conectados entre sí).")
)
Uno de los 2 clanes máximos (24 aeropuertos todos conectados entre sí).
Aeropuerto Ciudad País Continente Grado
ATL Atlanta United States Americas 217
ORD Chicago United States Americas 206
DFW Dallas-Fort Worth United States Americas 188
DEN Denver United States Americas 169
IAH Houston United States Americas 169
EWR Newark United States Americas 153
LAX Los Angeles United States Americas 148
YYZ Toronto Canada Americas 147
CLT Charlotte United States Americas 141
DTW Detroit United States Americas 135
LAS Las Vegas United States Americas 133
MSP Minneapolis United States Americas 131
PHL Philadelphia United States Americas 123
SFO San Francisco United States Americas 105
BOS Boston United States Americas 103
MCO Orlando United States Americas 102
FLL Fort Lauderdale United States Americas 97
DCA Washington United States Americas 91
PHX Phoenix United States Americas 90
STL St. Louis United States Americas 61
CLE Cleveland United States Americas 56
MCI Kansas City United States Americas 43
MSY New Orleans United States Americas 39
MKE Milwaukee United States Americas 30

El número clan de la red es \(\omega = 24\): existen 2 conjuntos de 24 aeropuertos (que difieren en un solo aeropuerto) con ruta directa entre todos ellos, formados íntegramente por hubs de Estados Unidos y Canadá. Este resultado complementa el de la coreness: el núcleo más denso de la red (coreness 31) lo forman los hubs europeos, mientras que el mayor conjunto de aeropuertos completamente interconectados corresponde a la red doméstica norteamericana, donde los grandes hubs (coreness 30) mantienen rutas directas entre todos ellos.

7.5 Componentes, puentes y puntos de articulación

En la red dirigida se distinguen las componentes débiles (conectadas ignorando la dirección) de las fuertes (donde cada aeropuerto puede llegar a cualquier otro y regresar respetando el sentido de las rutas). En la red no dirigida, un puente es una arista cuya eliminación desconecta la red y un punto de articulación es un aeropuerto cuya eliminación lo hace. Para cada uno se mide cuántos aeropuertos quedarían separados de la componente principal.

comp_debil  <- components(g_dir_gc, mode = "weak")
comp_fuerte <- components(g_dir_gc, mode = "strong")

kable(
  data.frame(
    Medida = c("Componentes débiles", "Componentes fuertes",
               "Aeropuertos en la mayor componente fuerte",
               "Proporción en la mayor componente fuerte"),
    Valor = c(comp_debil$no, fmt_num(comp_fuerte$no, 0),
              fmt_num(max(comp_fuerte$csize), 0),
              fmt_por(max(comp_fuerte$csize) / n36, 1))
  ),
  row.names = FALSE, align = c("l", "r"),
  caption = "Componentes débiles y fuertes de la componente gigante dirigida."
)
Componentes débiles y fuertes de la componente gigante dirigida.
Medida Valor
Componentes débiles 1
Componentes fuertes 42
Aeropuertos en la mayor componente fuerte 3.147
Proporción en la mayor componente fuerte 98,7 %
br_ids <- as.integer(bridges(g36))

# Tamaño del lado menor al eliminar cada puente
lado_menor <- vapply(br_ids, function(e) {
  min(components(delete_edges(g36, e))$csize)
}, numeric(1))

n_br         <- length(br_ids)
n_br_trivial <- sum(lado_menor == 1)

df_br <- data.frame(
  id = br_ids, a = el36[br_ids, 1], b = el36[br_ids, 2], separa = lado_menor,
  km = E(g36)$distance_km[br_ids], aerolineas = E(g36)$num_aerolineas[br_ids]
) %>%
  mutate(Ruta = paste(iata36[a], "–", iata36[b]),
         Ciudades = paste(V(g36)$city[a], "–", V(g36)$city[b])) %>%
  filter(separa > 1) %>%
  arrange(desc(separa))

kable(
  head(df_br, 10) %>%
    transmute(Ruta, Ciudades, `Aeropuertos separados` = separa,
              `Distancia (km)` = km, Aerolíneas = aerolineas),
  digits = 0, row.names = FALSE,
  format.args = list(decimal.mark = ",", big.mark = "."),
  caption = "Puentes no triviales que separan más aeropuertos de la componente principal."
)
Puentes no triviales que separan más aeropuertos de la componente principal.
Ruta Ciudades Aeropuertos separados Distancia (km) Aerolíneas
MAO – TFF Manaus – Tefe 5 521 1
YMO – YFA Moosonee – Fort Albany 4 126 1
ADQ – ANC Kodiak – Anchorage 4 406 2
ANC – CDB Anchorage – Cold Bay 4 997 2
BSB – MQH Brasilia – Minacu 4 260 1
GOH – RKV Godthaab – Reykjavik 4 1.430 1
YFA – ZKE Fort Albany – Kashechewan 3 9 1
TLV – ETH Tel-aviv – Elat 3 273 2
MAB – RDC Maraba – Redencao 3 311 1
RGN – TVY Yangon – Dawei 3 383 1
ap <- as.integer(articulation_points(g36))

separados_ap <- vapply(ap, function(v) {
  cc <- components(delete_vertices(g36, v))$csize
  sum(cc) - max(cc)
}, numeric(1))

df_ap <- data.frame(
  v = ap, separa = separados_ap,
  Aeropuerto = iata36[ap], Ciudad = V(g36)$city[ap], País = V(g36)$country[ap],
  Grado = k36[ap], intermediacion = df_med$intermediacion[ap]
) %>% arrange(desc(separa))

kable(
  head(df_ap, 12) %>%
    transmute(Aeropuerto, Ciudad, País, Grado,
              `Intermediación (%)` = 100 * intermediacion,
              `Aeropuertos que quedan aislados` = separa),
  digits = c(0, 0, 0, 0, 2, 0), row.names = FALSE,
  format.args = list(decimal.mark = ",", big.mark = "."),
  caption = "Puntos de articulación cuya eliminación separa más aeropuertos de la componente principal."
)
Puntos de articulación cuya eliminación separa más aeropuertos de la componente principal.
Aeropuerto Ciudad País Grado Intermediación (%) Aeropuertos que quedan aislados
ANC Anchorage United States 34 5,94 87
PPT Papeete French Polynesia 30 1,64 26
BOG Bogota Colombia 74 2,44 24
YZF Yellowknife Canada 18 1,21 19
BET Bethel United States 22 1,38 18
DEN Denver United States 169 2,64 18
SJO San Jose Costa Rica 38 1,06 16
AEP Buenos Aires Argentina 44 1,51 15
HIR Honiara Solomon Islands 19 0,94 15
FAI Fairbanks United States 18 1,22 15
HND Tokyo Japan 74 1,69 14
FIH Kinshasa Congo (Kinshasa) 20 0,99 13
rho_ap <- cor(df_ap$separa, df_ap$Grado, method = "spearman")

En la red dirigida, la componente gigante es débilmente conexa por construcción, y su mayor componente fuerte reúne el 98,7 % de los aeropuertos; los restantes solo pueden alcanzarse o abandonarse en un sentido porque algunas rutas solo se operan en un sentido.

La red no dirigida tiene 746 puentes (4,0 % de las aristas), pero 705 de ellos (94,5 %) son triviales: conectan un aeropuerto de grado 1 con el resto de la red, así que su pérdida solo aísla a ese aeropuerto. Los 41 puentes restantes sostienen a grupos de dos o más aeropuertos; el más crítico separa 5 aeropuertos (MAO – TFF, Manaus – Tefe).

Hay 329 puntos de articulación (10,3 % de los aeropuertos). Su número de aeropuertos dependientes se correlaciona moderadamente con el grado (Spearman de 0,32). Los casos con más aeropuertos dependientes (ANC, PPT, BOG) no son necesariamente los hubs de mayor grado, sino puertas de entrada a subredes regionales sin otra conexión con el resto del mundo, en concordancia con lo observado al comparar el grado con la intermediación.

7.6 Conectividad de vértices y aristas

La conectividad de vértices \(\kappa(G)\) y de aristas \(\lambda(G)\) son el mínimo número de aeropuertos o de rutas que hay que eliminar para desconectar la red. Por el teorema de Menger, coinciden con el mínimo número de caminos disjuntos entre pares de nodos, y siempre cumplen \(\kappa(G) \le \lambda(G) \le \delta(G)\), donde \(\delta(G)\) es el grado mínimo. Como la componente gigante tiene aeropuertos de grado 1 y puntos de articulación, ambas valen 1. Para obtener una medida informativa se calculan (i) la conectividad del núcleo más denso (el \(k\)-core de mayor orden) y (ii) el número de rutas disjuntas entre pares de los diez aeropuertos de mayor grado.

stopifnot(is_connected(g36))
kappa_g  <- vertex_connectivity(g36)
lambda_g <- edge_connectivity(g36)
delta_g  <- min(k36)

core_max_id <- which(coreness(g36) == max(coreness(g36)))
g_core      <- induced_subgraph(g36, core_max_id)
kappa_core  <- vertex_connectivity(g_core)
lambda_core <- edge_connectivity(g_core)
delta_core  <- min(degree(g_core))

# Rutas disjuntas en aristas entre los 10 aeropuertos de mayor grado
hubs10   <- order(-k36)[1:10]
pares_h  <- combn(hubs10, 2)
lambda_h <- apply(pares_h, 2, function(p) edge_connectivity(g36, source = p[1], target = p[2]))
cota_h   <- apply(pares_h, 2, function(p) min(k36[p]))

kable(
  data.frame(
    Red = c("Componente gigante",
            paste0("Núcleo de mayor coreness (", length(core_max_id), " aeropuertos)")),
    `Conectividad de vértices κ` = c(kappa_g, kappa_core),
    `Conectividad de aristas λ` = c(lambda_g, lambda_core),
    `Grado mínimo δ` = c(delta_g, delta_core),
    check.names = FALSE
  ),
  row.names = FALSE, align = c("l", "r", "r", "r"),
  caption = "Conectividad de vértices y aristas."
)
Conectividad de vértices y aristas.
Red Conectividad de vértices κ Conectividad de aristas λ Grado mínimo δ
Componente gigante 1 1 1
Núcleo de mayor coreness (93 aeropuertos) 31 31 31

En la componente gigante, \(\kappa = \lambda = 1\), como era de esperar. En el núcleo de mayor coreness (93 aeropuertos), en cambio, hay que eliminar 31 aeropuertos o 31 rutas para desconectarlo, frente a un grado mínimo de 31 dentro del núcleo. Entre los diez aeropuertos de mayor grado, el número de rutas disjuntas va de 177 a 244, y en el 29 % de los pares alcanza su máximo posible (el grado del menor de los dos). La red es, por tanto, muy robusta entre sus hubs y frágil en su periferia: la vulnerabilidad se concentra en los aeropuertos que dependen de una sola conexión o de un único punto de entrada.

7.7 Aristas críticas: intermediación, distancia y redundancia operacional

La intermediación de una arista es la fracción de caminos más cortos (en aristas) entre todos los pares de aeropuertos que pasan por ella. Se consideran de alta intermediación las aristas del 1 % superior, y estructuralmente críticas con baja redundancia operacional aquellas que, además, son operadas por una sola aerolínea.

eb <- edge_betweenness(g36, directed = FALSE, weights = NA)

df_eb <- data.frame(
  id = seq_len(m36), a = el36[, 1], b = el36[, 2],
  eb = eb, eb_pct = 100 * eb / choose(n36, 2),
  km = E(g36)$distance_km,
  aerolineas = E(g36)$num_aerolineas,
  registros = E(g36)$num_registros
) %>%
  mutate(
    puente = id %in% br_ids,
    intercont = cont36[a] != cont36[b] & cont36[a] != "Otros" & cont36[b] != "Otros",
    Ruta = paste(iata36[a], "–", iata36[b]),
    Ciudades = paste(V(g36)$city[a], "–", V(g36)$city[b]),
    aer_cat = factor(ifelse(aerolineas >= 6, "6+", as.character(aerolineas)),
                     levels = c(as.character(1:5), "6+"))
  )

umbral_eb <- quantile(df_eb$eb, 0.99)
df_eb$alta <- df_eb$eb >= umbral_eb
criticas   <- df_eb %>% filter(alta, aerolineas == 1) %>% arrange(desc(eb))

rho_eb_km  <- cor(df_eb$eb, df_eb$km, method = "spearman")
rho_eb_aer <- cor(df_eb$eb, df_eb$aerolineas, method = "spearman")

kable(
  data.frame(
    Grupo = c("Alta intermediación (1 % superior)", "Resto de aristas"),
    Aristas = c(sum(df_eb$alta), sum(!df_eb$alta)),
    `Distancia mediana (km)` = c(median(df_eb$km[df_eb$alta]), median(df_eb$km[!df_eb$alta])),
    `Aerolíneas (mediana)` = c(median(df_eb$aerolineas[df_eb$alta]), median(df_eb$aerolineas[!df_eb$alta])),
    `% con una sola aerolínea` = c(100 * mean(df_eb$aerolineas[df_eb$alta] == 1),
                                   100 * mean(df_eb$aerolineas[!df_eb$alta] == 1)),
    `% intercontinentales` = c(100 * mean(df_eb$intercont[df_eb$alta]),
                               100 * mean(df_eb$intercont[!df_eb$alta])),
    `% puentes` = c(100 * mean(df_eb$puente[df_eb$alta]), 100 * mean(df_eb$puente[!df_eb$alta])),
    check.names = FALSE
  ),
  digits = c(0, 0, 0, 1, 1, 1, 1), row.names = FALSE,
  format.args = list(decimal.mark = ",", big.mark = "."),
  caption = "Características de las aristas de alta intermediación frente al resto."
)
Características de las aristas de alta intermediación frente al resto.
Grupo Aristas Distancia mediana (km) Aerolíneas (mediana) % con una sola aerolínea % intercontinentales % puentes
Alta intermediación (1 % superior) 189 4.047 2 42,3 46,0 8,5
Resto de aristas 18.644 1.153 1 55,7 13,9 3,9
kable(
  head(criticas, 12) %>%
    transmute(Ruta, Ciudades, `Intermediación (% de pares)` = eb_pct,
              `Distancia (km)` = km,
              Intercontinental = ifelse(intercont, "Sí", "No"),
              Puente = ifelse(puente, "Sí", "No")),
  digits = c(0, 0, 2, 0, 0, 0, 0), row.names = FALSE,
  format.args = list(decimal.mark = ",", big.mark = "."),
  caption = "Aristas críticas: alta intermediación y operadas por una sola aerolínea."
)
Aristas críticas: alta intermediación y operadas por una sola aerolínea.
Ruta Ciudades Intermediación (% de pares) Distancia (km) Intercontinental Puente
SEA – FAI Seattle – Fairbanks 1,17 2.464 No No
OTZ – ANC Kotzebue – Anchorage 0,75 881 No No
OME – ANC Nome – Anchorage 0,75 866 No No
SFJ – CPH Sondrestrom – Copenhagen 0,73 3.433 Sí No
YYC – YZF Calgary – Yellowknife 0,67 1.264 No No
YUL – YVP Montreal – Quujjuaq 0,62 1.451 No No
GOH – KEF Godthaab – Keflavik 0,58 1.402 Sí No
HNL – ANC Honolulu – Anchorage 0,57 4.479 No No
CDG – BSB Paris – Brasilia 0,56 8.766 Sí No
JNU – SEA Juneau – Seattle 0,53 1.462 No No
POM – NRT Port Moresby – Tokyo 0,52 5.083 Sí No
YWG – YRL Winnipeg – Red Lake 0,45 276 No No
p_eb_aer <- ggplot(df_eb, aes(x = aer_cat, y = eb_pct)) +
  geom_boxplot(fill = "#2b5c8f", alpha = 0.5, outlier.alpha = 0.2) +
  scale_y_log10() +
  labs(title = "A. Intermediación según número de aerolíneas",
       x = "Aerolíneas que operan la ruta", y = "Intermediación de la arista (% de pares, log)") +
  theme_minimal() + theme(panel.grid.minor = element_blank(),
                          plot.title = element_text(face = "bold", size = 11))

p_eb_km <- ggplot(df_eb, aes(km, eb_pct)) +
  geom_point(aes(color = alta & aerolineas == 1), alpha = 0.4, size = 1) +
  scale_color_manual(values = c("FALSE" = "grey60", "TRUE" = "#e41a1c"),
                     labels = c("Resto", "Crítica (alta intermediación, 1 aerolínea)"), name = NULL) +
  scale_x_log10(labels = comma) + scale_y_log10() +
  labs(title = "B. Intermediación según distancia",
       x = "Distancia de la ruta (km, escala log)", y = "Intermediación de la arista (% de pares, log)") +
  theme_minimal() + theme(panel.grid.minor = element_blank(), legend.position = "bottom",
                          plot.title = element_text(face = "bold", size = 11))

p_eb_aer | p_eb_km

Las aristas del 1 % superior de intermediación son más largas que el resto (mediana de 4.047 km frente a 1.153 km), y el 46,0 % de ellas es intercontinental, frente al 13,9 % del resto. En el conjunto de la red, sin embargo, la intermediación de una arista apenas se relaciona con su distancia (Spearman de 0,15) y prácticamente no se relaciona con su número de aerolíneas (0,06).

En promedio, las aristas de alta intermediación tienen más redundancia operacional que el resto (el 42,3 % son operadas por una sola aerolínea, frente al 55,7 % del resto), pero esa redundancia no siempre acompaña a la importancia estructural: 80 de las 189 aristas de alta intermediación (42,3 %) son operadas por una sola aerolínea, y 13 de ellas son además puentes. Estas rutas son los enlaces más frágiles de la red: concentran una fracción importante de los caminos más cortos y su continuidad depende de un único operador. La más crítica es SEA – FAI (Seattle – Fairbanks), por la que pasa el 1,17 % de los caminos más cortos entre pares de aeropuertos.

8 Detección de comunidades

Se aplican dos algoritmos de optimización de la modularidad sobre la componente gigante. El algoritmo de Louvain agrupa los aeropuertos de forma jerárquica: primero mueve cada nodo a la comunidad vecina que más aumenta la modularidad y luego agrega las comunidades en supernodos, repitiendo el proceso hasta que la modularidad deja de aumentar. El algoritmo de Leiden sigue la misma lógica, pero añade una fase de refinamiento que garantiza que las comunidades estén bien conectadas internamente, algo que Louvain no asegura. Ambos se aplican primero sobre la red binaria, en la que todas las rutas pesan lo mismo, y después sobre la red ponderada por el número de aerolíneas que operan cada ruta, para evaluar cómo cambia la partición al incorporar la intensidad de servicio. Como ambos algoritmos son estocásticos, se fija una semilla antes de cada ejecución.

La calidad de cada partición se resume con la modularidad,

\[ Q = \frac{1}{2m}\sum_{i,j}\left(y_{ij} - \frac{k_i k_j}{2m}\right)\delta(c_i, c_j), \]

que compara la proporción de aristas dentro de las comunidades con la esperada si las conexiones se distribuyeran al azar conservando los grados; valores mayores indican comunidades más separadas entre sí. En la versión ponderada, \(y_{ij}\) y \(k_i\) se reemplazan por el número de aerolíneas de la ruta y la fuerza del aeropuerto. Como referencia se usa la partición por continentes.

8.1 Particiones de la red binaria

# Red binaria: se elimina el atributo de peso para que ningún algoritmo lo use
g_com <- delete_edge_attr(giant_component, "weight")
cont  <- V(g_com)$continent; cont[is.na(cont)] <- "Otros"

# Particiones (cada algoritmo con su propia semilla)
set.seed(123); m_louv <- membership(cluster_louvain(g_com))
set.seed(123); m_leid <- membership(cluster_leiden(g_com, objective_function = "modularity",
                                                   n_iterations = -1))

part <- list(Continentes = as.integer(factor(cont)), Louvain = m_louv, Leiden = m_leid)

# Mapas con layout geográfico (longitud, latitud); se colorean las 10 comunidades
# más grandes y el resto queda en gris. Los colores no se corresponden entre paneles.
par(mfrow = c(3, 1), mar = c(0.5, 0.5, 2.5, 0.5))
for (nm in names(part)) {
  m   <- as.integer(part[[nm]])
  top <- as.integer(names(sort(table(m), decreasing = TRUE)))[1:min(10, length(unique(m)))]
  col <- ifelse(m %in% top, hcl.colors(length(top), "Dynamic")[match(m, top)], "grey80")
  plot(g_com, layout = cbind(V(g_com)$longitude, V(g_com)$latitude), asp = 0,
       vertex.label = NA, vertex.size = 1.5, vertex.color = col, vertex.frame.color = NA,
       edge.color = adjustcolor("grey40", 0.05), edge.width = 0.2,
       main = paste0(nm, " (K = ", length(unique(m)), "; Q = ",
                     fmt_num(modularity(g_com, m), 3), ")"))
}

par(mfrow = c(1, 1))

Como los mapas solo permiten distinguir las comunidades más grandes, la composición de cada una se resume en la siguiente tabla: número de aeropuertos, continente predominante, países principales y aeropuerto de mayor grado (hub principal).

# Resume cada comunidad de una partición
describe_part <- function(m, nombre) {
  df <- data.frame(com = as.integer(m), pais = V(g_com)$country, cont = cont,
                   iata = V(g_com)$iata, grado = degree(g_com))
  do.call(rbind, lapply(split(df, df$com), function(d) {
    paises <- sort(table(d$pais), decreasing = TRUE)
    data.frame(
      Algoritmo = nombre,
      Comunidad = d$com[1],
      Aeropuertos = nrow(d),
      `Continente principal` = names(which.max(table(d$cont))),
      `% del continente principal` = 100 * max(table(d$cont)) / nrow(d),
      `Países principales` = paste0(names(paises)[1:min(3, length(paises))],
                                    " (", paises[1:min(3, length(paises))], ")", collapse = ", "),
      `Hub principal` = d$iata[which.max(d$grado)],
      check.names = FALSE
    )
  })) %>% arrange(desc(Aeropuertos))
}

cols <- c("Comunidad", "Aeropuertos", "Continente principal",
          "% del continente principal", "Países principales", "Hub principal")

tab_louv <- describe_part(m_louv, "Louvain")[, cols] %>% arrange(desc(Comunidad))
tab_leid <- describe_part(m_leid, "Leiden")[, cols]  %>% arrange(desc(Comunidad))

knitr::kables(
  list(
    kable(tab_louv, digits = 1, row.names = FALSE, valign = "t",
          format.args = list(decimal.mark = ",", big.mark = "."), caption = "Louvain"),
    kable(tab_leid, digits = 1, row.names = FALSE, valign = "t",
          format.args = list(decimal.mark = ",", big.mark = "."), caption = "Leiden")
  ),
  caption = "Composición de las comunidades detectadas por Louvain y Leiden (red binaria)."
)
Composición de las comunidades detectadas por Louvain y Leiden (red binaria).
Louvain
Comunidad Aeropuertos Continente principal % del continente principal Países principales Hub principal
17 131 Americas 100,0 United States (131) ANC
16 25 Americas 100,0 Bolivia (14), Chile (6), Peru (5) LPB
15 606 Asia 79,5 China (172), Australia (113), Indonesia (63) PEK
14 7 Americas 85,7 Panama (6), Turkey (1) PAC
13 3 Asia 100,0 Israel (3) ETH
12 161 Europe 72,0 Russia (103), Kazakhstan (17), Uzbekistan (11) DME
11 10 Africa 100,0 Kenya (10) WIL
10 26 Europe 100,0 Norway (26) TOS
9 160 Americas 100,0 Brazil (122), Argentina (31), Chile (4) GRU
8 7 Europe 100,0 United Kingdom (7) KOI
7 506 Asia 50,4 India (68), Iran (38), Turkey (36) ISL
6 751 Americas 99,5 United States (410), Canada (63), Mexico (56) ATL
5 47 Americas 100,0 Canada (47) YZF
4 93 Americas 100,0 Canada (93) YQB
3 486 Europe 92,2 France (54), United Kingdom (45), Spain (40) AMS
2 21 Americas 81,0 Greenland (17), Iceland (4) GOH
1 148 Oceania 97,3 French Polynesia (27), Vanuatu (27), New Zealand (25) AKL
Leiden
Comunidad Aeropuertos Continente principal % del continente principal Países principales Hub principal
17 131 Americas 100,0 United States (131) ANC
16 488 Asia 98,8 China (172), Indonesia (63), Japan (62) PEK
15 28 Oceania 96,4 French Polynesia (27), Chile (1) PPT
14 7 Americas 85,7 Panama (6), Turkey (1) PAC
13 3 Asia 100,0 Israel (3) ETH
12 161 Europe 72,0 Russia (103), Kazakhstan (17), Uzbekistan (11) DME
11 28 Europe 100,0 Norway (28) TOS
10 298 Americas 100,0 Brazil (122), Colombia (48), Argentina (36) GRU
9 7 Europe 100,0 United Kingdom (7) KOI
8 281 Asia 90,7 India (68), Iran (38), Turkey (36) ISL
7 214 Africa 100,0 Congo (Kinshasa) (18), South Africa (18), Kenya (16) CMN
6 638 Americas 99,4 United States (410), Canada (63), Mexico (56) ATL
5 47 Americas 100,0 Canada (47) YZF
4 93 Americas 100,0 Canada (93) YQB
3 505 Europe 88,3 France (54), United Kingdom (45), Spain (40) AMS
2 21 Americas 81,0 Greenland (17), Iceland (4) GOH
1 238 Oceania 98,3 Australia (113), Vanuatu (27), New Zealand (25) SYD

Louvain y Leiden identifican el mismo número de comunidades (17) y modularidades similares (0,653 y 0,660), ambas superiores a la de la partición por continentes (0,574), y coinciden en buena parte de su estructura (NMI = 0,879). Varias comunidades son idénticas en ambos algoritmos y revelan subdivisiones que la clasificación continental no captura: Alaska (131 aeropuertos, articulada por Anchorage) se separa del resto de Norteamérica; el norte de Canadá forma dos comunidades propias, una en torno a Yellowknife y otra en Quebec; Groenlandia e Islandia forman un grupo independiente de Europa; Rusia y Asia Central constituyen una comunidad postsoviética articulada por Moscú-Domodedovo, y el norte de Noruega forma una red regional propia. Ambos algoritmos detectan también pequeñas redes domésticas conectadas al resto a través de un único aeropuerto, como las de Panamá, Israel y las islas Orcadas.

Las diferencias se concentran en el sur del mundo. Leiden recupera a África como una sola comunidad (214 aeropuertos, con Casablanca como hub principal), mientras que Louvain la integra con Oriente Medio y el sur de Asia en una comunidad articulada por Estambul, en la que los aeropuertos asiáticos son apenas la mitad. En Oceanía, Leiden separa la Polinesia Francesa y agrupa a Australia con Nueva Zelanda, mientras que Louvain une Australia con China e Indonesia. En Sudamérica, Leiden forma una única comunidad, mientras que Louvain distingue un bloque Brasil–Argentina de otro andino (Bolivia, Chile y Perú). En conjunto, las comunidades reproducen la geografía a una escala más fina que los continentes y muestran que las regiones periféricas (Alaska, el Ártico canadiense, Groenlandia, la Polinesia) se organizan como subredes que dependen de uno o pocos aeropuertos de entrada, en coherencia con los puntos de articulación identificados en la sección de cohesión.

8.2 Particiones ponderadas por número de aerolíneas

Para incorporar la intensidad de servicio, ambos algoritmos se repiten usando como peso de cada ruta el número de aerolíneas que la operan. Las cuatro particiones (binaria y ponderada para cada algoritmo) se comparan con la partición por continentes mediante:

  • Número de comunidades y tamaño de la mayor.
  • Modularidad binaria y ponderada. Cada partición se evalúa con ambas versiones de \(Q\), para saber si la partición ponderada describe mejor la red con intensidad que la binaria.
  • Información mutua normalizada (NMI). Mide el acuerdo entre dos particiones completas: vale 1 cuando son idénticas y 0 cuando son independientes. Se reporta frente a los continentes y frente a la versión binaria del mismo algoritmo.
  • Porcentaje de rutas intercontinentales internas. Proporción de las rutas que unen aeropuertos de continentes distintos y que quedan dentro de una misma comunidad; indica cuánto de la conexión entre continentes queda absorbido por la estructura comunitaria.
# ---- Datos comunes -----------------------------------------------------------
el_c    <- as_edgelist(g_com, names = FALSE)
inter_e <- cont[el_c[, 1]] != cont[el_c[, 2]] &                  # ruta intercontinental
           cont[el_c[, 1]] != "Otros" & cont[el_c[, 2]] != "Otros"
eb_c    <- if (exists("df_eb")) df_eb$eb else edge_betweenness(g_com, directed = FALSE)
iata_c  <- V(g_com)$iata
w_aer   <- E(giant_component)$num_aerolineas                     # intensidad de servicio

# ---- Particiones ponderadas (cada algoritmo con su propia semilla) ------------
set.seed(123); m_louv_w <- membership(cluster_louvain(g_com, weights = w_aer))
set.seed(123); m_leid_w <- membership(cluster_leiden(g_com, objective_function = "modularity",
                                                     weights = w_aer, n_iterations = -1))

parts <- list(Continentes = part$Continentes,
              `Louvain binaria` = m_louv, `Louvain ponderada` = m_louv_w,
              `Leiden binaria`  = m_leid, `Leiden ponderada`  = m_leid_w)

# ---- Tabla comparativa ---------------------------------------------------------
kable(
  data.frame(
    Partición    = names(parts),
    Comunidades  = sapply(parts, function(m) length(unique(m))),
    `Mayor comunidad` = sapply(parts, function(m) max(table(m))),
    `Q binaria`   = sapply(parts, function(m) modularity(g_com, m)),
    `Q ponderada` = sapply(parts, function(m) modularity(g_com, m, weights = w_aer)),
    `NMI con continentes` = sapply(parts, compare, comm2 = parts$Continentes, method = "nmi"),
    `NMI con su versión binaria` = c(NA, 1, compare(m_louv, m_louv_w, method = "nmi"),
                                     1, compare(m_leid, m_leid_w, method = "nmi")),
    `% rutas intercont. internas` =
      sapply(parts, function(m) 100 * mean(m[el_c[inter_e, 1]] == m[el_c[inter_e, 2]])),
    check.names = FALSE, row.names = NULL
  ),
  digits = c(0, 0, 0, 3, 3, 3, 3, 1), format.args = list(decimal.mark = ",", big.mark = "."),
  caption = "Comparación de particiones: continentes, red binaria y red ponderada por número de aerolíneas."
)
Comparación de particiones: continentes, red binaria y red ponderada por número de aerolíneas.
Partición Comunidades Mayor comunidad Q binaria Q ponderada NMI con continentes NMI con su versión binaria % rutas intercont. internas
Continentes 6 1.228 0,574 0,570 1,000 — 0,0
Louvain binaria 17 751 0,653 0,643 0,641 1,000 38,0
Louvain ponderada 17 638 0,650 0,654 0,609 0,849 44,7
Leiden binaria 17 638 0,660 0,654 0,694 1,000 33,1
Leiden ponderada 16 632 0,655 0,656 0,643 0,920 42,4

8.3 Comunidades transcontinentales

Una comunidad se considera transcontinental cuando al menos dos continentes aportan cinco o más aeropuertos cada uno; este umbral evita clasificar como transcontinentales comunidades con un único aeropuerto de otro continente. Para cada una se identifican los elementos que sostienen su integración: los aeropuertos con más rutas intercontinentales dentro de la comunidad y, entre esas rutas, las de mayor intermediación. El análisis se presenta para Leiden, que obtiene la mayor modularidad, en sus versiones binaria y ponderada.

transcont <- function(m, nombre) {
  m <- as.integer(m)
  do.call(rbind, lapply(sort(unique(m)), function(cc) {
    v   <- which(m == cc)
    tab <- sort(table(cont[v][cont[v] != "Otros"]), decreasing = TRUE)
    if (length(tab) < 2 || tab[2] < 5) return(NULL)                   # criterio transcontinental
    e     <- which(m[el_c[, 1]] == cc & m[el_c[, 2]] == cc & inter_e)  # rutas intercont. internas
    n_v   <- sort(table(c(el_c[e, 1], el_c[e, 2])), decreasing = TRUE)[1:min(4, length(e))]
    top_e <- e[order(-eb_c[e])][1:min(3, length(e))]
    data.frame(
      Red = nombre, Aeropuertos = length(v),
      Composición = paste0(names(tab)[1:min(3, length(tab))], " (", tab[1:min(3, length(tab))], ")",
                           collapse = ", "),
      `Rutas intercont. internas` = length(e),
      `Aeropuertos que la integran` = paste0(iata_c[as.integer(names(n_v))], " (", n_v, ")",
                                             collapse = ", "),
      `Rutas con mayor intermediación` = paste(iata_c[el_c[top_e, 1]], iata_c[el_c[top_e, 2]],
                                               sep = "–", collapse = ", "),
      check.names = FALSE
    )
  }))
}

kable(rbind(transcont(m_leid, "Binaria"), transcont(m_leid_w, "Ponderada")), row.names = FALSE,
      caption = "Comunidades transcontinentales de Leiden (al menos dos continentes con 5 o más aeropuertos). Entre paréntesis, número de rutas intercontinentales internas de cada aeropuerto.")
Comunidades transcontinentales de Leiden (al menos dos continentes con 5 o más aeropuertos). Entre paréntesis, número de rutas intercontinentales internas de cada aeropuerto.
Red Aeropuertos Composición Rutas intercont. internas Aeropuertos que la integran Rutas con mayor intermediación
Binaria 505 Europe (446), Africa (49), Asia (10) 538 AYT (50), RAK (47), TLV (44), LCA (37) ALG–FRA, ALG–LHR, ALG–CDG
Binaria 281 Asia (255), Africa (25), Europe (1) 85 CAI (27), JED (13), HBE (13), KRT (11) KRT–DXB, ASW–JED, CND–ISL
Binaria 161 Europe (116), Asia (45) 228 DME (35), LED (23), TAS (22), EVN (22) LWN–DME, DME–GBB, DME–KOV
Ponderada 365 Oceania (261), Asia (100) 33 SIN (10), KUL (8), DPS (8), PER (6) POM–SIN, SIN–BNE, SIN–DRW
Ponderada 562 Europe (454), Africa (58), Asia (50) 755 ISL (87), AYT (51), TLV (48), SAW (48) OSL–ISL, ALG–FRA, ALG–CDG
Ponderada 241 Asia (217), Africa (24) 74 CAI (26), HBE (12), JED (11), KRT (10) KRT–DXB, ASW–JED, TIP–DXB
Ponderada 153 Europe (110), Asia (43) 206 DME (34), LED (23), TAS (21), DYU (20) LWN–DME, DME–GBB, DME–KOV

En la red binaria, Leiden identifica tres comunidades transcontinentales. La mayor es un bloque euromediterráneo de 505 aeropuertos (446 europeos, 49 africanos y 10 asiáticos), cuya integración descansa en aeropuertos turísticos y de la periferia mediterránea, como Antalya (AYT), Marrakech (RAK), Tel Aviv (TLV) y Larnaca (LCA), y en las rutas que unen Argel con los grandes hubs europeos (ALG–FRA, ALG–LHR, ALG–CDG), que son las de mayor intermediación dentro de la comunidad. La segunda reúne Oriente Medio y el noreste de África (255 aeropuertos asiáticos y 25 africanos), articulada por El Cairo (CAI), Yeda (JED), Alejandría (HBE) y Jartum (KRT), con rutas como Jartum–Dubái y Asuán–Yeda. La tercera es la comunidad postsoviética (116 aeropuertos europeos y 45 asiáticos), sostenida por Moscú-Domodedovo (DME), San Petersburgo (LED), Taskent (TAS) y Ereván (EVN), y por rutas radiales desde Moscú hacia el Cáucaso y Kazajistán.

Debe tenerse en cuenta que parte de este carácter transcontinental se debe a la clasificación de countrycode, que asigna Turquía, Israel y Chipre a Asia y Rusia a Europa; en términos geográficos, las tres comunidades corresponden a regiones de frontera entre continentes más que a conexiones de larga distancia. Las rutas intercontinentales largas, como las transatlánticas o las que unen Europa con Asia oriental, no quedan dentro de ninguna comunidad: conectan comunidades distintas, lo que confirma que la integración global descansa en los conectores intercontinentales identificados en la caracterización de los aeropuertos y no en la estructura comunitaria.

8.4 Comunidades que cambian al ponderar

El NMI resume en un solo número cuánto se parecen dos particiones completas, pero no indica qué comunidades cambian. Además, los números de comunidad son etiquetas arbitrarias: la comunidad 3 de la partición ponderada no corresponde necesariamente a la comunidad 3 de la binaria. Para emparejarlas se usa el coeficiente de Jaccard, que mide el solapamiento entre dos conjuntos de aeropuertos \(A\) (comunidad ponderada) y \(B\) (comunidad binaria):

\[ J(A, B) = \frac{|A \cap B|}{|A \cup B|}, \]

es decir, el número de aeropuertos que comparten dividido entre el total de aeropuertos que aparecen en alguna de las dos. Vale 1 cuando ambas comunidades contienen exactamente los mismos aeropuertos y 0 cuando no comparten ninguno; a diferencia de medir solo qué fracción de \(A\) proviene de \(B\), Jaccard disminuye tanto si la comunidad gana aeropuertos como si los pierde. Para cada comunidad ponderada se calcula \(J\) frente a todas las comunidades binarias y se toma la de mayor valor como su equivalente; se muestran solo las comunidades con \(J < 0{,}9\), es decir, las que cambian de forma apreciable.

# Para cada comunidad ponderada se busca la comunidad binaria con mayor Jaccard
cambios <- do.call(rbind, lapply(sort(unique(m_leid_w)), function(cc) {
  v_w <- which(m_leid_w == cc)
  jac <- sapply(unique(m_leid), function(b) {
    v_b <- which(m_leid == b)
    length(intersect(v_w, v_b)) / length(union(v_w, v_b))
  })
  b      <- unique(m_leid)[which.max(jac)]
  paises <- sort(table(V(g_com)$country[v_w]), decreasing = TRUE)
  data.frame(
    `Aeropuertos (ponderada)` = length(v_w),
    `Aeropuertos (binaria equivalente)` = sum(m_leid == b),
    Jaccard = max(jac),
    `Países principales` = paste0(names(paises)[1:min(3, length(paises))],
                                  " (", paises[1:min(3, length(paises))], ")", collapse = ", "),
    `Hub principal` = iata_c[v_w][which.max(degree(g_com)[v_w])],
    check.names = FALSE
  )
}))

kable(cambios %>% filter(Jaccard < 0.9) %>% arrange(Jaccard),
      digits = c(0, 0, 2, 0, 0), row.names = FALSE,
      format.args = list(decimal.mark = ",", big.mark = "."),
      caption = "Comunidades de Leiden que cambian al ponderar por número de aerolíneas (Jaccard < 0,9 con su comunidad binaria más parecida).")
Comunidades de Leiden que cambian al ponderar por número de aerolíneas (Jaccard < 0,9 con su comunidad binaria más parecida).
Aeropuertos (ponderada) Aeropuertos (binaria equivalente) Jaccard Países principales Hub principal
365 238 0,65 Australia (113), Indonesia (63), Malaysia (32) SIN
388 488 0,80 China (172), Japan (62), Philippines (37) PEK
241 281 0,84 India (68), Iran (38), Saudi Arabia (25) DXB
562 505 0,89 France (54), United Kingdom (44), Turkey (42) AMS

Ponderar por el número de aerolíneas modifica la partición solo de forma parcial. La modularidad apenas cambia: las particiones binarias obtienen prácticamente la misma modularidad ponderada que las ponderadas (0,654 frente a 0,656 en Leiden), y la coincidencia entre ambas versiones es alta (NMI de 0,920 en Leiden y 0,849 en Louvain, que resulta más sensible a la ponderación). La mayoría de las comunidades, incluidas Alaska, el norte de Canadá, Groenlandia, el norte de Noruega y la región postsoviética, no cambian.

Los cambios se concentran en cuatro comunidades y van en una misma dirección: al dar más peso a las rutas operadas por varias aerolíneas, las comunidades se reorganizan en torno a los grandes corredores comerciales, aunque estos crucen fronteras continentales. El cambio más importante ocurre en el sudeste asiático (\(J = 0{,}65\)): Indonesia, Malasia y Singapur dejan la comunidad de China y se integran con Australia y la Polinesia en una comunidad transcontinental de 365 aeropuertos articulada por Singapur (SIN), Kuala Lumpur (KUL), Bali (DPS) y Perth (PER), con rutas como Singapur–Brisbane y Singapur–Darwin. Turquía pasa de la comunidad de Oriente Medio a la europea, en la que Estambul (ISL) se convierte en el principal aeropuerto integrador (87 rutas intercontinentales internas), mientras que la comunidad de Oriente Medio y el sur de Asia pasa a estar articulada por Dubái (DXB). Por último, la Polinesia Francesa, que en la red binaria formaba una comunidad propia, se integra con Oceanía, ya que sus rutas, operadas por pocas aerolíneas, pierden peso frente a las del corredor Australia–Asia.

Como consecuencia, las particiones ponderadas se alejan de la clasificación continental (el NMI con los continentes baja de 0,694 a 0,643 en Leiden) y absorben una mayor proporción de las rutas intercontinentales (del 33,1 % al 42,4 %). La interpretación cambia, por tanto, según lo que se mida: la red binaria describe la estructura de acceso, en la que cada ruta cuenta igual y las comunidades siguen de cerca la geografía y las redes regionales, mientras que la red ponderada describe la estructura de servicio, en la que dominan los corredores con competencia entre aerolíneas y los hubs que concentran esa oferta, como Singapur, Estambul y Dubái, adquieren un papel integrador que la red binaria no revela.

9 Asortatividad y mezcla

La asortatividad mide si los aeropuertos tienden a conectarse con otros parecidos a ellos. Se estudian tres atributos: el grado (¿los aeropuertos grandes se conectan con aeropuertos grandes?), el país y el continente (¿las rutas tienden a quedarse dentro del mismo país o continente?). Todo el análisis se realiza sobre la componente gigante no dirigida y binaria.

Para el grado se usa el coeficiente de asortatividad de Newman, que es la correlación de Pearson entre los grados de los dos extremos de cada arista:

\[ r_k = \frac{\sum_{(u,v)\in E}(x_u - \mu_x)(k_v - \mu_x)}{\sum_{(u,v)\in E}(x_u - \mu_x)^2}, \]

donde \(\mu_x= \frac{1}{2m}\sum_{u,v}k_u\) es el grado medio de los extremos de las aristas. Para atributos categóricos (país y continente) se usa la versión nominal,

\[ r_c = \frac{\sum_c^M e_{cc} - \sum_c^M a_c b_c}{1 - \sum_c a_cb_c}, \]

donde \(e_{cd}\) es la proporción de extremos de aristas que unen la categoría \(c\) con la \(d\) (la matriz de mezcla (MxM) ) y \(a_c = \sum_d e_{cd}\) es la proporción de extremos en la categoría \(c\). El término \(\sum_c a_cb_c\) es la fracción de aristas internas que se esperaría si las conexiones no dependieran de la categorías. Ambos coeficientes valen 1 cuando todas las aristas unen nodos iguales, 0 cuando no hay preferencia y son negativos cuando predominan las conexiones entre nodos distintos.

9.1 Asortatividad observada y modelos nulos

Un coeficiente distinto de cero no basta para afirmar que existe segregación, porque parte de él puede deberse a la propia estructura de la red. Por eso se compara con dos modelos nulos que responden preguntas diferentes:

  • Permutación de etiquetas. Se mantiene fija la red (las mismas rutas) y se reasignan al azar las etiquetas de país o continente entre los aeropuertos, conservando cuántos aeropuertos hay en cada categoría. Responde a la pregunta: dada esta red, ¿la geografía de los aeropuertos está relacionada con quién se conecta con quién? Destruye cualquier relación entre la ubicación y la posición en la red, pero conserva toda la topología (grados, triángulos, comunidades).
  • Rewiring con grados preservados. Se mantienen fijas las etiquetas de cada aeropuerto y su grado, pero se redistribuyen las rutas al azar mediante intercambios de aristas. Responde a la pregunta: dado el número de rutas de cada aeropuerto, ¿se conectan más dentro de su continente de lo que harían al azar? Conserva los grados y la geografía, pero destruye toda la estructura de a quién se conecta cada aeropuerto.

La diferencia es conceptual: la permutación pregunta si las etiquetas están ordenadas sobre una red dada, mientras que el rewiring pregunta si las conexiones están ordenadas dadas las etiquetas y los grados. La permutación no se aplica al grado, porque el grado es una propiedad de la propia red y no una etiqueta externa. En el rewiring, en cambio, el grado sí tiene una referencia nula útil: en una red simple con hubs muy grandes, la restricción de no repetir aristas obliga a los hubs a conectarse con aeropuertos pequeños, lo que produce una asortatividad negativa incluso al azar. Por eso el coeficiente de grado debe compararse con su valor nulo y no con cero.

# ---- Red y atributos -----------------------------------------------------------
g_as    <- delete_edge_attr(giant_component, "weight")
k_as    <- degree(g_as)
cont_as <- V(g_as)$continent; cont_as[is.na(cont_as)] <- "Otros"
f_cont  <- as.integer(factor(cont_as))
f_pais  <- as.integer(factor(V(g_as)$country))

asort <- function(g, fp, fc) c(
  Grado      = assortativity_degree(g, directed = FALSE),
  País       = assortativity_nominal(g, fp, directed = FALSE),
  Continente = assortativity_nominal(g, fc, directed = FALSE)
)
r_obs <- asort(g_as, f_pais, f_cont)

# ---- Umbral de hubs y grilla de grados para el rich club -------------------------
umbral_hub_as <- as.numeric(quantile(k_as, 0.75) + 1.5 * IQR(k_as))   # criterio de Tukey
hubs_as       <- which(k_as >= umbral_hub_as)
k_grid        <- unique(round(exp(seq(log(2), log(max(k_as) - 5), length.out = 30))))

# Coeficiente de rich club: densidad del subgrafo de aeropuertos con grado > k
phi_k <- function(g) sapply(k_grid, function(x) edge_density(induced_subgraph(g, which(k_as > x))))
phi_obs  <- phi_k(g_as)
hub_obs  <- ecount(induced_subgraph(g_as, hubs_as))

# ---- Modelo nulo 1: permutación de etiquetas (red fija) ---------------------------
B_perm <- 1000
set.seed(123)
null_perm <- t(replicate(B_perm, c(
  País       = assortativity_nominal(g_as, sample(f_pais), directed = FALSE),
  Continente = assortativity_nominal(g_as, sample(f_cont), directed = FALSE)
)))

# ---- Modelo nulo 2: rewiring con grados preservados (se guarda en disco) ----------
B_rw        <- 200
archivo_mix <- file.path("cache_resultados", paste0("nulos_asortatividad_B", B_rw, ".rds"))

if (file.exists(archivo_mix)) {
  null_rw <- readRDS(archivo_mix)
} else {
  set.seed(123)
  null_rw <- lapply(seq_len(B_rw), function(b) {
    g_b <- sample_degseq(k_as, method = "edge.switching.simple")
    list(r   = asort(g_b, f_pais, f_cont),
         phi = phi_k(g_b),
         hub = ecount(induced_subgraph(g_b, hubs_as)))
  })
  saveRDS(null_rw, archivo_mix)
}
null_rw_r   <- do.call(rbind, lapply(null_rw, `[[`, "r"))
null_rw_phi <- do.call(rbind, lapply(null_rw, `[[`, "phi"))
null_rw_hub <- sapply(null_rw, `[[`, "hub")

# ---- Tabla de resultados -----------------------------------------------------------
res_perm <- resumen_nulo(r_obs[c("País", "Continente")], null_perm)
res_rw   <- resumen_nulo(r_obs, null_rw_r)

tab_as <- data.frame(
  Atributo  = names(r_obs),
  Observado = fmt_num(r_obs, 3),
  `Permutación: media nula` = c("—", fmt_num(res_perm$media_nula, 3)),
  `Permutación: Z` = c("—", fmt_num(res_perm$z, 1)),
  `Permutación: p` = c("—", fmt_pmc(res_perm$p_mc, B_perm)),
  `Rewiring: media nula` = fmt_num(res_rw$media_nula, 3),
  `Rewiring: Z` = fmt_num(res_rw$z, 1),
  `Rewiring: p` = fmt_pmc(res_rw$p_mc, B_rw),
  check.names = FALSE, row.names = NULL
)
kable(tab_as, align = c("l", rep("r", 7)),
      caption = paste0("Asortatividad observada frente a ", B_perm, " permutaciones de etiquetas y ",
                       B_rw, " redes con grados preservados."))
Asortatividad observada frente a 1000 permutaciones de etiquetas y 200 redes con grados preservados.
Atributo Observado Permutación: media nula Permutación: Z Permutación: p Rewiring: media nula Rewiring: Z Rewiring: p
Grado -0,018 — — — -0,100 16,0 ≤ 0,0100
País 0,444 -0,002 324,4 ≤ 0,0020 -0,002 299,4 ≤ 0,0100
Continente 0,800 -0,002 204,3 ≤ 0,0020 -0,003 186,3 ≤ 0,0100

La asortatividad por grado es prácticamente nula (\(r_k = -0,018\)), lo que a primera vista sugeriría que los aeropuertos se conectan sin preferencia por el tamaño de sus vecinos. Sin embargo, las redes rewireadas con los mismos grados producen un valor claramente negativo (-0,100), porque en una red simple con hubs muy grandes la imposibilidad de repetir aristas obliga a los hubs a conectarse con aeropuertos pequeños. Frente a esa referencia, la red real es mucho menos disortativa de lo que dictan sus grados (\(Z = 16,0\)): los aeropuertos grandes se conectan entre sí con más frecuencia de lo que permitiría el azar, un primer indicio de la estructura de rich club que se analiza más adelante.

La asortatividad por continente (0,800) y por país (0,444) es muy superior a la de ambos modelos nulos, cuyos valores medios son prácticamente cero (\(Z > 180\) en todos los casos). Que la asortatividad continental sea mayor que la nacional indica que la segregación más fuerte ocurre a escala continental: muchas rutas cruzan fronteras nacionales, sobre todo en Europa, pero muy pocas salen del continente. Los dos modelos nulos llevan a la misma conclusión por caminos distintos: la permutación muestra que, sobre esta red, la geografía de los aeropuertos no es aleatoria, y el rewiring muestra que, dado el número de rutas de cada aeropuerto, estas se concentran dentro del mismo continente mucho más de lo que explica el azar. La segregación geográfica no es, por tanto, un efecto de la distribución del grado ni del tamaño de cada continente.

9.2 Matriz de mezcla continental

La matriz de mezcla \(e_{cd}\) muestra qué proporción de los extremos de las rutas une cada par de continentes. Para separar la preferencia real del tamaño de cada continente, se divide entre el valor esperado si las rutas no dependieran del continente, \(a_c a_d\), que es la expectativa del modelo de configuración (rewiring con grados preservados). Un cociente mayor que 1 indica que el par de continentes está más conectado de lo esperado.

el_as <- as_edgelist(g_as, names = FALSE)
ok    <- cont_as[el_as[, 1]] != "Otros" & cont_as[el_as[, 2]] != "Otros"
lev   <- sort(setdiff(unique(cont_as), "Otros"))

# Cada arista se cuenta en ambos sentidos para obtener una matriz simétrica
c1 <- factor(c(cont_as[el_as[ok, 1]], cont_as[el_as[ok, 2]]), levels = lev)
c2 <- factor(c(cont_as[el_as[ok, 2]], cont_as[el_as[ok, 1]]), levels = lev)
e_mix  <- unclass(table(c1, c2)) / (2 * sum(ok))     # matriz de mezcla e_cd
a_mix  <- rowSums(e_mix)
ratio  <- e_mix / outer(a_mix, a_mix)                  # observado / esperado

heat_mix <- function(M, titulo, etiqueta, log_scale = FALSE) {
  d <- as.data.frame(as.table(M), stringsAsFactors = FALSE)
  names(d) <- c("x", "y", "v")
  d$x <- factor(d$x, levels = lev); d$y <- factor(d$y, levels = rev(lev))
  d$f <- if (log_scale) pmax(d$v, 1e-3) else d$v       # evita log(0) si un par no tiene rutas
  p <- ggplot(d, aes(x, y, fill = f)) +
    geom_tile(color = "white") +
    geom_text(aes(label = etiqueta(v)), size = 3.2) +
    labs(title = titulo, x = NULL, y = NULL) +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 30, hjust = 1), panel.grid = element_blank(),
          plot.title = element_text(face = "bold", size = 11))
  if (log_scale) {
    p + scale_fill_gradient2(low = "#e41a1c", mid = "white", high = "#2b5c8f",
                             midpoint = 0, trans = "log10", name = "Obs./esp.",
                             breaks = c(0.01, 0.1, 1, 10), labels = c("0,01", "0,1", "1", "10"))
  } else {
    p + scale_fill_gradient(low = "white", high = "#2b5c8f", name = "Proporción")
    
  }
}

heat_mix(e_mix, "A. Matriz de mezcla (% de extremos de rutas)",
         function(v) paste0(fmt_num(100 * v, 1), " %")) |
heat_mix(ratio, "B. Observado / esperado (a_c · a_d)",
         function(v) fmt_num(v, 2), log_scale = TRUE)

diag_obs <- sum(diag(e_mix))
diag_esp <- sum(a_mix^2)

El 85,7 % de los extremos de las rutas une aeropuertos del mismo continente, frente al 28,2 % que se esperaría si las conexiones no dependieran del continente. Todos los continentes están más conectados consigo mismos de lo esperado (panel B), aunque la magnitud del exceso depende de su tamaño: Oceanía (26,5 veces lo esperado) y África (11,0) concentran pocas rutas del total, por lo que incluso una red interna modesta supera ampliamente su expectativa, mientras que en América, Asia y Europa, que reúnen la mayor parte de las rutas, el exceso es de entre 2,4 y 3,2 veces. Fuera de la diagonal, todos los pares de continentes están menos conectados de lo esperado. Los vínculos relativamente más fuertes son África–Europa (0,69) y África–Asia (0,38), coherentes con la dependencia de África respecto de los hubs europeos y de Oriente Medio. Los más débiles son los que implican cruzar océanos sin un continente intermedio: América–África y América–Asia apenas alcanzan el 4–5 % de lo esperado, y entre Europa y Oceanía no existe ninguna ruta directa en los datos, de modo que todo el tráfico entre ambos debe hacer escala en Asia.

9.3 Conexiones entre aeropuertos de alto grado y rich club

La asortatividad por grado resume la correlación entre grados sobre todas las aristas, y está dominada por la gran masa de rutas que conectan aeropuertos pequeños con hubs. Una red puede ser disortativa en conjunto y, al mismo tiempo, tener a sus hubs densamente conectados entre sí. Esta segunda propiedad se mide con el coeficiente de rich club,

\[ \phi(k) = \frac{2 E_{>k}}{N_{>k}(N_{>k} - 1)}, \]

que es la densidad del subgrafo formado por los \(N_{>k}\) aeropuertos con grado mayor que \(k\), donde \(E_{>k}\) es el número de rutas entre ellos. Como \(\phi(k)\) crece con \(k\) incluso al azar (los hubs tienen más rutas y por tanto más probabilidad de conectarse entre sí), se normaliza con las redes rewireadas:

\[ \rho(k) = \frac{\phi_{obs}(k)}{\bar\phi_{nula}(k)}. \]

Valores de \(\rho(k) > 1\) indican que los aeropuertos de grado mayor que \(k\) están más interconectados de lo que explican sus grados: existe un rich club. Además, se compara directamente el número de rutas entre los hubs definidos por el criterio de Tukey (\(k \ge 19,5\)) con su distribución nula.

rho_k  <- phi_obs / colMeans(null_rw_phi)
rho_lo <- phi_obs / apply(null_rw_phi, 2, quantile, 0.975)
rho_hi <- phi_obs / apply(null_rw_phi, 2, quantile, 0.025)
n_club <- sapply(k_grid, function(x) sum(k_as > x))

ggplot(data.frame(k = k_grid, rho = rho_k, lo = rho_lo, hi = rho_hi), aes(k, rho)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), fill = "grey80", alpha = 0.7) +
  geom_line(color = "#2b5c8f", linewidth = 1) +
  geom_point(color = "#2b5c8f", size = 1.8) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  geom_vline(xintercept = umbral_hub_as, linetype = "dotted", color = "#e41a1c") +
  scale_x_log10() +
  labs(title = "Coeficiente de rich club normalizado",
       subtitle = "Banda: intervalo del 95 % de las redes rewireadas | Línea roja: umbral de hubs (Tukey)",
       x = "Grado k (escala log)", y = expression(rho(k) == phi[obs](k) / bar(phi)[nula](k))) +
  theme_minimal() + theme(panel.grid.minor = element_blank(),
                          plot.title = element_text(face = "bold"))

res_hub <- resumen_nulo(c(hub = hub_obs), matrix(null_rw_hub, ncol = 1))
k_max_rho <- k_grid[which.max(rho_k)]

Entre los 455 hubs hay 9.483 rutas, frente a 9.087 en promedio en las redes con grados preservados (1,04 veces lo esperado; \(Z = 14,2\), p ≤ 0,0100). El coeficiente normalizado \(\rho(k)\) se mantiene por encima de 1 y crece con el grado: es prácticamente 1 para valores bajos de \(k\), cuando el club incluye a casi toda la red, y aumenta hasta cerca de 1,15–1,20 entre los aeropuertos con más de 80 rutas, por encima de la banda del 95 % de las redes rewireadas. Existe, por tanto, una estructura de rich club moderada: los aeropuertos de mayor grado tienen entre un 5 % y un 20 % más de rutas entre sí de lo que explica su número de conexiones. Para los grados más altos (\(k > 150\)), el club se reduce a muy pocos aeropuertos y la variabilidad nula es tan grande que no permite concluir.

Este resultado no es equivalente a la asortatividad por grado; ambas miden cosas distintas. La asortatividad se calcula sobre todas las rutas y está dominada por las miles de conexiones entre hubs y aeropuertos pequeños, mientras que el rich club se fija solo en si los aeropuertos de mayor grado mantienen rutas entre sí. En esta red ambos análisis son complementarios: la asortatividad observada es casi nula, pero mucho mayor que la que impone la distribución del grado, y el rich club muestra dónde se origina esa diferencia, en el núcleo de hubs que se conectan entre sí más de lo esperado. La red combina así dos rasgos: cada hub distribuye conexiones hacia una periferia de aeropuertos pequeños y, al mismo tiempo, los hubs forman un núcleo interconectado, en coherencia con el núcleo de coreness máxima y los clanes de hubs identificados en la sección de cohesión. Este núcleo es el que sostiene la conectividad de mundo pequeño y explica la robustez observada entre los hubs.

10 Vulnerabilidad y robustez

Esta sección evalúa cómo se degrada la conectividad de la red cuando se eliminan aeropuertos. Se simulan ataques dirigidos, que eliminan primero los aeropuertos más importantes según una medida de centralidad, y se comparan con fallas aleatorias, que eliminan aeropuertos al azar. Todo el análisis se realiza sobre la componente gigante no dirigida y binaria.

10.1 Diseño de las simulaciones

Después de cada bloque de eliminaciones se registran cuatro medidas, todas relativas a la red original de \(n_0\) aeropuertos:

  • Tamaño relativo de la componente gigante, \(S = N_{GC}/n_0\): proporción de los aeropuertos originales que siguen en la mayor componente conexa.
  • Eficiencia global relativa, \(E/E_0\), donde \[ E = \frac{1}{n_0(n_0-1)}\sum_{i \neq j}\frac{1}{d_G(i,j)}, \] con \(d_G\) medida en aristas (vuelos) y \(1/d_G = 0\) si \(i\) y \(j\) no están conectados o si alguno fue eliminado. A diferencia de la distancia promedio, la eficiencia está bien definida cuando la red se fragmenta, y captura a la vez la pérdida de conexión y el alargamiento de los caminos.
  • Número de componentes: indica si la red se rompe en muchos fragmentos pequeños.
  • Alcance intercontinental: proporción de los pares de aeropuertos de continentes distintos que siguen conectados por algún camino. Permite distinguir la pérdida de conexión local de la pérdida de integración global.

Las estrategias de eliminación son:

Estrategia Orden de eliminación
Aleatoria Al azar; se repite 20 veces y se reporta la media y el intervalo del 95 %
Grado (estático) Según el grado en la red original
Grado (adaptativo) Según el grado, recalculado sobre la red que queda
Intermediación (estático) Según la intermediación en la red original
Intermediación (adaptativo) Según la intermediación, recalculada sobre la red que queda
Participación (estático) Según el coeficiente de participación continental \(P\), entre aeropuertos con al menos 5 rutas

En los ataques estáticos, el orden se fija al inicio; en los adaptativos, la centralidad se recalcula después de cada bloque de 10 eliminaciones (el 0,3 % de la red), de modo que el ataque responde a cómo se reorganizan los caminos a medida que la red se daña. Recalcular después de cada aeropuerto individual tendría un costo computacional muy alto para la intermediación y produce prácticamente las mismas curvas con bloques de este tamaño. Se eliminan hasta 320 aeropuertos (el 10 % de la red). Como resumen se usa el índice de robustez \(R = \frac{1}{Q}\sum_{q} S(q)\), el promedio de \(S\) a lo largo de las eliminaciones: valores menores indican un ataque más dañino.

# ---- Red, referencias y medidas -------------------------------------------------
g_rob    <- delete_edge_attr(giant_component, "weight")
n0       <- vcount(g_rob)
nm_rob   <- V(g_rob)$name
cont_rob <- V(g_rob)$continent; cont_rob[is.na(cont_rob)] <- "Otros"
E0       <- global_efficiency(g_rob)

# Pares de aeropuertos de distinto continente que están en la misma componente
pares_inter <- function(memb, ct) {
  ok  <- ct != "Otros"
  tab <- table(memb[ok], ct[ok])
  sum((rowSums(tab)^2 - rowSums(tab^2)) / 2)
}
inter0 <- pares_inter(rep(1, n0), cont_rob)

metricas <- function(g) {
  nr <- vcount(g)
  if (nr < 2) return(c(S = 0, E = 0, Componentes = nr, Inter = 0))
  cc <- components(g)
  ct <- V(g)$continent; ct[is.na(ct)] <- "Otros"
  c(S           = max(cc$csize) / n0,
    E           = global_efficiency(g) * nr * (nr - 1) / (n0 * (n0 - 1)) / E0,
    Componentes = cc$no,
    Inter       = pares_inter(cc$membership, ct) / inter0)
}

# ---- Ataques ------------------------------------------------------------------------
pasos <- seq(0, 320, by = 10)

ataque_estatico <- function(g, orden, pasos) {
  t(sapply(pasos, function(q) metricas(delete_vertices(g, orden[seq_len(q)]))))
}

ataque_adaptativo <- function(g, centralidad, pasos) {
  res <- matrix(NA_real_, length(pasos), 4, dimnames = list(NULL, c("S", "E", "Componentes", "Inter")))
  res[1, ] <- metricas(g)
  for (i in seq_along(pasos)[-1]) {
    cen    <- centralidad(g)
    quitar <- names(sort(cen, decreasing = TRUE))[seq_len(pasos[i] - pasos[i - 1])]
    g      <- delete_vertices(g, quitar)
    res[i, ] <- metricas(g)
  }
  res
}

# ---- Órdenes estáticos ----------------------------------------------------------------
k_rob     <- degree(g_rob)
b_rob     <- betweenness(g_rob, directed = FALSE)
P_rob     <- ifelse(df_med$grado >= 5 & !is.na(df_med$P), df_med$P, 0)
ord_grado <- nm_rob[order(-k_rob)]
ord_betw  <- nm_rob[order(-b_rob)]
ord_P     <- nm_rob[order(-P_rob, -k_rob)]

# ---- Hubs regionales y conectores intercontinentales (sección de caracterización) ----
es_hub  <- df_med$grado >= umbral_hub
idx_reg <- which(es_hub & !is.na(df_med$P) & df_med$P <  P_cut)
idx_con <- which(es_hub & !is.na(df_med$P) & df_med$P >= P_cut)
ord_reg <- nm_rob[idx_reg[order(-df_med$grado[idx_reg])]]
ord_con <- nm_rob[idx_con[order(-df_med$grado[idx_con])]]
pasos_h <- seq(0, floor(min(length(ord_reg), length(ord_con)) / 10) * 10, by = 10)

# ---- Escenario regional: todos los hubs de una región a la vez ---------------------------
region   <- "Europe"
v_region <- which(es_hub & cont_rob == region)
n_region <- length(v_region)

# ---- Simulaciones (se guardan en disco: la primera ejecución tarda varios minutos) -----
B_rand      <- 20
B_glob      <- 50
archivo_rob <- file.path("cache_resultados", "robustez.rds")

if (file.exists(archivo_rob)) {
  rob <- readRDS(archivo_rob)
} else {
  rob <- list()
  rob$grado_est <- ataque_estatico(g_rob, ord_grado, pasos)
  rob$betw_est  <- ataque_estatico(g_rob, ord_betw,  pasos)
  rob$P_est     <- ataque_estatico(g_rob, ord_P,     pasos)
  rob$grado_ad  <- ataque_adaptativo(g_rob, function(g) degree(g), pasos)
  rob$betw_ad   <- ataque_adaptativo(g_rob, function(g) betweenness(g, directed = FALSE), pasos)
  set.seed(123)
  rob$aleatorio <- lapply(seq_len(B_rand), function(b) ataque_estatico(g_rob, sample(nm_rob), pasos))
  rob$regional  <- ataque_estatico(g_rob, ord_reg, pasos_h)
  rob$conector  <- ataque_estatico(g_rob, ord_con, pasos_h)
  rob$region    <- metricas(delete_vertices(g_rob, v_region))
  set.seed(123)
  rob$global    <- t(replicate(B_glob, metricas(delete_vertices(g_rob, sample(which(es_hub), n_region)))))
  saveRDS(rob, archivo_rob)
}

# ---- Datos en formato largo para graficar ------------------------------------------------
a_largo <- function(M, estrategia, ps) data.frame(q = ps, f = ps / n0, M, Estrategia = estrategia)

alea_media <- Reduce(`+`, rob$aleatorio) / B_rand
alea_lo    <- apply(simplify2array(rob$aleatorio), c(1, 2), quantile, 0.025)
alea_hi    <- apply(simplify2array(rob$aleatorio), c(1, 2), quantile, 0.975)

df_ataques <- rbind(
  a_largo(alea_media,    "Aleatoria",                  pasos),
  a_largo(rob$grado_est, "Grado (estático)",           pasos),
  a_largo(rob$grado_ad,  "Grado (adaptativo)",         pasos),
  a_largo(rob$betw_est,  "Intermediación (estático)",  pasos),
  a_largo(rob$betw_ad,   "Intermediación (adaptativo)", pasos),
  a_largo(rob$P_est,     "Participación (estático)",   pasos)
)

10.2 Ataques dirigidos frente a fallas aleatorias

etiq_met <- c(S = "A. Tamaño relativo de la componente gigante",
              E = "B. Eficiencia global relativa",
              Componentes = "C. Número de componentes",
              Inter = "D. Alcance intercontinental")

largo <- df_ataques %>%
  pivot_longer(c(S, E, Componentes, Inter), names_to = "Métrica", values_to = "valor") %>%
  mutate(Métrica = factor(etiq_met[Métrica], levels = etiq_met))

# Banda del 95 % de las eliminaciones aleatorias
banda <- do.call(rbind, lapply(names(etiq_met), function(m)
  data.frame(f = pasos / n0, Métrica = etiq_met[[m]], lo = alea_lo[, m], hi = alea_hi[, m])))
banda$Métrica <- factor(banda$Métrica, levels = etiq_met)

col_est <- c("Aleatoria" = "grey40", "Grado (estático)" = "#2b5c8f", "Grado (adaptativo)" = "#2b5c8f",
             "Intermediación (estático)" = "#e41a1c", "Intermediación (adaptativo)" = "#e41a1c",
             "Participación (estático)" = "#4daf4a")
lin_est <- c("Aleatoria" = "solid", "Grado (estático)" = "solid", "Grado (adaptativo)" = "dashed",
             "Intermediación (estático)" = "solid", "Intermediación (adaptativo)" = "dashed",
             "Participación (estático)" = "solid")

ggplot(largo, aes(f, valor, color = Estrategia, linetype = Estrategia)) +
  geom_ribbon(data = banda, aes(x = f, ymin = lo, ymax = hi), inherit.aes = FALSE,
              fill = "grey80", alpha = 0.6) +
  geom_line(linewidth = 0.9) +
  facet_wrap(~ Métrica, scales = "free_y") +
  scale_x_continuous(labels = function(x) paste0(fmt_num(100 * x, 0), " %")) +
  scale_color_manual(values = col_est) + scale_linetype_manual(values = lin_est) +
  labs(x = "Proporción de aeropuertos eliminados", y = NULL, color = NULL, linetype = NULL,
       title = "Robustez de la red frente a ataques dirigidos y fallas aleatorias",
       subtitle = "Banda gris: intervalo del 95 % de las eliminaciones aleatorias") +
  theme_minimal() +
  theme(legend.position = "bottom", panel.grid.minor = element_blank(),
        strip.text = element_text(face = "bold"), plot.title = element_text(face = "bold"))

resumen_ataque <- function(M) {
  i5 <- which(pasos == 160); i10 <- length(pasos)
  c(R = mean(M[, "S"]), S5 = M[i5, "S"], S10 = M[i10, "S"],
    E10 = M[i10, "E"], C10 = M[i10, "Componentes"], I10 = M[i10, "Inter"])
}
estrategias <- list(`Aleatoria (media)` = alea_media,
                    `Grado (estático)` = rob$grado_est, `Grado (adaptativo)` = rob$grado_ad,
                    `Intermediación (estático)` = rob$betw_est,
                    `Intermediación (adaptativo)` = rob$betw_ad,
                    `Participación (estático)` = rob$P_est)
tab_rob <- as.data.frame(do.call(rbind, lapply(estrategias, resumen_ataque)))

kable(
  data.frame(Estrategia = names(estrategias),
             `Robustez R` = fmt_num(tab_rob$R, 3),
             `S (5 %)` = fmt_num(tab_rob$S5, 3),
             `S (10 %)` = fmt_num(tab_rob$S10, 3),
             `E/E0 (10 %)` = fmt_num(tab_rob$E10, 3),
             `Componentes (10 %)` = fmt_num(tab_rob$C10, 0),
             `Alcance intercont. (10 %)` = fmt_por(tab_rob$I10, 1),
             check.names = FALSE, row.names = NULL),
  align = c("l", rep("r", 6)),
  caption = "Resumen de las estrategias de eliminación: índice de robustez y estado de la red tras eliminar el 5 % y el 10 % de los aeropuertos."
)
Resumen de las estrategias de eliminación: índice de robustez y estado de la red tras eliminar el 5 % y el 10 % de los aeropuertos.
Estrategia Robustez R S (5 %) S (10 %) E/E0 (10 %) Componentes (10 %) Alcance intercont. (10 %)
Aleatoria (media) 0,932 0,932 0,865 0,749 83 74,9 %
Grado (estático) 0,733 0,772 0,256 0,039 958 6,4 %
Grado (adaptativo) 0,717 0,741 0,282 0,039 1.123 7,5 %
Intermediación (estático) 0,622 0,581 0,361 0,081 1.131 12,5 %
Intermediación (adaptativo) 0,437 0,297 0,072 0,017 1.003 0,2 %
Participación (estático) 0,807 0,862 0,458 0,218 498 10,7 %
peor <- names(estrategias)[which.min(tab_rob$R)]

Las fallas aleatorias apenas afectan a la red: tras eliminar al azar el 10 % de los aeropuertos, la componente gigante conserva el % de los aeropuertos originales, cerca del 90 % máximo posible, y la eficiencia se mantiene en el % de su valor inicial, porque la mayoría de los aeropuertos tiene muy pocas rutas y una eliminación al azar casi siempre afecta a la periferia. Los ataques dirigidos, en cambio, degradan la red con rapidez: con el 10 % de los aeropuertos eliminados, la eficiencia cae por debajo del 10 % de su valor inicial en los ataques por grado e intermediación, y la red se divide en alrededor de mil componentes. Esta asimetría entre la robustez frente a fallas aleatorias y la fragilidad frente a ataques dirigidos es característica de las redes con distribución del grado muy heterogénea.

La estrategia más dañina es la intermediación adaptativa (\(R = 0,437\)). Su curva muestra un colapso abrupto en torno al 2,7 % de eliminaciones: la componente gigante cae del 72 % al 45 % de los aeropuertos y el alcance intercontinental, de cerca del 55 % a menos del 20 %, en apenas unas decenas de eliminaciones. Al recalcular la intermediación después de cada bloque, el ataque identifica los puentes que quedan tras eliminar los anteriores, hasta separar los continentes entre sí. Al 10 %, la componente gigante reúne solo el % de los aeropuertos y prácticamente ningún par de continentes sigue conectado ( %). La intermediación estática es más dañina que el grado en las primeras eliminaciones (componente gigante del % frente al % al 5 %), porque elimina primero los aeropuertos puente, aunque su efecto se agota antes, porque la intermediación original deja de identificar los nuevos puentes. La diferencia entre estático y adaptativo es, en cambio, pequeña para el grado (0,733 frente a 0,717), porque eliminar un hub apenas altera el orden de los grados de los demás. El ataque por grado mantiene la componente gigante durante más tiempo, pero reduce la eficiencia casi tanto como el de intermediación y presenta una caída brusca cerca del 10 %, cuando se eliminan los últimos hubs del núcleo y la componente gigante pasa del 55 % al %.

El ataque por participación continental es el menos dañino de los dirigidos (\(R = 0,807\)) y genera muchas menos componentes, pero produce un tipo de daño distinto: alrededor del 9 % de eliminaciones, cuando se eliminan los principales conectores intercontinentales, el alcance entre continentes se desploma del 38 % al 12 %. Al 10 %, la componente gigante todavía conserva el % de los aeropuertos y el % de la eficiencia (los valores más altos de los ataques dirigidos), pero solo el % de los pares entre continentes sigue conectado. La red no se pulveriza, sino que se separa en grandes bloques regionales que siguen funcionando internamente, pero sin conexión entre sí.

10.3 Hubs regionales frente a conectores intercontinentales

Para analizar si los dos tipos de hubs identificados en la caracterización de los aeropuertos sostienen la red de manera distinta, se eliminan por separado los hubs regionales y los conectores intercontinentales, en ambos casos en orden de grado decreciente.

df_tipo <- rbind(a_largo(rob$regional, "Hubs regionales", pasos_h),
                 a_largo(rob$conector, "Conectores intercontinentales", pasos_h)) %>%
  pivot_longer(c(S, E, Componentes, Inter), names_to = "Métrica", values_to = "valor") %>%
  mutate(Métrica = factor(etiq_met[Métrica], levels = etiq_met))

ggplot(df_tipo, aes(q, valor, color = Estrategia)) +
  geom_line(linewidth = 1) +
  facet_wrap(~ Métrica, scales = "free_y", nrow = 1) +
  scale_color_manual(values = c("Hubs regionales" = "#2b5c8f",
                                "Conectores intercontinentales" = "#e41a1c")) +
  labs(x = "Número de hubs eliminados", y = NULL, color = NULL,
       title = "Eliminación de hubs regionales frente a conectores intercontinentales") +
  theme_minimal() +
  theme(legend.position = "bottom", panel.grid.minor = element_blank(),
        strip.text = element_text(face = "bold", size = 8), plot.title = element_text(face = "bold"))

ult <- nrow(rob$regional)

Tras eliminar 220 hubs de cada tipo, la pérdida de hubs regionales deja la red en 311 componentes, con una componente gigante del 79,4 % y un alcance intercontinental del 65,0 %, mientras que la pérdida de conectores intercontinentales produce 529 componentes, una componente gigante del 66,2 % y un alcance intercontinental del 37,7 %. Los dos tipos de hubs se asocian a mecanismos de fragmentación distintos: los hubs regionales son la puerta de entrada de su periferia, de modo que su pérdida desprende grupos de aeropuertos pequeños y aumenta el número de componentes; los conectores intercontinentales, en cambio, sostienen los caminos cortos entre regiones, de modo que su pérdida alarga los trayectos y reduce la eficiencia y el alcance entre continentes, aunque buena parte de los aeropuertos siga perteneciendo a la componente gigante.

10.4 Pérdida simultánea de una región

Por último, se simula el cierre simultáneo de todos los hubs de una misma región, un escenario comparable a los cierres del espacio aéreo europeo por la erupción del Eyjafjallajökull en 2010, y se compara con la eliminación del mismo número de hubs elegidos al azar entre todos los hubs del mundo (50 repeticiones).

glob_media <- colMeans(rob$global)
glob_lo    <- apply(rob$global, 2, quantile, 0.025)
glob_hi    <- apply(rob$global, 2, quantile, 0.975)

fmt_m <- list(S = function(x) fmt_por(x, 1), E = function(x) fmt_por(x, 1),
              Componentes = function(x) fmt_num(x, 0), Inter = function(x) fmt_por(x, 1))
kable(
  data.frame(
    Medida = c("Tamaño relativo de la componente gigante", "Eficiencia global relativa",
               "Número de componentes", "Alcance intercontinental"),
    `Hubs de la región` = mapply(function(f, x) f(x), fmt_m, rob$region),
    `Hubs distribuidos: media` = mapply(function(f, x) f(x), fmt_m, glob_media),
    `IC 95 %` = paste0("[", mapply(function(f, x) f(x), fmt_m, glob_lo), "; ",
                       mapply(function(f, x) f(x), fmt_m, glob_hi), "]"),
    check.names = FALSE, row.names = NULL
  ),
  align = c("l", "r", "r", "r"),
  caption = paste0("Eliminación de los ", n_region, " hubs de ", region,
                   " frente a la eliminación del mismo número de hubs distribuidos globalmente.")
)
Eliminación de los 165 hubs de Europe frente a la eliminación del mismo número de hubs distribuidos globalmente.
Medida Hubs de la región Hubs distribuidos: media IC 95 %
Tamaño relativo de la componente gigante 85,2 % 84,0 % [79,6 %; 87,8 %]
Eficiencia global relativa 67,9 % 66,1 % [60,9 %; 71,8 %]
Número de componentes 219 249 [186; 323]
Alcance intercontinental 68,0 % 70,6 % [63,8 %; 76,6 %]

Eliminar los 165 hubs de Europe deja la componente gigante en el 85,2 % de los aeropuertos y el alcance intercontinental en el 68,0 %, valores que caen dentro del intervalo del 95 % obtenido al eliminar el mismo número de hubs distribuidos por el mundo (84,0 % y 70,6 % en promedio). La pérdida simultánea de toda una región no resulta, por tanto, más dañina para la conectividad global que una pérdida equivalente repartida. Esto refleja la alta redundancia de la red europea: sus aeropuertos no clasificados como hubs mantienen numerosas rutas entre sí y con aeropuertos del norte de África y Oriente Medio, de modo que la mayoría sigue conectada al resto de la red a pesar del cierre de todos sus hubs. Además, como Europa concentra una parte importante de los hubs del mundo, la eliminación distribuida también incluye, en promedio, una fracción considerable de hubs europeos. El resultado podría ser distinto en regiones con menos redundancia interna, como Oceanía o África, cuya conectividad depende de muy pocos aeropuertos de entrada.

11 Sensibilidad a la definición de la red

Los resultados anteriores se obtuvieron con una red no dirigida y binaria: dos aeropuertos están conectados si existe al menos una ruta entre ellos, en cualquier sentido, y todas las rutas pesan lo mismo. Esta sección evalúa qué conclusiones dependen de esas decisiones, comparando con (i) la red dirigida, (ii) la red ponderada por número de aerolíneas y (iii) una red que exige reciprocidad, es decir, que solo conserva las rutas operadas en ambos sentidos.

11.1 Red dirigida frente a red no dirigida

# Red dirigida binaria y simple sobre los mismos aeropuertos de la componente gigante
g_d <- induced_subgraph(g_directed, which(cl$membership == largest_cluster_id))
g_d <- igraph::simplify(g_d, remove.multiple = TRUE, remove.loops = TRUE, edge.attr.comb = "ignore")

k_in  <- degree(g_d, mode = "in")
k_out <- degree(g_d, mode = "out")
c_deb <- components(g_d, mode = "weak")
c_fue <- components(g_d, mode = "strong")

kable(
  data.frame(
    Medida = c("Arcos dirigidos", "Reciprocidad por arcos", "Reciprocidad por díadas conectadas",
               "Correlación de Spearman entre grado de entrada y de salida",
               "Aeropuertos con grado de entrada igual al de salida",
               "Componentes débiles", "Componentes fuertes",
               "Aeropuertos en la mayor componente fuerte"),
    Valor = c(fmt_num(ecount(g_d), 0),
              fmt_num(reciprocity(g_d, mode = "default"), 3),
              fmt_num(reciprocity(g_d, mode = "ratio"), 3),
              fmt_num(cor(k_in, k_out, method = "spearman"), 3),
              fmt_por(mean(k_in == k_out), 1),
              fmt_num(c_deb$no, 0), fmt_num(c_fue$no, 0),
              paste0(fmt_num(max(c_fue$csize), 0), " (", fmt_por(max(c_fue$csize) / vcount(g_d), 1), ")"))
  ),
  row.names = FALSE, align = c("l", "r"),
  caption = "Propiedades de la red dirigida (componente gigante)."
)
Propiedades de la red dirigida (componente gigante).
Medida Valor
Arcos dirigidos 36.860
Reciprocidad por arcos 0,978
Reciprocidad por díadas conectadas 0,957
Correlación de Spearman entre grado de entrada y de salida 0,987
Aeropuertos con grado de entrada igual al de salida 82,8 %
Componentes débiles 1
Componentes fuertes 42
Aeropuertos en la mayor componente fuerte 3.147 (98,7 %)

La red dirigida es casi simétrica. De los 36.860 arcos, el 97,8 % tiene su ruta de regreso (reciprocidad por arcos, \(2M/(2M + A)\), donde \(M\) es el número de pares mutuos y \(A\) el de pares asimétricos), y el 95,7 % de los pares de aeropuertos conectados lo están en ambos sentidos.

El grado de entrada y el de salida están casi perfectamente correlacionados (Spearman de 0,987), y en el 82,8 % de los aeropuertos son exactamente iguales. La red tiene una sola componente débil, y la mayor componente fuerte, en la que cada aeropuerto puede llegar a cualquier otro y volver respetando el sentido de las rutas, reúne el 98,7 % de los aeropuertos. Las 41 componentes fuertes restantes son aeropuertos aislados que solo pueden alcanzarse o abandonarse en un sentido, típicamente extremos de rutas triangulares.

Para determinar si esto se traduce en las conclusiones sobre centralidad y comunidades, se recalculan las medidas que dependen de la dirección y se comparan con sus equivalentes no dirigidas. Para las comunidades se usa el algoritmo Infomap, que, a diferencia de Louvain y Leiden en igraph, admite redes dirigidas: modela un paseo aleatorio que sigue el sentido de las rutas y busca la partición que describe ese paseo con la menor longitud de código.

# ---- Centralidad: dirigida vs. no dirigida ------------------------------------------
k_top   <- ceiling(0.01 * vcount(g_d))
top_sol <- function(x, y) length(intersect(order(-x)[1:k_top], order(-y)[1:k_top])) / k_top
iata_d  <- V(g_d)$iata
top5    <- function(x) paste(iata_d[order(-x)[1:5]], collapse = ", ")

pr_d <- page_rank(g_d, directed = TRUE)$vector
bt_d <- betweenness(g_d, directed = TRUE)

comparaciones <- list(
  `Grado de entrada vs. grado`   = list(df_med$grado,          k_in),
  `Grado de salida vs. grado`    = list(df_med$grado,          k_out),
  `Intermediación`               = list(df_med$intermediacion, bt_d),
  `PageRank`                     = list(df_med$pagerank,       pr_d)
)

kable(
  data.frame(
    Medida = names(comparaciones),
    `Spearman` = sapply(comparaciones, function(p) fmt_num(cor(p[[1]], p[[2]], method = "spearman"), 3)),
    `Superposición top-k` = sapply(comparaciones, function(p) fmt_por(top_sol(p[[1]], p[[2]]), 0)),
    `Top 5 no dirigida` = sapply(comparaciones, function(p) top5(p[[1]])),
    `Top 5 dirigida` = sapply(comparaciones, function(p) top5(p[[2]])),
    check.names = FALSE, row.names = NULL
  ),
  align = c("l", "r", "r", "l", "l"),
  caption = paste0("Centralidad en la red no dirigida frente a la dirigida (superposición de los ", k_top, " primeros aeropuertos).")
)
Centralidad en la red no dirigida frente a la dirigida (superposición de los 32 primeros aeropuertos).
Medida Spearman Superposición top-k Top 5 no dirigida Top 5 dirigida
Grado de entrada vs. grado 0,987 100 % AMS, FRA, CDG, ISL, ATL FRA, CDG, AMS, ISL, ATL
Grado de salida vs. grado 0,985 100 % AMS, FRA, CDG, ISL, ATL FRA, CDG, AMS, ISL, ATL
Intermediación 0,956 100 % CDG, LAX, ANC, DXB, FRA CDG, LAX, DXB, ANC, FRA
PageRank 0,976 97 % ATL, ORD, ISL, DFW, DEN ATL, ORD, ISL, DEN, DFW
# ---- Comunidades: Infomap dirigido vs. no dirigido -----------------------------------
set.seed(123); m_inf_d <- membership(cluster_infomap(g_d, nb.trials = 5))
set.seed(123); m_inf_u <- membership(cluster_infomap(g_com, nb.trials = 5))

kable(
  data.frame(
    Comparación = c("Infomap dirigido vs. Infomap no dirigido",
                    "Infomap no dirigido vs. Leiden (binaria)",
                    "Infomap dirigido vs. Leiden (binaria)"),
    `Comunidades` = c(paste(length(unique(m_inf_d)), "vs.", length(unique(m_inf_u))),
                      paste(length(unique(m_inf_u)), "vs.", length(unique(m_leid))),
                      paste(length(unique(m_inf_d)), "vs.", length(unique(m_leid)))),
    NMI = fmt_num(c(compare(m_inf_d, m_inf_u, method = "nmi"),
                    compare(m_inf_u, m_leid,  method = "nmi"),
                    compare(m_inf_d, m_leid,  method = "nmi")), 3),
    check.names = FALSE, row.names = NULL
  ),
  align = c("l", "r", "r"),
  caption = "Comparación de particiones al conservar la dirección de las rutas."
)
Comparación de particiones al conservar la dirección de las rutas.
Comparación Comunidades NMI
Infomap dirigido vs. Infomap no dirigido 153 vs. 162 0,977
Infomap no dirigido vs. Leiden (binaria) 162 vs. 17 0,719
Infomap dirigido vs. Leiden (binaria) 153 vs. 17 0,721

Conservar la dirección tampoco modifica las conclusiones sobre centralidad y comunidades. Los grados de entrada y de salida ordenan a los aeropuertos prácticamente igual que el grado no dirigido, y la intermediación y el PageRank dirigidos mantienen correlaciones altas con sus versiones no dirigidas. La partición de Infomap dirigida es casi idéntica a la no dirigida (153 frente a 162 comunidades; NMI = 0,977). Infomap identifica muchas más comunidades que Leiden (162 frente a 17), porque su criterio, la compresión de un paseo aleatorio, separa cada pequeña subred regional en la que el paseo tiende a quedar atrapado, mientras que la modularidad favorece bloques más grandes. Aun así, ambas particiones coinciden en lo esencial (NMI = 0,719), y esa coincidencia es la misma con la versión dirigida (0,721). Las diferencias entre particiones se deben, por tanto, al algoritmo y no a la dirección de las rutas.

11.2 Red ponderada frente a red binaria

En la red ponderada, el peso de cada ruta es el número de aerolíneas que la operan, \(a_{ij}\). Este peso mide intensidad de servicio: más aerolíneas implican más frecuencias, más opciones y menor dependencia de un solo operador. Los algoritmos de caminos mínimos, sin embargo, interpretan el peso como un costo que se suma a lo largo del camino. Usar \(a_{ij}\) directamente invertiría su significado, porque las rutas mejor servidas parecerían “más largas”. Por eso, para los caminos mínimos se usa la transformación

\[ c_{ij} = \frac{1}{a_{ij}}, \]

que es positiva (como exige el algoritmo de Dijkstra) y decreciente en la intensidad: una ruta con una sola aerolínea tiene costo 1 y una con cuatro aerolíneas, costo 0,25. Con esta transformación, un camino de dos tramos bien servidos puede ser “más corto” que un tramo directo operado por una sola aerolínea, lo que representa la idea de que la conexión efectiva depende tanto del número de escalas como de la disponibilidad de servicio. Se compara la intermediación así obtenida con la intermediación binaria (en número de tramos).

costo  <- 1 / E(giant_component)$num_aerolineas
bt_w   <- betweenness(giant_component, directed = FALSE, weights = costo)
iata_u <- V(giant_component)$iata

top_w <- function(x) paste(iata_u[order(-x)[1:8]], collapse = ", ")
kable(
  data.frame(
    Red = c("Binaria (tramos)", "Ponderada (costo 1 / aerolíneas)"),
    `Ocho aeropuertos con mayor intermediación` = c(top_w(df_med$intermediacion), top_w(bt_w)),
    check.names = FALSE, row.names = NULL
  ),
  align = c("l", "l"),
  caption = paste0("Intermediación binaria frente a ponderada. Spearman = ",
                   fmt_num(cor(df_med$intermediacion, bt_w, method = "spearman"), 3),
                   "; superposición del top-", k_top, " = ",
                   fmt_por(top_sol(df_med$intermediacion, bt_w), 0), ".")
)
Intermediación binaria frente a ponderada. Spearman = 0,796; superposición del top-32 = 84 %.
Red Ocho aeropuertos con mayor intermediación
Binaria (tramos) CDG, LAX, ANC, DXB, FRA, AMS, PEK, ORD
Ponderada (costo 1 / aerolíneas) LHR, LAX, ATL, JFK, SIN, MIA, ORD, ANC
sube <- order(-(rank(-df_med$intermediacion) - rank(-bt_w)))[1:5]

La ponderación sí modifica qué aeropuertos actúan como puentes. La intermediación ponderada mantiene una correlación de 0,796 con la binaria, y el 84 % de los 32 aeropuertos más intermediarios coincide, pero el orden en la cima cambia. En la red binaria dominan los grandes hubs europeos y asiáticos (París, Fráncfort, Ámsterdam, Dubái, Pekín) junto con los puentes regionales como Anchorage. Al ponderar por número de aerolíneas, en cambio, pasan al frente Londres-Heathrow, Atlanta, Nueva York-JFK, Singapur y Miami, aeropuertos que concentran rutas operadas por muchas aerolíneas. Los caminos de menor costo se desplazan hacia los corredores con más servicio, como el transatlántico (Londres–Nueva York), el de Miami hacia América Latina y el de Singapur hacia Oceanía, aunque impliquen algún tramo adicional. Este resultado es coherente con la detección de comunidades, donde la ponderación reorganizó la partición en torno a los mismos corredores (NMI entre las particiones binaria y ponderada de Leiden = 0,920), y con ella Singapur pasó a articular una comunidad que une Asia y Oceanía. La red binaria describe la estructura de acceso (qué aeropuertos pueden alcanzarse y con cuántas escalas), mientras que la ponderada describe la estructura de servicio (por dónde fluye la oferta). Las conclusiones sobre cuáles son los puentes críticos dependen, por tanto, de cuál de las dos se considere: Anchorage y Los Ángeles son críticos en ambas, mientras que los hubs europeos de mayor grado pierden relevancia cuando se toma en cuenta la intensidad del servicio.

11.3 Red con reciprocidad exigida

Finalmente, se construye una red que solo conserva las rutas operadas en ambos sentidos: dos aeropuertos están conectados si existen tanto la ruta \(A \to B\) como la \(B \to A\). Las rutas asimétricas pueden corresponder a vuelos triangulares (A → B → C → A), a rutas estacionales registradas en un solo sentido o a errores de registro, y no garantizan un viaje de ida y vuelta. Se analiza qué rutas desaparecen y cómo cambian la conectividad y la vulnerabilidad.

# ---- Red recíproca -------------------------------------------------------------------
g_mut_d <- subgraph.edges(g_d, E(g_d)[which_mutual(g_d)], delete.vertices = FALSE)
g_rec   <- as.undirected(g_mut_d, mode = "collapse", edge.attr.comb = "ignore")

clave <- function(g) {
  el <- as_edgelist(g, names = TRUE)
  paste(pmin(el[, 1], el[, 2]), pmax(el[, 1], el[, 2]), sep = "|")
}
perdida <- !(clave(giant_component) %in% clave(g_rec))       # rutas no recíprocas

# ---- Características de las rutas que desaparecen ----------------------------------------
el_u   <- as_edgelist(giant_component, names = FALSE)
cont_u <- V(giant_component)$continent; cont_u[is.na(cont_u)] <- "Otros"
inter_u <- cont_u[el_u[, 1]] != cont_u[el_u[, 2]]
km_u    <- E(giant_component)$distance_km
aer_u   <- E(giant_component)$num_aerolineas
eb_u    <- if (exists("eb_c")) eb_c else edge_betweenness(g_rob, directed = FALSE)

kable(
  data.frame(
    Grupo = c("Rutas no recíprocas (desaparecen)", "Rutas recíprocas (se conservan)"),
    Rutas = c(sum(perdida), sum(!perdida)),
    `Distancia mediana (km)` = c(median(km_u[perdida]), median(km_u[!perdida])),
    `% con una sola aerolínea` = c(100 * mean(aer_u[perdida] == 1), 100 * mean(aer_u[!perdida] == 1)),
    `% intercontinentales` = c(100 * mean(inter_u[perdida]), 100 * mean(inter_u[!perdida])),
    `Intermediación mediana (% de pares)` = c(100 * median(eb_u[perdida]) / choose(n0, 2),
                                              100 * median(eb_u[!perdida]) / choose(n0, 2)),
    check.names = FALSE, row.names = NULL
  ),
  digits = c(0, 0, 0, 1, 1, 4), format.args = list(decimal.mark = ",", big.mark = "."),
  caption = "Rutas que desaparecen al exigir reciprocidad frente a las que se conservan."
)
Rutas que desaparecen al exigir reciprocidad frente a las que se conservan.
Grupo Rutas Distancia mediana (km) % con una sola aerolínea % intercontinentales Intermediación mediana (% de pares)
Rutas no recíprocas (desaparecen) 806 1.009 86,8 22,5 0,0089
Rutas recíprocas (se conservan) 18.027 1.159 54,2 14,0 0,0059
top_perd <- which(perdida)[order(-eb_u[perdida])][1:8]
kable(
  data.frame(
    Ruta = paste(iata_u[el_u[top_perd, 1]], "–", iata_u[el_u[top_perd, 2]]),
    Ciudades = paste(V(giant_component)$city[el_u[top_perd, 1]], "–",
                     V(giant_component)$city[el_u[top_perd, 2]]),
    `Distancia (km)` = km_u[top_perd],
    Aerolíneas = aer_u[top_perd],
    `Intermediación (% de pares)` = 100 * eb_u[top_perd] / choose(n0, 2),
    check.names = FALSE, row.names = NULL
  ),
  digits = c(0, 0, 0, 0, 3), format.args = list(decimal.mark = ",", big.mark = "."),
  caption = "Rutas no recíprocas con mayor intermediación en la red original."
)
Rutas no recíprocas con mayor intermediación en la red original.
Ruta Ciudades Distancia (km) Aerolíneas Intermediación (% de pares)
BNE – DFW Brisbane – Dallas-Fort Worth 13.380 2 0,651
JED – THR Jeddah – Teheran 1.956 1 0,430
SYD – DFW Sydney – Dallas-Fort Worth 13.824 2 0,392
LPB – MIA La Paz – Miami 4.888 2 0,315
AMS – JRO Amsterdam – Kilimanjaro 6.909 3 0,296
BSB – MQH Brasilia – Minacu 260 1 0,251
AMS – SFB Amsterdam – Sanford 7.266 1 0,209
GRP – MQH Gurupi – Minacu 226 1 0,188
# ---- Conectividad y vulnerabilidad: red original vs. recíproca -------------------------------
m_orig <- metricas(g_rob)
m_rec  <- metricas(g_rec)

archivo_rec <- file.path("cache_resultados", "robustez_reciproca.rds")
if (file.exists(archivo_rec)) {
  rob_rec <- readRDS(archivo_rec)
} else {
  rob_rec <- list(
    grado = ataque_estatico(g_rec, V(g_rec)$name[order(-degree(g_rec))], pasos),
    betw  = ataque_estatico(g_rec, V(g_rec)$name[order(-betweenness(g_rec, directed = FALSE))], pasos)
  )
  saveRDS(rob_rec, archivo_rec)
}

kable(
  data.frame(
    Medida = c("Rutas", "Componentes", "Tamaño relativo de la componente gigante",
               "Eficiencia global relativa", "Alcance intercontinental",
               "Robustez R frente a ataque por grado", "Robustez R frente a ataque por intermediación"),
    `Red original` = c(fmt_num(ecount(g_rob), 0), fmt_num(m_orig[["Componentes"]], 0),
                       fmt_por(m_orig[["S"]], 1), fmt_por(m_orig[["E"]], 1), fmt_por(m_orig[["Inter"]], 1),
                       fmt_num(mean(rob$grado_est[, "S"]), 3), fmt_num(mean(rob$betw_est[, "S"]), 3)),
    `Red recíproca` = c(fmt_num(ecount(g_rec), 0), fmt_num(m_rec[["Componentes"]], 0),
                        fmt_por(m_rec[["S"]], 1), fmt_por(m_rec[["E"]], 1), fmt_por(m_rec[["Inter"]], 1),
                        fmt_num(mean(rob_rec$grado[, "S"]), 3), fmt_num(mean(rob_rec$betw[, "S"]), 3)),
    check.names = FALSE, row.names = NULL
  ),
  align = c("l", "r", "r"),
  caption = "Conectividad y vulnerabilidad de la red original frente a la red con reciprocidad exigida (medidas relativas a la red original)."
)
Conectividad y vulnerabilidad de la red original frente a la red con reciprocidad exigida (medidas relativas a la red original).
Medida Red original Red recíproca
Rutas 18.833 18.027
Componentes 1 114
Tamaño relativo de la componente gigante 100,0 % 96,4 %
Eficiencia global relativa 100,0 % 93,2 %
Alcance intercontinental 100,0 % 93,1 %
Robustez R frente a ataque por grado 0,733 0,704
Robustez R frente a ataque por intermediación 0,622 0,601

Al exigir reciprocidad desaparecen 806 rutas (4,3 % del total), que corresponden exactamente a los pares asimétricos de la red dirigida. Estas rutas tienen una distancia mediana de 1.009 km, frente a 1.159 km en las que se conservan, y el 86,8 % está operado por una sola aerolínea, lo que es coherente con rutas triangulares, estacionales o de baja frecuencia. Su pérdida tiene un efecto limitado y localizado: la red pasa de una a 114 componentes, pero la componente gigante conserva el 96,4 % de los aeropuertos, la eficiencia global el 93,2 % y el alcance intercontinental el 93,1 %. Los aeropuertos que quedan desconectados, alrededor del 3,6 %, son periféricos y dependían de una ruta en un solo sentido.

La vulnerabilidad del sistema tampoco cambia de forma sustancial. El índice de robustez de la red recíproca es ligeramente menor que el de la original, tanto frente al ataque por grado (0,704 frente a 0,733) como frente al ataque por intermediación (0,601 frente a 0,622), pero la diferencia, de unos 0,02–0,03, es del mismo orden que la pérdida inicial de aeropuertos al exigir reciprocidad, ya que todas las medidas se expresan respecto a la red original. Los ataques dirigidos degradan ambas redes al mismo ritmo, de modo que las conclusiones de la sección de vulnerabilidad se mantienen: la red sigue siendo robusta frente a fallas aleatorias y frágil frente a la pérdida de sus hubs y puentes.

12 Discusión

12.1 Estructura global y eficiencia de la red

La red aérea mundial es globalmente conectada pero localmente dispersa: la componente gigante reúne el 41,4 % de los aeropuertos, pero concentra el 99,9 % de las rutas, con una densidad de apenas 0,37 %. Pese a esa baja densidad, exhibe propiedades de mundo pequeño (\(L \approx 3,96\); el 90 % de los pares se conecta con 5 vuelos o menos), gracias a una distribución del grado muy asimétrica en la que unos pocos hubs concentran la mayor parte de las conexiones.

Esta eficiencia en número de vuelos se acompaña, en la mayoría de los casos, de eficiencia en kilómetros: en la mediana, el camino más corto sobre la red recorre solo un 10,5 % más que la línea recta, y el camino de menos escalas casi nunca se aleja mucho del más corto en distancia (penalización mediana de 0,4 %). El conflicto entre ambos criterios se concentra en un pequeño grupo de pares periféricos, sobre todo del Caribe y Centroamérica, donde minimizar escalas puede multiplicar la distancia hasta por 4,1, y donde aceptar una sola escala adicional suele bastar para volver al corredor geográfico natural.

12.2 Restricción espacial

La red está fuertemente condicionada por la geografía, y esta restricción no se explica por la secuencia de grados: al conservar exactamente los grados pero eliminar la información espacial, la probabilidad de conexión depende mucho menos de la distancia y las rutas se alargan más de cuatro veces en promedio (1.772 km frente a 7.504 km). La red concentra su conectividad en el rango regional, hasta unos 3.500 km, y solo el 14,3 % de sus rutas es intercontinental, frente al 72,0 % esperado sin restricción espacial. Los hubs mantienen rutas más largas que los aeropuertos pequeños (unos 1.845 km de media frente a 485 km en los de grado 1), un patrón opuesto al que produciría el azar, pero incluso ellos recorren solo cerca del 25 % de la longitud esperada sin restricción.

La geografía también organiza la estructura local y comunitaria de la red. El número de triángulos es 2,9 veces el esperado dados los grados, y la asortatividad por continente alcanza 0,800, frente a valores nulos cercanos a cero en ambos modelos de referencia: el 85,7 % de los extremos de las rutas une aeropuertos del mismo continente. Las comunidades detectadas por Louvain y Leiden reproducen la geografía a una escala más fina que los continentes, separando subregiones como Alaska, el Ártico canadiense, Groenlandia, el norte de Noruega o la región postsoviética, y las pocas comunidades transcontinentales corresponden a regiones de frontera (el Mediterráneo, Oriente Medio y el noreste de África, el espacio postsoviético) más que a conexiones de larga distancia. La proximidad geográfica y la topológica coinciden, por tanto, en gran medida: aeropuertos cercanos tienden a estar conectados, a compartir vecinos y a pertenecer a la misma comunidad, y la integración global descansa en un número reducido de rutas largas operadas desde los grandes hubs.

12.3 Hubs, puentes y núcleo de la red

La importancia de un aeropuerto no es unidimensional. Los rankings de centralidad coinciden en el panorama general, pero discrepan en el extremo superior, donde se distinguen tres perfiles: la conectividad y pertenencia al núcleo denso (grado, fuerza, coreness), dominada por los hubs europeos; la cercanía al resto de la red (armónica y propia), y el papel de puente (intermediación), que destaca a aeropuertos de tamaño medio como Anchorage, Papeete o los del norte de Canadá, que son la única entrada a subredes regionales. La participación continental precisa esta distinción: los hubs regionales articulan el tráfico dentro de su continente, mientras que los conectores intercontinentales concentran la intermediación global. Anchorage muestra, sin embargo, que un aeropuerto puede ser un puente crítico sin conectar continentes, si da acceso a una subregión aislada.

En el nivel global, la red combina una periferia en estrella con un núcleo muy cohesionado. El agrupamiento local cae con el grado (\(\hat\beta = -0,26\)), lo que indica una organización jerárquica en la que los hubs conectan grupos regionales que no se conectan entre sí. Al mismo tiempo, la asortatividad por grado es prácticamente nula pero muy superior a la que imponen los grados (-0,018 frente a -0,100), y los hubs tienen un 4 % más de rutas entre sí de lo esperado: existe un rich club moderado, coherente con el núcleo de coreness máxima (93 aeropuertos que requieren eliminar 31 para desconectarse) y con los clanes de 24 hubs norteamericanos completamente interconectados.

12.4 Vulnerabilidad

La red es robusta frente a fallas aleatorias y frágil frente a ataques dirigidos. Eliminar al azar el 10 % de los aeropuertos deja una componente gigante con el % de los aeropuertos originales, mientras que los ataques dirigidos la reducen a menos de la mitad y dividen la red en alrededor de mil componentes. El ataque más dañino es el adaptativo por intermediación (\(R = 0,437\) frente a 0,932 en las fallas aleatorias), que provoca un colapso abrupto cerca del 3 % de eliminaciones al ir eliminando los puentes que surgen tras cada pérdida. Esto confirma que una evaluación de la robustez basada solo en el número de conexiones subestima la criticidad de los nodos puente.

Los dos tipos de hubs se asocian a mecanismos de fragmentación distintos. La pérdida de hubs regionales desprende periferias y aumenta el número de componentes, mientras que la de conectores intercontinentales afecta sobre todo a la integración entre continentes: tras eliminar 220 hubs de cada tipo, el alcance intercontinental cae al 65,0 % en el primer caso y al 37,7 % en el segundo. El ataque por participación continental ilustra el caso extremo: la red conserva grandes bloques regionales que funcionan internamente, pero pierde casi toda la conexión entre continentes. La pérdida simultánea de todos los hubs europeos, en cambio, no resulta más dañina que una pérdida equivalente repartida por el mundo, lo que refleja la alta redundancia de la red europea.

A estas fuentes de vulnerabilidad se suman los puntos de articulación (329 aeropuertos, encabezados por Anchorage, cuya pérdida aísla 87 aeropuertos) y la baja redundancia operacional: más de la mitad de las rutas son operadas por una sola aerolínea, al exigir al menos dos operadores la componente mayor se reduce al 57,5 % de los aeropuertos, y 80 de las rutas de mayor intermediación dependen de un único operador.

12.5 Sensibilidad a la definición de la red

Las conclusiones no dependen de la dirección ni de la reciprocidad de las rutas. La red dirigida es casi simétrica (reciprocidad por arcos de 0,978), y la centralidad y las comunidades obtenidas con ella son prácticamente idénticas a las de la red no dirigida (NMI entre las particiones de Infomap de 0,977). Exigir reciprocidad elimina 806 rutas, en su mayoría operadas por una sola aerolínea, pero la componente gigante conserva el 96,4 % de los aeropuertos y la vulnerabilidad frente a ataques no cambia de forma sustancial.

La intensidad de servicio, en cambio, sí modifica la interpretación. Al ponderar por número de aerolíneas, las comunidades se reorganizan en torno a los grandes corredores comerciales, como el que une Australia con el sudeste asiático a través de Singapur, y los puentes más importantes pasan a ser aeropuertos con muchas aerolíneas, como Londres-Heathrow, Atlanta, Nueva York-JFK, Singapur y Miami. La red binaria describe, por tanto, la estructura de acceso, organizada por la geografía, mientras que la red ponderada describe la estructura de servicio, organizada por los corredores con competencia entre aerolíneas.

12.6 Limitaciones

  • Antigüedad de los datos. La base de rutas de OpenFlights dejó de actualizarse en junio de 2014, de modo que los resultados describen la red aérea de ese momento. Desde entonces han aparecido rutas relevantes para las conclusiones, como las conexiones directas entre Europa y Oceanía, y la pandemia de 2020 modificó profundamente la red.
  • Pesos aproximados. El número de aerolíneas y de registros ruta–aerolínea aproxima la intensidad del servicio, pero no corresponde al número real de vuelos, asientos o pasajeros; la red es, además, estática y no recoge la estacionalidad de las rutas.
  • Clasificación geográfica. La agrupación continental de countrycode trata a “Americas” como un único continente y asigna Turquía, Israel y Chipre a Asia y Rusia a Europa, lo que afecta a la participación continental, a la asortatividad y a la identificación de comunidades transcontinentales.
  • Modelos nulos y convenciones. El rewiring con grados preservados no conserva la ubicación de los hubs ni la estructura continental, por lo que no separa la restricción espacial de otros factores como la demanda o la regulación; como extensión, podría usarse un modelo de gravedad. La clasificación entre hubs regionales y conectores depende de umbrales convencionales (criterio de Tukey y mediana de \(P\)).
  • Simulaciones. Los ataques adaptativos recalculan la centralidad cada 10 eliminaciones y no después de cada una, y las fallas aleatorias y los modelos nulos usan un número limitado de réplicas (entre 20 y 1000), lo que fija la precisión de los intervalos y de los valores \(p\).

13 Referencias

  • Barthélemy, M. (2011). Spatial networks. Physics Reports, 499(1–3), 1–101.
  • Guimerà, R., Mossa, S., Turtschi, A. y Amaral, L. A. N. (2005). The worldwide air transportation network: Anomalous centrality, community structure, and cities’ global roles. Proceedings of the National Academy of Sciences, 102(22), 7794–7799.
  • Newman, M. E. J. (2018). Networks (2.ª ed.). Oxford University Press.
  • OpenFlights. Airport, airline and route data. https://openflights.org/data.html