Theoretical model

---

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")
invisible(fpackage.check(packages))
rm(packages)

---

2 Simulation

Set up parameter space:

1. Scale-free networks using the configuration model (e.g., Newman et al., 2001).
2. Small-world networks using the Watts-Strogatz (1998) model.

# full factorial design
n <- 100
conf <- expand.grid(group_size = n, minority_prop = c(0.05, 0.1, 0.15), min_deg = c(2), max_deg = c(2 *
    (sqrt(n)), n - 1), dist = c("power-law", "log-normal"), alpha = c(2.1, 2.5, 3), r_kk = seq(-0.4,
    0.1, length.out = 6), rho_kx = seq(0, 0.5, length.out = 6), influence = c("strong", "weak"), choice_rule = c("deterministic",
    "probabilistic"))

# apply filters; we don't need to simulate the whole parameter space for weak influence...
conf <- conf %>%
    filter(!(influence == "weak" & minority_prop == 0.05)) %>%
    filter(!(influence == "weak" & minority_prop > 0.05 & (rho_kx < 0.2 | r_kk > -0.2)))

# parameter space 2: 'small-world' networks using WS-model
ws <- expand.grid(group_size = n, min_deg = c(3, 4), minority_prop = c(0.05, 0.1, 0.15), p = c(0.01,
    0.05, 0.1, 0.25, 1), r_kk = seq(-0.4, 0.1, length.out = 6), rho_kx = seq(0, 0.5, length.out = 6),
    influence = c("weak", "strong"), choice_rule = c("deterministic", "probabilistic"))

# apply filters; we don't need to simulate the whole parameter space for weak influence...
ws <- ws %>%
    filter(!(influence == "weak" & minority_prop == 0.05)) %>%
    filter(!(influence == "weak" & (rho_kx < 0.2 | r_kk > -0.2)))

# also, we don't manipulate degree-trait correlation and disassortativity when the network is
# nearly regular
ws <- ws %>%
    filter(!(p <= 0.05 & (rho_kx != 0 | r_kk != 0)))

fshowdf(fdesign(conf), caption = "Configuration model design space")

Configuration model design space

parameter	n_levels	levels
group_size	1	100
minority_prop	3	0.05, 0.1, 0.15
min_deg	1	2
max_deg	2	20, 99
dist	2	power-law, log-normal
alpha	3	2.1, 2.5, 3
r_kk	6	-0.4, -0.3, -0.2, -0.1, 0, 0.1
rho_kx	6	0, 0.1, 0.2, 0.3, 0.4, 0.5
influence	2	strong, weak
choice_rule	2	deterministic, probabilistic

fshowdf(fdesign(ws), caption = "Watts-Strogatz model design space")

Watts-Strogatz model design space

parameter	n_levels	levels
group_size	1	100
min_deg	2	3, 4
minority_prop	3	0.05, 0.1, 0.15
p	5	0.01, 0.05, 0.1, 0.25, 1
r_kk	6	-0.4, -0.3, -0.2, -0.1, 0, 0.1
rho_kx	6	0, 0.1, 0.2, 0.3, 0.4, 0.5
influence	2	weak, strong
choice_rule	2	deterministic, probabilistic

Simulate norm evolution across N seeds for all target networks:

2.1 Configuration model

# first a test:

# 1 create a network generate degree sequence
degseq <- fdegseq(n = 50, alpha = 2.5, k_min = 2, k_max = 25, dist = "power-law", seed = 123)

# construct network from degree sequence
network <- sample_degseq(degseq, method = "vl")

# assign roles to nodes
min_prop = 0.2
V(network)$role <- sample(c(rep("trendsetter", floor(50 * min_prop)), rep("conformist", 50 - floor(50 *
    min_prop))))

# simulate

fabm(network = network, params = list(s = 15, e = 10, w = 40, z = 50, lambda1 = 5, lambda2 = 1.8), max_rounds = 50,
    mi_threshold = 0.49, choice_rule = "deterministic")

