R Notebook

# 1. Install required packages
if (!require("mgm")) install.packages("mgm")

Loading required package: mgm
Warning: there is no package called ‘mgm’WARNING: Rtools is required to build R packages but is not currently installed. Please download and install the appropriate version of Rtools before proceeding:

https://cran.rstudio.com/bin/windows/Rtools/
Installing package into ‘C:/Users/user/AppData/Local/R/win-library/4.4’
(as ‘lib’ is unspecified)
trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.4/mgm_1.2-15.zip'
Content type 'application/zip' length 924996 bytes (903 KB)
downloaded 903 KB

package ‘mgm’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\user\AppData\Local\Temp\RtmpcXsMWX\downloaded_packages

if (!require("qgraph")) install.packages("qgraph")

Loading required package: qgraph
Warning: package ‘qgraph’ was built under R version 4.4.3Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

if (!require("igraph")) install.packages("igraph")

Loading required package: igraph
Warning: package ‘igraph’ was built under R version 4.4.3
Attaching package: ‘igraph’

The following objects are masked from ‘package:stats’:

    decompose, spectrum

The following object is masked from ‘package:base’:

    union

library(mgm)

Warning: package ‘mgm’ was built under R version 4.4.3This is mgm 1.2-15
Please report issues on Github: https://github.com/jmbh/mgm/issues

library(qgraph)
library(igraph)

library(readr)
library(dplyr)

# =========================================================================
# STEP 1.1: Load the CSV Dataset
# =========================================================================
file_path <- "C:/Users/user/Downloads/GE2324/Project/Dataset- With Variable context explained/06399-Offenders-Data.csv"
df_raw <- read_csv(file_path, show_col_types = FALSE)

# =========================================================================
# STEP 1.2: Robust Cleaning & Handling Structural Missingness
# =========================================================================
clean_offender_csv_data <- function(df) {
  # 1. Recode standard ICPSR missing/confidentiality codes to NA
  df_recoded <- df %>%
    mutate(
      across(where(is.numeric), ~na_if(., 97)),
      across(where(is.numeric), ~na_if(., 98)),
      across(where(is.numeric), ~na_if(., 99)),
      across(where(is.numeric), ~na_if(., 9998)),
      across(where(is.numeric), ~na_if(., 9999))
    )
  
  # 2. Identify and drop columns that are 100% missing (e.g., confidentiality-blanked dates)
  missing_pct <- sapply(df_recoded, function(x) mean(is.na(x)))
  cols_to_drop <- names(missing_pct[missing_pct == 1])
  
  cat("Dropping 100% missing/blanked columns:\n")
  print(cols_to_drop)
  
  df_filtered_cols <- df_recoded %>%
    select(-all_of(cols_to_drop))
  
  # 3. Construct chronological time index from BOOKYEAR for tvmgm() timepoints
  df_with_time <- df_filtered_cols %>%
    mutate(
      time_index = as.numeric(BOOKYEAR) - min(as.numeric(BOOKYEAR), na.rm = TRUE) + 1
    )
  
  # 4. Retain complete cases on the remaining operational columns
  df_final <- df_with_time %>%
    filter(complete.cases(.))
  
  return(df_final)
}

df_clean <- clean_offender_csv_data(df_raw)

Dropping 100% missing/blanked columns:
[1] "INJDTE"   "DEATHDTE"

# =========================================================================
# STEP 1.3: Extract Feature Matrix & Define Variable Taxonomy for mgm
# =========================================================================
feature_df <- df_clean %>%
  select(-any_of(c("HOMINEW", "OFFND", "time_index")))

feature_matrix <- as.matrix(feature_df)
p <- ncol(feature_matrix)

# Define variable types vector ("g" = Gaussian, "c" = Categorical, "p" = Poisson)
var_types <- rep("g", p)
var_levels <- rep(1, p)

time_vector <- df_clean$time_index

# Verification checks
cat("\nSuccessfully preprocessed", nrow(feature_matrix), "complete observations across", p, "variables.\n")

Successfully preprocessed 10845 complete observations across 74 variables.

