在R中用Monte Carlo/Bootstrapping实现网络抽样与平行边权重分析

Got it, let's walk through how to solve this problem with igraph and base R—no extra packages required unless you want to streamline some steps. I’ve tackled similar network sampling tasks before, so here’s a practical, reproducible approach:

Step 1: Prepare a Test Network (for validation) #

First, let's simulate a large multi-graph (with parallel edges) to test our code. This mirrors the structure you’re working with:

library(igraph)

set.seed(123) # For reproducibility
large_graph <- make_empty_graph(directed = FALSE) %>%
  add_vertices(1000) %>%
  # Add intentional parallel edges between specific node pairs
  add_edges(c(
    rep(c(1,2), 5),  # 5 parallel edges between 1 & 2
    rep(c(3,4), 3),  # 3 parallel edges between 3 & 4
    rep(c(5,6), 7),  # 7 parallel edges between 5 & 6
    # Add random non-parallel edges to mimic a real large network
    sample(1:1000, 2000, replace = TRUE)
  )) %>%
  # Assign random weights to edges
  set_edge_attr("weight", value = rnorm(ecount(.), mean = 5, sd = 2))

Step 2: Define a Sampling & Calculation Function #

We’ll create a reusable function that handles one round of sampling, extracts the subgraph, and computes the standard deviation of parallel edge weights. This function accounts for both directed and undirected graphs:

sample_and_calculate_sd <- function(graph, sample_size, directed = FALSE) {
  # Convert proportion to node count if needed
  if (sample_size < 1) {
    sample_size <- round(sample_size * vcount(graph))
  }
  
  # Randomly sample nodes (without replacement)
  sampled_nodes <- sample(V(graph), size = sample_size, replace = FALSE)
  
  # Extract the induced subgraph from sampled nodes
  subgraph <- induced_subgraph(graph, sampled_nodes)
  
  # Create unique identifiers for edge pairs
  if (directed) {
    # For directed graphs, (1->2) != (2->1)
    edge_pairs <- paste(as.integer(head_of(subgraph, E(subgraph))),
                        as.integer(tail_of(subgraph, E(subgraph))),
                        sep = "-")
  } else {
    # For undirected graphs, standardize pair order to avoid duplicate groups
    edge_pairs <- apply(cbind(
      as.integer(head_of(subgraph, E(subgraph))),
      as.integer(tail_of(subgraph, E(subgraph)))
    ), 1, function(x) paste(sort(x), collapse = "-"))
  }
  
  # Extract edge weights and compute SD for each parallel edge pair
  weights <- E(subgraph)$weight
  pair_sds <- tapply(weights, edge_pairs, sd, na.rm = TRUE)
  
  # Filter out single-edge pairs (SD is NA for these)
  pair_sds <- pair_sds[!is.na(pair_sds)]
  
  # Return both pair-level SDs and an overall SD (to measure global fluctuation)
  list(
    pair_level_sds = pair_sds,
    overall_sd = ifelse(length(pair_sds) > 0, sd(pair_sds), NA)
  )
}

Step 3: Run Multiple Samples & Analyze Fluctuation #

Now we’ll run the function multiple times, collect results, and analyze how the parallel edge weight SDs vary across samples:

# Configure sampling parameters
n_samples <- 100          # Number of random samples to run
sample_proportion <- 0.2  # Sample 20% of nodes each time

# Run all samples and store results
results <- replicate(n_samples,
                     sample_and_calculate_sd(large_graph, sample_proportion),
                     simplify = FALSE)

# Extract the overall SD from each sample for analysis
overall_sds <- sapply(results, function(x) x$overall_sd)

# Visualize the distribution of overall SDs
hist(overall_sds,
     main = "Distribution of Parallel Edge Weight SDs Across Samples",
     xlab = "Overall SD of Parallel Edge Weights",
     col = "lightblue")

# Compute summary statistics for the fluctuation
cat("Mean of overall SDs across samples:", mean(overall_sds, na.rm = TRUE), "\n")
cat("SD of overall SDs (measures fluctuation):", sd(overall_sds, na.rm = TRUE), "\n")

Additional Notes #

• Edge Sampling Instead of Node Sampling: If you want to sample edges directly instead of nodes, replace the node sampling step with sampled_edges <- sample(E(graph), size = edge_sample_size) and use subgraph.edges(graph, sampled_edges) to extract the subgraph.
• Handling Large Networks: For extremely large graphs, induced_subgraph is efficient, but if memory is an issue, you can use delete_vertices(graph, V(graph)[-sampled_nodes]) as an alternative.
• Missing Weights: If your network has missing weight values, ensure na.rm = TRUE is included in the sd call to avoid invalid results.

内容的提问来源于stack exchange，提问作者FilipeTeixeira