同图绘制多个 bipartite 网络度分布（含幂律线）及批量分析需求

Hey there! Let's work through your bipartite network analysis needs step by step. I'll split the solution into two core parts: creating that log-log degree distribution plot with multiple networks and a scale-free comparison, plus calculating and exporting the Pearson correlations and slopes for all 500 networks to Excel.

1. Log-Log Degree Distribution Plot (10+ Networks + Scale-Free Comparison) #

First, we'll use igraph for network handling and ggplot2 for visualization—they play nicely together for this task. Let's start by loading the required packages and confirming your network generation code works (I'll reuse your exact seed to keep consistency):

library(igraph)
library(ggplot2)
library(dplyr)
library(writexl) # For Excel export later

Step 1: Generate your 500 bipartite networks (your code, unchanged) #

set.seed(1) 
gs1 <- list() 
for (x in seq_len(500L)) { 
  gs1[[x]] <- sample_bipartite(358, 27, type = "gnm", m = 827, directed = TRUE) 
}

Step 2: Helper function to extract log-transformed degree distributions #

This function takes a single network, calculates its degree distribution, and converts values to log10 (since we need a log-log plot):

get_log_degree_dist <- function(g) {
  # Calculate node degrees
  node_degrees <- degree(g)
  # Convert to a frequency table, then clean and log-transform
  dist_df <- as.data.frame(table(node_degrees)) %>%
    mutate(
      degree = as.integer(as.character(node_degrees)),
      frequency = as.integer(Freq),
      log_degree = log10(degree),
      log_frequency = log10(frequency)
    ) %>%
    filter(log_degree != -Inf, log_frequency != -Inf) # Drop zero values that break log transforms
  return(dist_df)
}

Step 3: Extract data for 10 networks + a scale-free comparison network #

Let's pick the first 10 networks (you can swap this for a random sample with sample(gs1, 10) if you prefer) and add a Barabási-Albert scale-free network for comparison:

# Grab 10 networks from your list
selected_networks <- gs1[1:10]
# Batch extract their degree distributions and label each network
network_distributions <- lapply(seq_along(selected_networks), function(net_id) {
  dist_data <- get_log_degree_dist(selected_networks[[net_id]])
  dist_data$network_label <- paste0("Bipartite Net ", net_id)
  return(dist_data)
})
# Combine into one data frame
all_dist_data <- do.call(rbind, network_distributions)

# Add a scale-free network for comparison
scale_free_net <- barabasi.game(n = vcount(selected_networks[[1]]), m = 2, directed = TRUE)
sf_dist_data <- get_log_degree_dist(scale_free_net)
sf_dist_data$network_label <- "Scale-Free (BA)"
all_dist_data <- rbind(all_dist_data, sf_dist_data)

Step 4: Plot the log-log degree distributions with power-law lines #

This plot will show all 10 bipartite networks, the scale-free comparison, a dashed linear fit (power-law approximation) for each, and a solid black reference power-law line (slope = -2, a common value for scale-free networks):

ggplot(all_dist_data, aes(x = log_degree, y = log_frequency, color = network_label)) +
  geom_point(alpha = 0.6, size = 1.5) + # Plot raw data points
  geom_smooth(method = "lm", se = FALSE, linetype = "dashed") + # Add linear fit (power-law line) for each network
  geom_abline(slope = -2, intercept = max(all_dist_data$log_frequency), 
              color = "black", linetype = "solid", size = 1) + # Reference power-law line
  labs(
    x = "log₁₀(Degree)",
    y = "log₁₀(Frequency)",
    title = "Log-Log Degree Distributions: Bipartite vs Scale-Free Networks",
    color = "Network"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

2. Batch Calculate Pearson Correlations & Slopes for All 500 Networks #

Next, we'll compute the Pearson correlation between log(degree) and log(frequency), plus the slope of the linear fit (which is the power-law exponent) for every network, then export to Excel.

Step 1: Helper function to compute stats for a single network #

compute_network_stats <- function(g) {
  dist_data <- get_log_degree_dist(g)
  # Skip networks with too few data points to fit a line
  if(nrow(dist_data) < 2) {
    return(data.frame(pearson_cor = NA, slope = NA))
  }
  # Fit linear model to log-log data
  power_law_fit <- lm(log_frequency ~ log_degree, data = dist_data)
  # Extract slope and Pearson correlation
  slope_val <- coef(power_law_fit)[["log_degree"]]
  pearson_val <- cor(dist_data$log_degree, dist_data$log_frequency, method = "pearson")
  return(data.frame(pearson_cor = pearson_val, slope = slope_val))
}

Step 2: Batch compute stats for all 500 networks #

# Run the function on every network in your list
all_network_stats <- lapply(gs1, compute_network_stats)
# Combine into a single data frame and add network IDs
stats_df <- do.call(rbind, all_network_stats)
stats_df$network_id <- seq_len(500)
# Reorder columns for readability
stats_df <- stats_df[, c("network_id", "pearson_cor", "slope")]

Step 3: Export to Excel #

write_xlsx(stats_df, "bipartite_network_degree_stats.xlsx")

Quick Notes #

• If some networks have very sparse degree distributions (only 1 unique degree value), the code returns NA for their stats to avoid errors—you can filter these out later if needed.
• The reference power-law line (slope = -2) is a common benchmark for scale-free networks; compare your bipartite network slopes to see how they differ.
• Feel free to adjust the number of networks plotted, the reference slope, or sampling method to fit your needs.

内容的提问来源于stack exchange，提问作者The Last Word