cat("Dimensions of feature_matrix:", dim(feature_matrix), "\n")

Dimensions of feature_matrix: 10845 74 

cat("Length of time_vector:", length(time_vector), "\n")

Length of time_vector: 10845 

Ready for proceeding to tvmgm().

Why This Output Is Correct

1. Resolution of Structural Missingness: The script successfully identified and dropped the two 100% blanked/confidentiality columns (INJDTE and DEATHDTE), avoiding the previous fatal error where 100% of the dataset was filtered out.
2. Accurate Feature Dimensionality Math: The feature matrix correctly contains 74 variables, which matches the exact column accounting:

• 78 original columns in the raw dataset
• - 2 dropped structural columns (INJDTE, DEATHDTE)
• - 2 excluded identifier columns (HOMINEW, OFFND)
• = 74 operational feature variables

3. Satisfaction of mgm Constraints: The mgm package strictly requires complete data with zero missing values (NA). Retaining 10,845 complete observations means the filtering logic successfully cleaned the data while preserving a large, robust sample size for estimation.
4. Temporal Alignment: The time_vector length of 10,845 matches the row count of feature_matrix (\(10,845 \times 74\)), ensuring seamless synchronization when passing time points and features into bwSelect() and tvmgm().

———————

Phase 2: Time-feature extraction+Specifying Variable Taxonomy Before running the time-varying mixed graphical model, we must extract all the time features, and define the data types and levels for each variable in your dataset so the algorithm knows whether to treat them as continuous, categorical, or count data.

# =========================================================================
# STEP 2.4: Chronological Sorting & Continuous Time Vector Creation
# =========================================================================

# 1. Sort the data chronologically to satisfy tvmgm assumptions
df_clean <- df_clean %>%
  arrange(BOOKYEAR, INJMONTH, INJTIME)

# 2. Create a continuous, strictly non-decreasing time vector
df_clean <- df_clean %>%
  mutate(
    # Normalize month (1-12) to a fraction of the year
    month_frac = ifelse(is.na(INJMONTH), 0, (INJMONTH - 1) / 12),
    # Normalize military time (0-2400) to a fraction of a month
    time_frac = ifelse(is.na(INJTIME), 0, (INJTIME / 2400) * (1/12)),
    # Combine into a smooth temporal index
    continuous_time = BOOKYEAR + month_frac + time_frac
  )

# 3. Extract the final ordered time vector and normalize it to [0,1]
time_vector <- df_clean$continuous_time
time_vector <- (time_vector - min(time_vector)) / (max(time_vector) - min(time_vector))


# =========================================================================
# STEP 2.5 (THE TEMPORAL DENSITY FILTER): Catching Local Extinction
# =========================================================================

# 1. Drop Time Variables, Ultra-High Cardinality Variables, AND our new time construction vars!
drop_vars <- c(
  "BOOKYEAR", "INJYEAR", "INJMONTH", "INJDAY", "DEATHYR", "DEATHMON", 
  "CENTRACT", "COMAREA", "DISTRICT", "LOCATION", 
  "CAUSFACT", "CAUSFAC2", "VREL", "OREL", "WCLUB", "WOTHER", "WHANDGUN", "WKNIFE",
  "HOMINEW", "OFFND", "time_index", "month_frac", "time_frac", "continuous_time" # Newly added drops
)

feature_df_subset <- df_clean %>%
  select(-any_of(drop_vars))

p_subset <- ncol(feature_df_subset)
actual_cols_subset <- colnames(feature_df_subset)

# 2. Re-assign Types for the Subset
var_types_subset <- rep("c", p_subset)
var_types_subset[actual_cols_subset %in% c("INJTIME", "INJ2DTH", "DEATHTIM")] <- "g"
var_types_subset[actual_cols_subset %in% c("NUMVIC", "NUMOFF")] <- "p"

feature_matrix_subset <- as.matrix(feature_df_subset)
N_ROWS <- nrow(feature_matrix_subset)

