MetaWeave: rationale, methods, and simulation

Rationale

Species distributions tell us which interactions are possible in a place, but not which interactions are likely. A regional interaction model tells us which species may interact, but not where they co-occur. MetaWeave combines these two sources of information to reconstruct local ecological networks.

The package separates three questions:

1. Where are species present? Distribution layers define local assemblages.
2. How are co-occurring species connected? An inference model assigns interaction probabilities.
3. What might an observed sample look like? Simulation converts a probability network into interaction counts.

This separation allows the ecological model to change without rewriting the spatial workflow.

Core workflow

MetaWeave uses a small set of objects and transformations:

The main pipeline is:

distributions -> assemblages -> inference model -> local networks
              -> summaries -> spatial rasters

Inference methods

MetaWeave currently provides two general inference backends.

Simulated example

The following example creates a small landscape with six producer species and five consumer species. The distributions are artificial, so the full analysis is reproducible without external data.

library(metaweave)
library(terra)
#> terra 1.9.50

set.seed(42)

grid <- create_standard_grid(
  extent = c(0, 6, 0, 4),
  resolution = 1,
  crs = "EPSG:4326"
)

make_range <- function(template, center_x, center_y, spread) {
  xy <- terra::crds(template)
  distance <- sqrt(
    (xy[, 1] - center_x)^2 +
      (xy[, 2] - center_y)^2
  )

  layer <- terra::rast(template)
  terra::values(layer) <- ifelse(distance <= spread, 1, NA)
  layer
}

producer_parameters <- data.frame(
  species = paste0("producer_", 1:6),
  x = c(1.0, 2.0, 3.0, 4.0, 5.0, 3.0),
  y = c(1.0, 3.0, 2.0, 3.0, 1.0, 3.5),
  spread = c(2.0, 2.2, 1.8, 2.0, 1.8, 2.4)
)

consumer_parameters <- data.frame(
  species = paste0("consumer_", 1:5),
  x = c(1.5, 2.5, 4.0, 5.0, 3.0),
  y = c(2.5, 1.0, 2.0, 3.0, 3.5),
  spread = c(2.0, 1.8, 2.2, 1.8, 2.4)
)

producers <- terra::rast(lapply(seq_len(nrow(producer_parameters)), function(i) {
  make_range(
    grid,
    producer_parameters$x[i],
    producer_parameters$y[i],
    producer_parameters$spread[i]
  )
}))
names(producers) <- producer_parameters$species

consumers <- terra::rast(lapply(seq_len(nrow(consumer_parameters)), function(i) {
  make_range(
    grid,
    consumer_parameters$x[i],
    consumer_parameters$y[i],
    consumer_parameters$spread[i]
  )
}))
names(consumers) <- consumer_parameters$species

producer_lookup <- data.frame(
  species = names(producers),
  guild = rep(c("producer_A", "producer_B"), each = 3)
)

consumer_lookup <- data.frame(
  species = names(consumers),
  guild = c("consumer_A", "consumer_A", "consumer_B",
            "consumer_B", "consumer_B")
)

theta <- matrix(
  c(
    0.80, 0.20,
    0.25, 0.75
  ),
  nrow = 2,
  byrow = TRUE,
  dimnames = list(
    c("producer_A", "producer_B"),
    c("consumer_A", "consumer_B")
  )
)

model <- block_model(
  row_lookup = producer_lookup,
  column_lookup = consumer_lookup,
  theta = theta,
  row_group = "producers",
  column_group = "consumers"
)

result <- run_spatial_inference(
  distributions = list(
    producers = producers,
    consumers = consumers
  ),
  model = model,
  min_species = 1
)

head(result$result$index)
#>   cell_id   x   y row_richness column_richness possible_links expected_links
#> 1       1 0.5 3.5            1               1              1            0.8
#> 2       2 1.5 3.5            2               2              4            2.0
#> 3       3 2.5 3.5            4               3             12            5.9
#> 4       4 3.5 3.5            4               3             12            5.7
#> 5       5 4.5 3.5            2               3              6            4.5
#> 6       6 5.5 3.5            1               2              2            1.5
#>   mean_probability
#> 1        0.8000000
#> 2        0.5000000
#> 3        0.4916667
#> 4        0.4750000
#> 5        0.7500000
#> 6        0.7500000

plot(
  result$spatial$data[[
    c("row_richness", "column_richness", "expected_links")
  ]],
  main = c(
    "Producer richness",
    "Consumer richness",
    "Expected links"
  )
)

Simulating interaction observations

