Experiment in SRDA computer-lab

---

1 Getting started

To copy the code, click the button in the upper right corner of the code-chunks.

1.1 clean up

rm(list = ls())
gc()

1.2 custom functions

We defined a number custom functions, at Download custom_functions.R.

source("./custom_functions.R")

1.3 necessary packages

• tidyverse: data wrangling
• igraph: generate and visualize graphs
• parallel: parallel computing to speed up simulation
• foreach: looping in parallel
• doParallel: parallel backend for foreach
• ggplot2: data visualization
• ggh4x: hacks for ggplot2
• ggpubr: make visualizations publication-ready

packages = c("tidyverse", "igraph", "ggplot2", "parallel", "doParallel", "foreach", "ggh4x", "ggpubr",
    "plotly", "RColorBrewer", "grid", "gridExtra", "patchwork", "ggplotify", "ggraph", "gganimate", "RColorBrewer",
    "ggtext", "magick", "jsonlite")

invisible(fpackage.check(packages))
rm(packages)

---

2 Experimental conditions

# pick one configuration that likely leads to an unpopular norm, and explore multiple 'seeds':
run_one_seed <- function(
  i,
  base_seed = 12532812,
  params = list(s = 15, e = 10, w = 40, z = 50, lambda1 = 5, lambda2 = 1.8),
  
  # tweak network
  k_min = 3, 
  k_max = 15,
  alpha = 2.4,
  rho = 0.4,
  r = -0.1,
  
  # retrieve network
  return_network = FALSE
  
) {
  # derived seed for this run
  seed_i <- base_seed + i
  set.seed(seed_i)

  # --- network creation ---
  degseq <- fdegseq(
    n = 20,
    alpha  = alpha,
    k_min  = k_min,
    k_max  = k_max,
    dist   = "log-normal", #use log-normal 
    seed   = seed_i
  )

  network <- sample_degseq(degseq, method = "vl")
  
  
   V(network)$role <- sample(
    c(rep("trendsetter", 3), rep("conformist", 17))
  )
  
  rewired_network <- frewire_r(network, r, verbose = FALSE, max_iter = 1e5)
  final_network   <- fswap_rho(rewired_network, rho, verbose = FALSE, max_iter = 1e4)
  
  # --- stats ---
  stats <- list(
    run           = i,
    seed          = seed_i,
    num_nodes     = vcount(final_network),
    num_edges     = ecount(final_network),
    avg_degree    = mean(degree(final_network)),
    sd_degree     = sd(degree(final_network)),
    net_density   = edge_density(final_network),
    net_diameter  = diameter(final_network, directed = FALSE, unconnected = TRUE),
    avg_path_len  = average.path.length(final_network, directed = FALSE),
    clust_coeff   = transitivity(final_network, type = "global"),
    assort_deg    = assortativity_degree(final_network),
    deg_trait_cor = fdegtraitcor(final_network)$cor,
    components    = components(final_network)$no
  )
  
  fplot_graph(final_network, layout = layout_with_fr(final_network)) 
  
  
  # --- initial actions ---
  V(final_network)$action <- ifelse(V(final_network)$role == "trendsetter", 1, 0)
  
  # --- deterministic simulation ---
  sim_det <- fabm(
    network      = final_network,
    params       = params,
    max_rounds   = 35,
    mi_threshold = 0.49,
    choice_rule  = "deterministic",
    plot         = TRUE,
    histories    = TRUE
  )
  
  # generate the gif for the current network
  gif_filename <- paste0("./figures/animation_network_", seed_i, ".gif")
  gif_path <- fnetworkgif(final_network, sim_det$decision_history, rounds = sim_det$equilibrium$round, output_dir = "./figures")
  # rename the gif to match the naming pattern
  file.rename(gif_path, gif_filename)

  if (!is.null(sim_det$plot)) {
    print(sim_det$plot)
  }
  
  # --- probabilistic simulation ---
  sim_prob <- fabm(
    network                = final_network,
    params                 = params,
    max_rounds             = 50,
    mi_threshold           = 0.49,
    choice_rule            = "probabilistic",
    stable_window          = 8,   # the length of the window of adoption values
    required_stable_rounds = 20, # number of windows needed to declare equilibrium
    plot                   = TRUE
  )
  if (!is.null(sim_prob$plot)) {
    print(sim_prob$plot)
  }
  
   result <- list(
    segregation_det  = sim_det$equilibrium$segregation,
    segregation_prob = sim_prob$equilibrium$segregation,
    stats            = stats
  )
  
  if (return_network) {
    result$network <- final_network
  }
  
  result
}