# 3. Proportional Lumping (Global)
MIN_PROP <- 0.05  
for (i in seq_len(p_subset)) {
  if (var_types_subset[i] == "c") {
    repeat {
      cat_counts <- table(feature_matrix_subset[, i])
      cat_props <- cat_counts / N_ROWS
      rare_cats <- as.numeric(names(cat_props[cat_props < MIN_PROP]))
      
      if (length(rare_cats) == 0) break 
      
      mode_cat <- as.numeric(names(sort(cat_counts, decreasing = TRUE)[1]))
      is_too_rare <- feature_matrix_subset[, i] %in% rare_cats
      feature_matrix_subset[is_too_rare, i] <- mode_cat
    }
    feature_matrix_subset[, i] <- as.numeric(as.factor(feature_matrix_subset[, i]))
  }
}

# 4. TEMPORAL DENSITY CHECK (The glmnet 10^-5 Fix)
# We simulate the exact tvmgm Gaussian kernel to ensure no category vanishes locally.
time_norm <- time_vector # Already normalized to [0,1] in Step 2.4 above
estpoints_test <- seq(0, 1, length = 30)
bw_test <- 0.1  # We test the narrowest bandwidth from our bwSeq

var_levels_subset <- numeric(p_subset)
cols_to_keep <- logical(p_subset)

for (i in seq_len(p_subset)) {
  if (var_types_subset[i] == "c") {
    var_levels_subset[i] <- length(unique(feature_matrix_subset[, i]))
    
    # Must have >= 2 levels globally
    if(var_levels_subset[i] < 2) {
      cols_to_keep[i] <- FALSE
      next
    }
    
    cats <- unique(feature_matrix_subset[, i])
    is_temporally_safe <- TRUE
    
    # Check probability at every local estimation window
    for (ep in estpoints_test) {
      # Replicate the tvmgm Equation 13 Gaussian Kernel
      w <- exp(-0.5 * ((time_norm - ep) / bw_test)^2)
      w <- w / max(w) 
      total_weight_sum <- sum(w)
      
      for (c_val in cats) {
        cat_weight_sum <- sum(w[feature_matrix_subset[, i] == c_val])
        
        # If the weighted probability drops below glmnet's threshold, flag for removal
        if ((cat_weight_sum / total_weight_sum) <= 1e-4) {
          is_temporally_safe <- FALSE
          break
        }
      }
      if (!is_temporally_safe) break
    }
    
    cols_to_keep[i] <- is_temporally_safe
    
  } else {
    var_levels_subset[i] <- 1
    # Check that continuous/Poisson vars have variance
    cols_to_keep[i] <- var(feature_matrix_subset[, i], na.rm = TRUE) > 0
  }
}

# 5. Create Final Guaranteed-Safe Data Objects
feature_matrix_final <- feature_matrix_subset[, cols_to_keep]
var_types_final <- var_types_subset[cols_to_keep]
var_levels_final <- var_levels_subset[cols_to_keep]
actual_cols_final <- actual_cols_subset[cols_to_keep]

cat("Temporal Density Check Complete.\n")

Temporal Density Check Complete.

cat("Retained", sum(cols_to_keep), "temporally stable features out of", p_subset, "for tvmgm.\n")

Retained 37 temporally stable features out of 56 for tvmgm.

# Fix Step 3 compute long
# =========================================================================
# STEP 2.6: The Missing Standardization 
# =========================================================================

cat("Applying standard normal scaling to continuous variables...\n")

g_indices <- which(var_types_final == "g")

for (i in g_indices) {
  # Force as numeric and scale to mean=0, sd=1
  feature_matrix_final[, i] <- scale(as.numeric(feature_matrix_final[, i]))
}

cat("Standardization complete. Matrix is mathematically flat.\n")

Step : Estimating the Time-Varying Mixed Graphical Model

This step first performs a time-stratified cross-validation to select the optimal bandwidth for the Gaussian kernel smoothing, and then fits the sequence of local models across the 30-year span.

The following Does not compute Optimal Bandwidth from scratch, which is computationally demanding. A substitue of 0.2 is chosen

#Fast track version of Step 3
# =========================================================================
# STEP 3.2 (FAST-TRACK RESTART): Fit the tvmgm Model 
# =========================================================================
set.seed(123)