# number of seeds
nIter = 3

# set up parallel backend to increase efficiency
ncores <- detectCores() - 1 
cl <- makeCluster(ncores)
registerDoParallel(cl)

# make folder to store simulations in
if (!dir.exists("./sims")) dir.create("./sims")

# parallel processing using foreach
system.time({
  foreach(i = 1:nrow(conf), .combine = 'c', .packages = c("igraph", "tidyverse")) %dopar% {
    
    cfg <- conf[i, ] # get configuration from full factorial
    results <- list()  # temporary storage for all iterations of this config
    
    for (iter in 1:nIter) {
      seed <- 123 + iter 
      set.seed(seed)
      
      results[[iter]] <- tryCatch({
        
        # generate degree sequence
        degseq <- fdegseq(n = cfg$group_size, 
                          alpha = cfg$alpha, 
                          k_min = cfg$min_deg, 
                          k_max = cfg$max_deg, 
                          dist = cfg$dist, 
                          seed = seed)
 
        # construct network
        network <- sample_degseq(degseq, method = "vl")
      
        # assign roles
        V(network)$role <- sample(
          c(
            rep("trendsetter", floor(cfg$group_size * cfg$minority_prop)),
            rep("conformist", cfg$group_size - floor(cfg$group_size * cfg$minority_prop))
          )
        )
        #fplot_graph(network)
      
        # rewire and swap
        rewired_network <- frewire_r(network, cfg$r_kk, verbose = FALSE, max_iter = 1e5)
        final_network <- fswap_rho(rewired_network, cfg$rho_kx, verbose = FALSE, max_iter = 1e4)
        actual_r <- assortativity_degree(rewired_network)
        final_rho <- fdegtraitcor(final_network)$cor
      
        # set initial action
        V(final_network)$action <- ifelse(V(final_network)$role == "trendsetter", 1, 0)
      
        # calculate global majority illusion
        # get threshold based on influence strength
        thresh <- ifelse(cfg$influence=="strong", .49, .50)
        mi <- fcalculate_majority_illusion(final_network, threshold = thresh)
      
        params <- if (cfg$influence == "strong") {
          list(s = 15, e = 10, w = 40, z = 50, lambda1 = 5, lambda2 = 1.8)
        } else {
          list(s = 15, e = 10, w = 40, z = 50, lambda1 = 3, lambda2 = 1.8)
        }
      
        # run simulation
        sim <- fabm(
          network = final_network,
          params = params,
          max_rounds = 50,
          mi_threshold = thresh,
          choice_rule = cfg$choice_rule
        )
        
        # package result
        list(
          config_id = i,
          dist = cfg$dist,
          alpha = cfg$alpha,
          target_r = cfg$r_kk,
          actual_r = actual_r,
          target_rho = cfg$rho_kx,
          actual_rho = final_rho,
          mi = mi,
          seed = seed,
          choice_rule = cfg$choice_rule,
          influence = cfg$influence,
          min_deg = cfg$min_deg,
          max_deg = cfg$max_deg,
          minority_prop = cfg$minority_prop,
          sim = list(
            outcomes = sim$outcomes,
            equilibrium = sim$equilibrium
          )
        )
        }, error = function(e) {
        message(sprintf("Iteration %d failed for config %d: %s", iter, i, e$message))
        NULL
      })
    }
    
    results <- Filter(Negate(is.null), results)
    
    # save all iterations for this configuration as a single file
    saveRDS(results, file = paste0("./sims/results_config_", i, ".rds"))
  }
})
#for 3168 confs and 50 iterations, about 18 hours.

# now load in the results
results_dir <- "./sims/"
files <- list.files(results_dir, pattern = "^results_config_\\d+\\.rds$", full.names = TRUE)