#and one "random" network
run_one_sw <- function(
  i,
  base_seed = 12532812,
  model = "watts-strogatz",
  beta = 0,
  nei = 3,
  clique_size = 5,
  pmin = 0.1,
  
  params = list(s = 15, e = 10, w = 40, z = 50, lambda1 = 5, lambda2 = 1.8),
  
  # retrieve network
  return_network = FALSE
  
) {
  # derived seed for this run
  seed_i <- base_seed + i
  set.seed(seed_i)

  if (model == "watts-strogatz") {
    network <- sample_smallworld(dim = 1, size = 20, nei = 3, p = beta) }
  else if (model == "caveman") {
    network <- simulate_caveman(n = 20, clique_size = clique_size) 
  }

  V(network)$role <- sample(
    c(rep("trendsetter", 3), rep("conformist", 17))
  )
  
  
  final_network <- network

  
  # --- stats ---
  stats <- list(
    run           = i,
    seed          = seed_i,
    num_nodes     = vcount(final_network),
    num_edges     = ecount(final_network),
    avg_degree    = mean(degree(final_network)),
    sd_degree     = sd(degree(final_network)),
    net_density   = edge_density(final_network),
    net_diameter  = diameter(final_network, directed = FALSE, unconnected = TRUE),
    avg_path_len  = average.path.length(final_network, directed = FALSE),
    clust_coeff   = transitivity(final_network, type = "global"),
    assort_deg    = assortativity_degree(final_network),
    deg_trait_cor = fdegtraitcor(final_network)$cor,
    components    = components(final_network)$no
  )
  
  stats_df <- data.frame(
    Metric = names(stats),
    Value  = unlist(stats),
    row.names = NULL
  )
  
  #print(stats_df)
  fplot_graph(final_network, layout = layout.kamada.kawai(final_network))
  

  # --- initial actions ---
  V(final_network)$action <- ifelse(V(final_network)$role == "trendsetter", 1, 0)
  
  # --- deterministic simulation ---
  sim_det <- fabm(
    network      = final_network,
    params       = params,
    max_rounds   = 35,
    mi_threshold = 0.49,
    choice_rule  = "deterministic",
    plot         = TRUE,
    histories    = TRUE
  )

  # generate the gif for the current network
  gif_filename <- paste0("./figures/animation_network_", seed_i, ".gif")
  gif_path <- fnetworkgif(final_network, sim_det$decision_history, rounds = sim_det$equilibrium$round, output_dir = "./figures")
  # rename the gif to match the naming pattern
  file.rename(gif_path, gif_filename)
  
  
  if (!is.null(sim_det$plot)) {
    print(sim_det$plot)
  }
  
  # --- probabilistic simulation ---
  sim_prob <- fabm(
    network                = final_network,
    params                 = params,
    max_rounds             = 50,
    mi_threshold           = 0.49,
    choice_rule            = "probabilistic",
    stable_window          = 8,   # the length of the window of adoption values
    required_stable_rounds = 20, # number of windows needed to declare equilibrium
    plot                   = TRUE
  )
  if (!is.null(sim_prob$plot)) {
    print(sim_prob$plot)
  }
  
  
  result <- list(
    segregation_det  = sim_det$equilibrium$segregation,
    segregation_prob = sim_prob$equilibrium$segregation,
    stats            = stats
  )
  
  if (return_network) {
    result$network <- final_network
  }
  
  result
}

