rm(list=ls())
library(readxl)
library(bipartite)
library(dplyr)

# Read the matrix from Excel (we use Sheet2 as you said)
matriz <- read_excel("data-METANETORKS-QS.xlsx", sheet = "parana") #change for others

# If it has row names in the first column
rownames(matriz) <- matriz[[1]]
matriz <- matriz[,-1]

# Convert to numeric matrix
numerica <- as.matrix(matriz)

# Run DIRT_LPA_wb_plus 100 times
max_modularity <- -Inf
best_result <- NULL

for (i in 1:100) {
  res <- DIRT_LPA_wb_plus(numerica)
  modularity_value <- res$modularity
  
  if (modularity_value > max_modularity) {
    max_modularity <- modularity_value
    best_result <- res
  }
}

# Convert the result to a moduleWeb object
modulos <- convert2moduleWeb(numerica, best_result)

# Plot the modular network (optional, but may fail if there's something wrong with the object)
# If you get an error here, you can ignore it for now
suppressWarnings(plotModuleWeb(modulos))

# Extract c and z values
cz_values <- czvalues(modulos, weighted = TRUE)
cz_table <- cz_values$cz

# Check
head(cz_table)

cz_table <- data.frame(
  interaction = names(cz_values$c),
  c = cz_values$c,
  z = cz_values$z
)
head(cz_table)
library(openxlsx)

# Create data frame with c and z values
cz_df <- data.frame(
  Interaccion = names(cz_values$c),
  C = cz_values$c,
  Z = cz_values$z
)

# View the table in RStudio (optional)
View(cz_df)

# Export as Excel
write.xlsx(cz_df, "cz_values_interaccionesparana.xlsx", row.names = FALSE)