# Clear the workspace
rm(list = ls())

# Load required libraries
library(bipartite)
library(readxl)
library(vegan)  # For beta diversity analysis if needed

# Read the matrix from the Excel file
file_path <- "data-METANETWORKS-QS.xlsx"  # Change to your file path
data <- read_excel(file_path, sheet = "aggregated")  # Change "aggregated" if needed

# If the first column contains row names
rownames(data) <- data[[1]]
data <- data[, -1]

# Compute observed network-level index
networklevel(data)

# Convert to matrix
data <- as.matrix(data)

# Calculate observed index value
observed_index <- networklevel(data, index = "connectance")  # Change metric if needed

# Generate null models using the "r2dtable" method
N <- 1000  # Number of null models to generate
null_models <- nullmodel(web = data, N = N, method = "r2d")  # Change null model method if needed

# Calculate the index for each null model
null_indices <- unlist(lapply(null_models, networklevel, index = "connectance"))

# Compute Z-score
null_mean <- mean(null_indices)  # Mean of null values
null_sd <- sd(null_indices)      # Standard deviation of null values
z_score <- (observed_index - null_mean) / null_sd

# Compute p-value
raw_p <- sum(null_indices > observed_index) / N
p_value <- ifelse(raw_p > 0.5, 1 - raw_p, raw_p)  # Two-tailed p-value

# Display results
cat("Observed index:", observed_index, "\n")
cat("Null mean:", null_mean, "\n")
cat("Null standard deviation:", null_sd, "\n")
cat("Z-score:", z_score, "\n")
cat("P-value:", p_value, "\n")
)