---

2.1 SRDA computer-lab demonstation

2.1.1 WG1: heterogenous degree distribution; centralized “fanatics”:

test <- run_one_seed(38, k_min = 4, k_max = 16, alpha = 2.3, rho = 1, r = -0.2, return_network = TRUE)

table(degree(test$network))

#> 
#>  4  5  6  8  9 16 
#> 12  1  1  1  3  2

cbind(degree(test$network), V(test$network)$role)

#>       [,1] [,2]         
#>  [1,] "4"  "conformist" 
#>  [2,] "16" "trendsetter"
#>  [3,] "4"  "conformist" 
#>  [4,] "4"  "conformist" 
#>  [5,] "4"  "conformist" 
#>  [6,] "4"  "conformist" 
#>  [7,] "4"  "conformist" 
#>  [8,] "5"  "conformist" 
#>  [9,] "16" "trendsetter"
#> [10,] "9"  "conformist" 
#> [11,] "9"  "trendsetter"
#> [12,] "8"  "conformist" 
#> [13,] "4"  "conformist" 
#> [14,] "4"  "conformist" 
#> [15,] "4"  "conformist" 
#> [16,] "6"  "conformist" 
#> [17,] "4"  "conformist" 
#> [18,] "9"  "conformist" 
#> [19,] "4"  "conformist" 
#> [20,] "4"  "conformist"

base = 12532812
seed = base + 38
knitr::include_graphics(paste0("./figures/animation_network_", seed, ".gif"))

# use this as the network structure for an otree session:
# cbind(degree(test$network),V(test$network)$role)

# convert to adjacency matrix
adj_matrix <- as.matrix(as_adjacency_matrix(test$network))

# get roles
role_vector <- ifelse(V(test$network)$role == "trendsetter", 1, 0)
# create a list to store the network data
net <- list(adj_matrix = adj_matrix, role_vector = role_vector)
# save the list as a JSON file
write_json(net, "SRDA_n20_heterogenous.json")

2.1.2 WG2: homogenous degree distribution–“random network”

test <- run_one_sw(i = 513, beta = 1, pmin = 0.1, return_network = TRUE)

table(degree(test$network))

#> 
#> 2 3 4 5 6 7 8 9 
#> 1 1 2 4 2 6 3 1

cbind(degree(test$network), V(test$network)$role)

#>       [,1] [,2]         
#>  [1,] "4"  "conformist" 
#>  [2,] "8"  "conformist" 
#>  [3,] "7"  "conformist" 
#>  [4,] "7"  "conformist" 
#>  [5,] "8"  "conformist" 
#>  [6,] "7"  "conformist" 
#>  [7,] "7"  "conformist" 
#>  [8,] "3"  "trendsetter"
#>  [9,] "4"  "conformist" 
#> [10,] "5"  "conformist" 
#> [11,] "8"  "conformist" 
#> [12,] "6"  "conformist" 
#> [13,] "6"  "conformist" 
#> [14,] "5"  "conformist" 
#> [15,] "2"  "trendsetter"
#> [16,] "5"  "conformist" 
#> [17,] "7"  "conformist" 
#> [18,] "5"  "conformist" 
#> [19,] "9"  "conformist" 
#> [20,] "7"  "trendsetter"

base_seed = 12532812
seed = base_seed + 513

knitr::include_graphics(paste0("./figures/animation_network_", seed, ".gif"))

# use this as the network structure for an otree session:
# cbind(degree(test$network),V(test$network)$role)

# convert to adjacency matrix
adj_matrix <- as.matrix(as_adjacency_matrix(test$network))

# get roles
role_vector <- ifelse(V(test$network)$role == "trendsetter", 1, 0)
# create a list to store the network data
net <- list(adj_matrix = adj_matrix, role_vector = role_vector)
# save the list as a JSON file
write_json(net, "SRDA_n20_random.json")

---

Copyright © Rob Franken