# to a long dataframe
all_results <- foreach(file = files, .combine = bind_rows, .packages = c("tibble", "dplyr")) %dopar% {
  config_results <- readRDS(file)
  summaries <- lapply(config_results, function(res) {
    eq <- res$sim$equilibrium
    seg <- res$sim$segregation
    tibble(
      config_id = res$config_id,
      alpha = res$alpha,
      dist = res$dist,
      target_r = res$target_r,
      actual_r = res$actual_r,
      target_rho = res$target_rho,
      actual_rho = res$actual_rho,
      mi = res$mi,
      influence = res$influence,
      choice_rule = res$choice_rule,
      minority_prop = res$minority_prop,
      min_deg = res$min_deg,
      max_deg = res$max_deg,
      seed = res$seed,
      reached_equilibrium = eq$reached,
      rounds_to_equilibrium = eq$round,
      final_adoption_rate = eq$prop_follow_trend,
      beh_assort = eq$segregation$r_behavior,
      L1CC_share = eq$segregation$L1CC_share
      
    )
  })
  
  bind_rows(summaries)
}

# stop cluster
stopCluster(cl)

# order by config_id
data <- all_results[order(all_results$config_id), ]

# and save the  dataframe
fsave(data, "sims_conf.Rda")

2.2 Watts-Strogatz model

#if (!dir.exists("./sims2")) dir.create("./sims2")

system.time({
  foreach(i = 1:10, .combine = 'c', .packages = c("igraph", "tidyverse")) %dopar% {
    
    cfg <- ws[i, ] # get configuration from full factorial
    results <- list()  # temporary storage for all iterations of this config
    
    for (iter in 1:nIter) {
      seed <- 123 + iter 
      set.seed(seed)
      
      results[[iter]] <- tryCatch({
        
        # construct network
        network <- sample_smallworld(dim = 1, size = cfg$group_size, nei = cfg$min_deg, p = cfg$p)
        
        # remove isolates that may arise due to rewiring
        isolates <- which(degree(network)==0)
        if (length(isolates) > 0) {
          network <- delete_vertices(network, isolates)
        }
        
        new_n <- vcount(network)
        
        # assign roles
        V(network)$role <- sample(
          c(
            rep("trendsetter", floor(new_n * cfg$minority_prop)),
            rep("conformist", new_n - floor(new_n * cfg$minority_prop))
          )
        )

        # rewire and swap
        rewired_network <- frewire_r(network, cfg$r_kk, verbose = TRUE, max_iter = 1e5)
        final_network <- fswap_rho(rewired_network, cfg$rho_kx, verbose = FALSE, max_iter = 1e4)
        actual_r <- assortativity_degree(rewired_network)
        final_rho <- fdegtraitcor(final_network)$cor
        
        # set initial action
        V(final_network)$action <- ifelse(V(final_network)$role == "trendsetter", 1, 0)
        
        #fplot_graph(final_network)
        
        # calculate global majority illusion
        # get threshold based on influence strength
        thresh <- ifelse(cfg$influence=="strong", .49, .50)
        mi <- fcalculate_majority_illusion(final_network, threshold = thresh)
        
        params <- if (cfg$influence == "strong") {
          list(s = 15, e = 10, w = 40, z = 50, lambda1 = 5, lambda2 = 1.8)
        } else {
          list(s = 15, e = 10, w = 40, z = 50, lambda1 = 3, lambda2 = 1.8)
        }
        
        # run simulation
        sim <- fabm(
          network = final_network,
          params = params,
          max_rounds = 50,
          mi_threshold = thresh,
          choice_rule = cfg$choice_rule
        )
        
        # package result
        list(
          config_id = i,
          min_deg = cfg$min_deg,
          p = cfg$p,
          target_r = cfg$r_kk,
          actual_r = actual_r,
          target_rho = cfg$rho_kx,
          actual_rho = final_rho,
          mi = mi,
          seed = seed,
          choice_rule = cfg$choice_rule,
          influence = cfg$influence,
          group_size = new_n,
          minority_prop = cfg$minority_prop,
          sim = list(
            outcomes = sim$outcomes,
            equilibrium = sim$equilibrium
          )
          
        )
      }, error = function(e) {
        message(sprintf("Iteration %d failed for config %d: %s", iter, i, e$message))
        NULL
      })
    }
    
    results <- Filter(Negate(is.null), results)
    
    # save all iterations for this configuration as a single file
    saveRDS(results, file = paste0("./sims2/results_config_", i, ".rds"))
  }
})