optimal_bw <- 0.2  

cat("Fitting Fast-Track tvmgm Model (10 Estimation Points)...\n")

fit_tvmgm <- tvmgm(
  data = feature_matrix_final,
  type = var_types_final,
  level = var_levels_final,
  timepoints = time_vector,
  estpoints = seq(0, 1, length = 10), 
  k = 2,
  bandwidth = optimal_bw,        
  threshold = "none",
  ruleReg = "AND",
  pbar = TRUE
)

cat("Estimation Complete!\n")

But original flow involves so

# =========================================================================
# STEP 3.1 (OPTIMIZED): Select Optimal Bandwidth 
# =========================================================================
set.seed(123)

cat("Starting bandwidth selection. Progress bar will appear shortly...\n")

bw_tvmgm <- bwSelect(
  data = feature_matrix_final,      
  type = var_types_final,           
  level = var_levels_final,         
  bwSeq = c(0.1, 0.2),              # Reduced search grid for faster execution
  bwFolds = 3,                      # Reduced to 3-fold CV to speed up computation
  bwFoldsize = 3,
  timepoints = time_vector,         
  modeltype = "mgm", 
  k = 2,                            
  threshold = "none",            
  ruleReg = "AND",                  
  pbar = TRUE
)

# Inspect the prediction error
cat("\nMean Errors for candidate bandwidths:\n")
print(round(bw_tvmgm$meanError, 3))

optimal_bw <- bw_tvmgm$bwSeq[which.min(bw_tvmgm$meanError)]
cat("\nOptimal bandwidth selected:", optimal_bw, "\n")

Updated version for csv output

# =========================================================================
# STEP 4: Compute Predictions & Nodewise Errors (Corrected Tabular CSV Export)
# =========================================================================
cat("Computing predictions and nodewise errors over time...\n")

Computing predictions and nodewise errors over time...

pred_tvmgm <- predict(
  object = fit_tvmgm,
  data = feature_matrix_final,
  tvMethod = "weighted",        # Computes a weighted average across time
  errorCon = c("RMSE", "R2"),   # Error metrics for Gaussian/Poisson
  errorCat = c("CC", "nCC")     # Error metrics for Categorical (Accuracy)
)

cat("Predictions extracted successfully!\n")

Predictions extracted successfully!

# 1. Export Aggregate Nodewise Errors (Overall Time-Series Summary)
error_matrix <- as.data.frame(pred_tvmgm$errors)
# Insert variable names cleanly
error_matrix$Variable <- actual_cols_final

output_file_aggregate <- "tvmgm_nodewise_aggregate_errors.csv"
write.csv(error_matrix, output_file_aggregate, row.names = FALSE)
cat("Successfully exported aggregate nodewise errors to:", output_file_aggregate, "\n")

Successfully exported aggregate nodewise errors to: tvmgm_nodewise_aggregate_errors.csv 

# 2. Export Time-Specific Errors Across Each Estimation Point (from tverrors)
if (!is.null(pred_tvmgm$tverrors)) {
  tverrors_list <- pred_tvmgm$tverrors
  tverror_dfs <- list()
  
  for (tp in seq_along(tverrors_list)) {
    tp_data <- as.data.frame(tverrors_list[[tp]])
    tp_data$Estimation_Point <- tp
    tp_data$Variable <- actual_cols_final
    tverror_dfs[[tp]] <- tp_data
  }
  
  master_tverror_table <- bind_rows(tverror_dfs)
  output_file_tverrors <- "tvmgm_nodewise_errors_per_estimation_point.csv"
  write.csv(master_tverror_table, output_file_tverrors, row.names = FALSE)
  cat("Successfully exported time-specific nodewise errors to:", output_file_tverrors, "\n")
}

Successfully exported time-specific nodewise errors to: tvmgm_nodewise_errors_per_estimation_point.csv 

# =========================================================================
# STEP 5: Visualizing the Time-Varying Networks
# =========================================================================
library(qgraph)

# Extract the final sanitized variable names
var_names <- actual_cols_final