A probability network describes ecological expectation. Field data are finite samples from that expectation. simulate_interaction_counts() distributes a fixed number of interaction events among species pairs in proportion to their inferred probabilities.

We select the local network with the largest number of possible links.

focal_id <- which.max(result$result$index$possible_links)
focal_network <- result$result$networks[[focal_id]]

round(focal_network$matrix, 2)
#>            consumer_1 consumer_2 consumer_3 consumer_5
#> producer_2       0.80       0.80       0.20       0.20
#> producer_3       0.80       0.80       0.20       0.20
#> producer_4       0.25       0.25       0.75       0.75
#> producer_6       0.25       0.25       0.75       0.75

Simulate a small and a large survey:

counts_50 <- simulate_interaction_counts(
  focal_network,
  size = 50,
  seed = 1
)

counts_500 <- simulate_interaction_counts(
  focal_network,
  size = 500,
  seed = 1
)

sum(counts_50)
#> [1] 50
sum(counts_500)
#> [1] 500

Larger samples usually represent the underlying relative interaction weights more closely. We can demonstrate this across repeated simulations.

target <- focal_network$matrix / sum(focal_network$matrix)

sampling_error <- function(sample_size, replicates = 100) {
  errors <- vapply(seq_len(replicates), function(i) {
    observed <- simulate_interaction_counts(
      focal_network,
      size = sample_size,
      seed = 1000 + i
    )
    observed <- observed / sum(observed)
    mean(abs(observed - target))
  }, numeric(1))

  data.frame(
    sample_size = sample_size,
    mean_absolute_error = mean(errors),
    error_sd = stats::sd(errors)
  )
}

simulation_results <- do.call(
  rbind,
  lapply(c(25, 50, 100, 250, 500), sampling_error)
)

simulation_results
#>   sample_size mean_absolute_error    error_sd
#> 1          25          0.03950156 0.007059190
#> 2          50          0.02677969 0.005688817
#> 3         100          0.01905000 0.004034705
#> 4         250          0.01201438 0.002603380
#> 5         500          0.00829375 0.001767667

plot(
  simulation_results$sample_size,
  simulation_results$mean_absolute_error,
  type = "b",
  log = "x",
  xlab = "Number of sampled interactions",
  ylab = "Mean absolute error"
)

This analysis measures sampling error, not whether the inference model is biologically correct. Increasing sample size reduces observation noise but does not correct a misspecified interaction model.

Standardizing a network metric

Networks sampled with different effort are difficult to compare. One simple approach is to simulate each network at the same sample size and then calculate a metric. Here we use observed connectance.

observed_connectance <- function(count_matrix) {
  sum(count_matrix > 0) / length(count_matrix)
}

connectance_100 <- rarefied_network_metric(
  prob_mat = focal_network$matrix,
  target_size = 100,
  metric = observed_connectance,
  n_per_level = 100,
  seed = 99
)

connectance_100
#> [1] 0.973125

This value is the average connectance expected after 100 sampled interaction events. The same procedure can be applied to every local network or replaced by a specialist rarefaction method.

Using a species-level probability matrix

Some models return probabilities for individual species pairs rather than a single block partition. MetaWeave accepts these predictions directly.

species_probabilities <- theta[
  producer_lookup$guild,
  consumer_lookup$guild,
  drop = FALSE
]

dimnames(species_probabilities) <- list(
  producer_lookup$species,
  consumer_lookup$species
)

matrix_model <- probability_matrix_model(
  probability_matrix = species_probabilities,
  row_group = "producers",
  column_group = "consumers",
  name = "species-level predictions"
)

matrix_result <- run_spatial_inference(
  distributions = list(
    producers = producers,
    consumers = consumers
  ),
  model = matrix_model,
  min_species = 1
)

Because this example expands the same block model, its results match the block backend. In practice, the matrix could contain ensemble, trait-based, or machine-learning predictions.

Interpretation and limits

MetaWeave reconstructs plausible local networks conditional on its inputs. A predicted link is not a confirmed observation. Results inherit uncertainty and bias from:

• species distributions and their spatial resolution;
• incomplete regional interaction data;
• the chosen inference model and fitted parameters; and
• variation in sampling effort.

The current simulation functions represent finite sampling from a fixed local probability matrix. They do not yet propagate uncertainty in species ranges, block membership, model parameters, or environmental responses. Those sources should be handled through repeated model inputs or future uncertainty-aware backends.

The main value of MetaWeave is therefore comparative and modular: the same spatial assemblages can be evaluated under different ecological hypotheses, and their consequences can be compared across space using a shared workflow.