# now load in the results
results_dir <- "./sims2/"
files <- list.files(results_dir, pattern = "^results_config_\\d+\\.rds$", full.names = TRUE)

# to a long dataframe
all_results <- foreach(file = files, .combine = bind_rows, .packages = c("tibble", "dplyr")) %dopar% {
  config_results <- readRDS(file)
  summaries <- lapply(config_results, function(res) {
    eq <- res$sim$equilibrium
    
    tibble(
      config_id = res$config_id,
      size = res$group_size,
      min_deg = res$min_deg,
      p = res$p,
      target_r = res$target_r,
      actual_r = res$actual_r,
      target_rho = res$target_rho,
      actual_rho = res$actual_rho,
      mi = res$mi,
      minority_prop = res$minority_prop,
      influence = res$influence,
      choice_rule = res$choice_rule,
      seed = res$seed,
      reached_equilibrium = eq$reached,
      rounds_to_equilibrium = eq$round,
      final_adoption_rate = eq$prop_follow_trend,
      beh_assort = eq$segregation$r_behavior,
      L1CC_share = eq$segregation$L1CC_share
    )
  })
  
  bind_rows(summaries)
}

stopCluster(cl)

# order by config_id
data <- all_results[order(all_results$config_id), ]

#fix(data)
fsave(data, "sims_sw.Rda")

---

3 Results

We use a heatmap to visualize the probability with which each target network configuration drives negative equilibrium, across seeds.

3.1 Configuration model

# import data
today = "20250627"  # date on which the data was saved
data <- fload(paste0("./data/processed/", today, "sims_conf.Rda"))
# str(data)

# define equilibrium definition depends on deterministic/stochastic choice-rule:
data$unpop <- NA
data$unpop[data$choice_rule == "deterministic"] <- ifelse(data$final_adoption_rate[data$choice_rule ==
    "deterministic"] == 1, 1, 0)  # full adoption

# for stochastic choice-rule, use alternative definition: where proportion taking up B >= p +
# (1-p)(1-e)
data$unpop[data$choice_rule == "probabilistic"] <- ifelse(data$final_adoption_rate[data$choice_rule ==
    "probabilistic"] >= data$minority_prop[data$choice_rule == "probabilistic"] + (1 - data$minority_prop[data$choice_rule ==
    "probabilistic"]) * (1 - 0.1), 1, 0)

3.1.1 strong influence

3.1.1.1 truncated

# make plots: 1 determinstic (main)
fwrapper(data = data, choice_rule = "deterministic", influence = "strong", kmin = 2, kmax = 2 * sqrt(n),
    minority_props = c(0.05, 0.1, 0.15))

# ggsave('./figures/det_strong_trunc.png', height=4, width = 14)

# 2 stochastic
fwrapper(data = data, choice_rule = "probabilistic", influence = "strong", kmin = 2, kmax = 2 * sqrt(n),
    minority_props = c(0.05, 0.1, 0.15))

# ggsave('./figures/stoch_strong_trunc.png', height=4, width = 14)

3.1.1.2 untruncated

fwrapper(data = data, choice_rule = "deterministic", influence = "strong", kmin = 2, kmax = n - 1, minority_props = c(0.05,
    0.1, 0.15))

# ggsave('./figures/det_strong_untrunc.png', height=4, width = 14)

