Experiment

---

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 = 1253281123,
  params = list(s = 15, e = 10, w = 40, z = 50, lambda1 = 5, lambda2 = 1.8),
  
  # tweak network
  k_min = 2, 
  k_max = 9,
  alpha = 2.1,
  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      = 10,
    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", 2), rep("conformist", 8))
  )
  

  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   = 50,
    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             = 100,
    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
}

####

run_one_sw <- function(
  i,
  base_seed = 1253281123,
  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 = 10, nei = 3, p = beta) }
  else if (model == "caveman") {
    network <- simulate_caveman(n = 10, clique_size = clique_size) 
  }

  V(network)$role <- sample(
    c(rep("trendsetter", 2), rep("conformist", 8))
  )
  
  
  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   = 50,
    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             = 100,
    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
}

test <- run_one_seed(1, k_min = 2, k_max = 9, alpha = 2.1, rho = 0.7, r = -0.1, return_network = TRUE)

base = 1253281123
seed = base + 1

# 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, "network_test_n10.json")

# random
test <- run_one_sw(i = 15, beta = 1, pmin = 0.15, return_network = TRUE)

# fplot_graph(test$network)

# 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, "network_test_n10_random.json")

---

Copyright © Rob Franken