# Set up the plotting window for side-by-side comparison
par(mfrow = c(1, 2))

# Plot 1: Early Era (Estimation Point 2 out of 10)
cat("Rendering Early Era Network...\n")
qgraph(
  fit_tvmgm$pairwise$wadj[, , 2], 
  edge.color = fit_tvmgm$pairwise$edgecolor[, , 2],
  layout = "spring",
  labels = var_names,
  theme = "colorblind",
  vsize = 5,
  cut = 0.1,  # Hides very weak edges to declutter the visual
  title = "Chicago Homicide Network\nEarly Era (Est. Point 2)"
)

# Plot 2: Late Era (Estimation Point 9 out of 10)
cat("Rendering Late Era Network...\n")
qgraph(
  fit_tvmgm$pairwise$wadj[, , 9], 
  edge.color = fit_tvmgm$pairwise$edgecolor[, , 9],
  layout = "spring",
  labels = var_names,
  theme = "colorblind",
  vsize = 5,
  cut = 0.1,
  title = "Chicago Homicide Network\nLate Era (Est. Point 9)"
)

# Reset plotting window back to normal
par(mfrow = c(1, 1))

Updated version for Network description

# =========================================================================
# STEP 5: Visualizing the Time-Varying Networks & Exporting Edge Tables
# =========================================================================
library(qgraph)
library(dplyr)

# Extract the final sanitized variable names
var_names <- actual_cols_final
n_estpoints <- dim(fit_tvmgm$pairwise$wadj)[3]

# -------------------------------------------------------------------------
# Tabular Extraction: Export Edge Lists Across All Estimation Points
# -------------------------------------------------------------------------
edge_list_list <- list()

for (tp in 1:n_estpoints) {
  wadj_mat <- fit_tvmgm$pairwise$wadj[, , tp]
  rownames(wadj_mat) <- var_names
  colnames(wadj_mat) <- var_names
  
  for (i in 1:(nrow(wadj_mat) - 1)) {
    for (j in (i + 1):ncol(wadj_mat)) {
      w <- wadj_mat[i, j]
      if (!is.na(w) && w != 0) {
        edge_list_list[[length(edge_list_list) + 1]] <- data.frame(
          Estimation_Point = tp,
          Node_1           = var_names[i],
          Node_2           = var_names[j],
          Weight           = w,
          Abs_Weight       = abs(w)
        )
      }
    }
  }
}

master_edge_table <- bind_rows(edge_list_list)
output_file_edges <- "tvmgm_pairwise_edges_all_timepoints.csv"
write.csv(master_edge_table, output_file_edges, row.names = FALSE)
cat("Successfully exported master edge list table to:", output_file_edges, "\n")

Successfully exported master edge list table to: tvmgm_pairwise_edges_all_timepoints.csv 

# -------------------------------------------------------------------------
# Visualization: Side-by-Side qgraph Rendering
# -------------------------------------------------------------------------
par(mfrow = c(1, 2))

# Plot 1: Early Era (Estimation Point 2 out of 10)
cat("Rendering Early Era Network...\n")

Rendering Early Era Network...

qgraph(
  fit_tvmgm$pairwise$wadj[, , 2], 
  edge.color = fit_tvmgm$pairwise$edgecolor[, , 2],
  layout = "spring",
  labels = var_names,
  theme = "colorblind",
  vsize = 5,
  cut = 0.1,  
  title = "Chicago Homicide Network\nEarly Era (Est. Point 2)"
)

# Plot 2: Late Era (Estimation Point 9 out of 10)
cat("Rendering Late Era Network...\n")

Rendering Late Era Network...

qgraph(
  fit_tvmgm$pairwise$wadj[, , 9], 
  edge.color = fit_tvmgm$pairwise$edgecolor[, , 9],
  layout = "spring",
  labels = var_names,
  theme = "colorblind",
  vsize = 5,
  cut = 0.1,
  title = "Chicago Homicide Network\nLate Era (Est. Point 9)"
)

par(mfrow = c(1, 1))

---

Network Statistics

Tabular extraction of the Network Results

library(igraph)
library(dplyr)