fwrapper(data = data, choice_rule = "probabilistic", influence = "strong", kmin = 2, kmax = n - 1, minority_props = c(0.05,
    0.1, 0.15))

# ggsave('./figures/stoch_strong_untrunc.png', height=4, width = 14)

3.1.2 weak influence

3.1.2.1 truncated

# make plots: 1 determinstic (main)
fwrapper(data = data, choice_rule = "deterministic", influence = "weak", kmin = 2, kmax = 2 * sqrt(n),
    minority_props = c(0.1, 0.15))

# 2 stochastic
fwrapper(data = data, choice_rule = "probabilistic", influence = "weak", kmin = 2, kmax = 2 * sqrt(n),
    minority_props = c(0.1, 0.15))

3.1.2.2 untruncated

fwrapper(data = data, choice_rule = "deterministic", influence = "weak", kmin = 2, kmax = n - 1, minority_props = c(0.1,
    0.15))

fwrapper(data = data, choice_rule = "probabilistic", influence = "strong", kmin = 2, kmax = n - 1, minority_props = c(0.1,
    0.15))

3.2 Watts-Strogatz model

# import data
today = "20250627"  # date on which the data was saved
data <- fload(paste0("./data/processed/", today, "sims_sw.Rda"))
# str(data)

# define equilibrium definition depends on deterministic/stochastic choice-rule:
data$unpop <- NA
data$unpop[data$choice_rule == "deterministic"] <- ifelse(data$final_adoption_rate[data$choice_rule ==
    "deterministic"] == 1, 1, 0)  # full adoption

# for stochastic choice-rule, use alternative definition: where proportion taking up B >= p +
# (1-p)(1-e)
data$unpop[data$choice_rule == "probabilistic"] <- ifelse(data$final_adoption_rate[data$choice_rule ==
    "probabilistic"] >= data$minority_prop[data$choice_rule == "probabilistic"] + (1 - data$minority_prop[data$choice_rule ==
    "probabilistic"]) * (1 - 0.1), 1, 0)

3.2.1 strong influence

# make plots: 1 determinstic (main)
fwrapper(data = data, fplot = fcreate_heatmap_sw, choice_rule = "deterministic", influence = "strong",
    minority_props = c(0.05, 0.1, 0.15))

# ggsave('./figures/det_strong_sw.png', height=4, width = 18)

# stochastic
fwrapper(data = data, fplot = fcreate_heatmap_sw, choice_rule = "probabilistic", influence = "strong",
    minority_props = c(0.05, 0.1, 0.15))

# =ggsave('./figures/stoch_strong_sw.png', height=4, width = 18)

3.2.2 weak influence

fwrapper(data = data, fplot = fcreate_heatmap_sw, choice_rule = "deterministic", influence = "weak",
    minority_props = c(0.1, 0.15))

# ggsave('./figures/det_weak_sw.png', height=4, width = 14)

fwrapper(data = data, fplot = fcreate_heatmap_sw, choice_rule = "probabilistic", influence = "weak",
    minority_props = c(0.1, 0.15))

# ggsave('./figures/stoch_weak_sw.png', height=4, width = 14)

explore clustering (auto-correlation) of norms in the network; that is, dyadic similarity (behavioral assortativity). but beyond that ‘pocketing’ (e.g., largest component size/share); or segregation at distance-1 and distance-k>1 or distance decay.

# all_results$beh_assort plot(all_results$beh_assort, all_results$p)

---

References

Newman, M. E. J., S. H. Strogatz, and D. J. Watts. 2001. “Random Graphs with Arbitrary Degree Distributions and Their Applications.” Physical Review E 64 (2): 026118. https://doi.org/10.1103/PhysRevE.64.026118.

Watts, Duncan J., and Steven H. Strogatz. 1998. “Collective Dynamics of ‘Small-World’ Networks.” Nature 393 (6684): 440–42. https://www.nature.com/articles/30918/figu-.

---

Copyright © Rob Franken