# Extract variable names
var_names <- actual_cols_final
n_estpoints <- dim(fit_tvmgm$pairwise$wadj)[3]

# Initialize a master list to store metrics for all time points
network_metrics_list <- list()

for (tp in 1:n_estpoints) {
  # 1. Extract weighted adjacency matrix for estimation point `tp`
  wadj_matrix <- fit_tvmgm$pairwise$wadj[, , tp]
  rownames(wadj_matrix) <- var_names
  colnames(wadj_matrix) <- var_names
  
  # 2. Convert to an igraph object (treating weights as absolute values)
  g <- graph_from_adjacency_matrix(
    abs(wadj_matrix), 
    mode = "undirected", 
    weighted = TRUE, 
    diag = FALSE
  )
  
  # 3. Compute Centrality & Network Statistics
  node_stats <- data.frame(
    Estimation_Point = tp,
    Variable         = var_names,
    Degree           = degree(g),
    Strength         = strength(g),                    # Sum of absolute edge weights
    Betweenness      = betweenness(g, normalized = TRUE),
    Closeness        = closeness(g, normalized = TRUE)
  )
  
  network_metrics_list[[tp]] <- node_stats
}

# Combine all estimation points into one comprehensive master table
master_network_table <- bind_rows(network_metrics_list)

# View top most central nodes across the timeline
print("Top 10 Most Central Variables by Strength:")

[1] "Top 10 Most Central Variables by Strength:"

print(
  master_network_table %>% 
    arrange(desc(Strength)) %>% 
    head(10)
)

NA

# Select a specific time slice, e.g., Estimation Point 9 (Late Era)
wadj_late <- fit_tvmgm$pairwise$wadj[, , 9]
rownames(wadj_late) <- var_names
colnames(wadj_late) <- var_names

g_late <- graph_from_adjacency_matrix(abs(wadj_late), mode = "undirected", weighted = TRUE, diag = FALSE)

# A. Community Detection (Walktrap or Louvain algorithm for weighted networks)
communities <- cluster_walktrap(g_late)
cat("\n--- Community Detection Structure (Late Era) ---\n")

--- Community Detection Structure (Late Era) ---

print(membership(communities))

 INJTIME  INJ2DTH DEATHTIM   NUMVIC   NUMOFF    V1SEX 
       7        8        7        9        1        4 
   V1AGE   V1RACE  PRIORV1 V1SXRACE     AREA    PLACE 
       1        4        1        4        1        1 
   PHOME PVEHICLE  PSTREET   CIRCUM MOTIVROB     GANG 
       1        1        1        2        2        2 
SYNDROME VICCRIME RELATION INTXUSED   LIQUOR DRGSUSED 
       2        1        1        3        3        6 
    DRUG INTOXTOT  DRUGTOT DRUGRELA   WEAPON   WHANDS 
       6        3        3        2        1        1 
WAUTOMAT     OSEX     OAGE    ORACE  PRIOROF  OSXRACE 
       1        5        1        5        1        5 
 CALIBER 
       1 

# B. Core-Periphery Decomposition (Finding the structural core of crimes)
coreness_values <- coreness(g_late)
cat("\n--- K-Core Decomposition (Structural Core) ---\n")

--- K-Core Decomposition (Structural Core) ---

print(coreness_values[coreness_values > 2])

  NUMOFF    V1SEX    V1AGE   V1RACE  PRIORV1 V1SXRACE 
       8        8        9        5        9        3 
    AREA    PLACE    PHOME  PSTREET   CIRCUM MOTIVROB 
       9        9        8        8        9        7 
    GANG SYNDROME VICCRIME RELATION INTXUSED   LIQUOR 
       4        9        9        9        9        7 
INTOXTOT  DRUGTOT DRUGRELA   WEAPON WAUTOMAT     OSEX 
       8        6        7        9        9        3 
    OAGE    ORACE  PRIOROF  OSXRACE  CALIBER 
       8        3        6        3        9 

# C. Global Network Statistics
cat("\n--- Global Network Metrics (Late Era) ---\n")

--- Global Network Metrics (Late Era) ---

cat("Global Clustering Coefficient (Transitivity):", transitivity(g_late, type = "global"), "\n")

Global Clustering Coefficient (Transitivity): 0.5660377 

cat("Average Path Length:", mean_distance(g_late, directed = FALSE), "\n")

Average Path Length: 0.8623305 

Output as tables

library(igraph)
library(dplyr)

# Extract variable names
var_names <- actual_cols_final
n_estpoints <- dim(fit_tvmgm$pairwise$wadj)[3]

# Initialize a master list to store metrics for all time points
network_metrics_list <- list()

for (tp in 1:n_estpoints) {
  # 1. Extract weighted adjacency matrix for estimation point `tp`
  wadj_matrix <- fit_tvmgm$pairwise$wadj[, , tp]
  rownames(wadj_matrix) <- var_names
  colnames(wadj_matrix) <- var_names
  
  # 2. Convert to an igraph object (treating weights as absolute values)
  g <- graph_from_adjacency_matrix(
    abs(wadj_matrix), 
    mode = "undirected", 
    weighted = TRUE, 
    diag = FALSE
  )
  
  # 3. Compute Centrality & Network Statistics
  node_stats <- data.frame(
    Estimation_Point = tp,
    Variable         = var_names,
    Degree           = degree(g),
    Strength         = strength(g),                    # Sum of absolute edge weights
    Betweenness      = betweenness(g, normalized = TRUE),
    Closeness        = closeness(g, normalized = TRUE)
  )
  
  network_metrics_list[[tp]] <- node_stats
}

# Combine all estimation points into one comprehensive master table
master_network_table <- bind_rows(network_metrics_list)

# -------------------------------------------------------------------------
# OUTPUT TO CSV: Node-Level Metrics Across All Estimation Points
# -------------------------------------------------------------------------
output_file_metrics <- "tvmgm_node_centrality_metrics.csv"
write.csv(master_network_table, output_file_metrics, row.names = FALSE)
cat("Successfully exported node metrics table to:", output_file_metrics, "\n")

# Select a specific time slice, e.g., Estimation Point 9 (Late Era)
target_tp <- 9
wadj_late <- fit_tvmgm$pairwise$wadj[, , target_tp]
rownames(wadj_late) <- var_names
colnames(wadj_late) <- var_names

g_late <- graph_from_adjacency_matrix(abs(wadj_late), mode = "undirected", weighted = TRUE, diag = FALSE)

# A. Community Detection (Walktrap algorithm)
communities <- cluster_walktrap(g_late)
community_df <- data.frame(
  Variable  = var_names,
  Community = membership(communities)
)

# B. Core-Periphery Decomposition (K-Core)
coreness_values <- coreness(g_late)
core_df <- data.frame(
  Variable = var_names,
  K_Core   = coreness_values
)

# Merge Community and Core dataframes together for the target time point
structural_summary_df <- community_df %>%
  left_join(core_df, by = "Variable") %>%
  left_join(
    master_network_table %>% filter(Estimation_Point == target_tp) %>% select(-Estimation_Point),
    by = "Variable"
  )

# -------------------------------------------------------------------------
# OUTPUT TO CSV: Structural Communities and Cores for Target Era
# -------------------------------------------------------------------------
output_file_structural <- paste0("tvmgm_structural_summary_estpoint_", target_tp, ".csv")
write.csv(structural_summary_df, output_file_structural, row.names = FALSE)
cat("Successfully exported structural summary table to:", output_file_structural, "\n")

# C. Global Network Statistics (Exported as a separate summary row/file)
global_stats_df <- data.frame(
  Estimation_Point            = target_tp,
  Global_Clustering_Coef      = transitivity(g_late, type = "global"),
  Average_Path_Length         = mean_distance(g_late, directed = FALSE),
  Edge_Density                = edge_density(g_late),
  Total_Edges                 = gsize(g_late)
)

output_file_global <- paste0("tvmgm_global_network_stats_estpoint_", target_tp, ".csv")
write.csv(global_stats_df, output_file_global, row.names = FALSE)
cat("Successfully exported global network statistics to:", output_file_global, "\n")