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", "lubridate", "ggtext")

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

---

2 Experimental conditions

2.1 1. heavy-tailed network with central minorities

Let’s take a likely-case for an unpopular norm, generate multiple networks of this target, and simulate the emergence of the norm following a deterministic and probabilistic choice model:

# pick one configuration that likely leads to an unpopular norm, and explore multiple 'seeds':
run_one_seed <- function(
  i,
  base_seed = 2253281,
  params = list(s = 15, e = 10, w = 40, z = 50, lambda1 = 5, lambda2 = 1.8),
  
  # tweak network
  k_min = 2, 
  k_max = 20,
  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      = 50,
    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", 5), rep("conformist", 45))
  )
  

  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
}

test <- run_one_seed(2, k_min = 2, k_max = 20, alpha = 2.4, rho = 0.4, r = -0.1, return_network = TRUE)

base = 2253281
seed = base + 2

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

# let's check the distribution of distances to 'trendsetters', among 'conformists'
g <- test$network
trend <- which(V(g)$role == "trendsetter")
conf <- which(V(g)$role == "conformist")

# distances: rows = sources (trendsetters), cols = all vertices
D <- distances(g, v = trend, to = V(g), mode = "all")
Dc <- D[, conf]  #keep just conformists

# calculate (a) the distance to the nearest trendsetter and (b) the average distance to
# trendsetters, across conformists
dists <- data.frame(shortest = apply(Dc, 2, min), average = apply(Dc, 2, function(x) mean(x[is.finite(x)])))
dists_long <- pivot_longer(dists, cols = c(shortest, average), names_to = "type", values_to = "distance_to_seed")

# plot the distribution
ggplot(dists_long, aes(x = distance_to_seed, fill = type)) + geom_histogram(bins = 20, alpha = 0.5, position = "identity") +
    theme_minimal()

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

test <- run_one_seed(3, k_min = 2, k_max = 20, alpha = 2.4, rho = 0.4, r = -0.1)

seed = seed + 1

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

test <- run_one_seed(4, k_min = 2, k_max = 20, alpha = 2.4, rho = 0.4, r = -0.1, return_network = TRUE)

base = 2253281
seed = base + 4

# let's check the distribution of distances to 'trendsetters', among 'conformists'
g <- test$network
trend <- which(V(g)$role == "trendsetter")
conf <- which(V(g)$role == "conformist")

# distances: rows = sources (trendsetters), cols = all vertices
D <- distances(g, v = trend, to = V(g), mode = "all")
Dc <- D[, conf]  #keep just conformists

# calculate (a) the distance to the nearest trendsetter and (b) the average distance to
# trendsetters, across conformists
dists <- data.frame(shortest = apply(Dc, 2, min), average = apply(Dc, 2, function(x) mean(x[is.finite(x)])))
dists_long <- pivot_longer(dists, cols = c(shortest, average), names_to = "type", values_to = "distance_to_seed")

# plot the distribution
ggplot(dists_long, aes(x = distance_to_seed, fill = type)) + geom_histogram(bins = 20, alpha = 0.5, position = "identity") +
    theme_minimal()

test$stats

#> $run
#> [1] 4
#> 
#> $seed
#> [1] 2253285
#> 
#> $num_nodes
#> [1] 50
#> 
#> $num_edges
#> [1] 89
#> 
#> $avg_degree
#> [1] 3.56
#> 
#> $sd_degree
#> [1] 3.091529
#> 
#> $net_density
#> [1] 0.07265306
#> 
#> $net_diameter
#> [1] 6
#> 
#> $avg_path_len
#> [1] 2.923265
#> 
#> $clust_coeff
#> [1] 0.1428571
#> 
#> $assort_deg
#> [1] -0.09244437
#> 
#> $deg_trait_cor
#> [1] 0.5489375
#> 
#> $components
#> [1] 1

table(degree(test$network))

#> 
#>  2  3  4  5  6  7  8  9 10 19 
#> 29  8  3  3  1  1  1  1  2  1

# identify trendsetters
ids <- which(V(test$network)$role == "trendsetter")

# check their centrality
sort(as.numeric(cbind(degree(test$network), V(test$network)$role)[ids]))

#> [1]  3  5  6 10 19

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

2.1.1 increased density

test <- run_one_seed(5, k_min = 2, k_max = 49, alpha = 2.4, rho = 0.4, r = -0.1, return_network = TRUE)

seed = seed + 1

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

# let's check the distribution of distances to 'trendsetters', among 'conformists'
g <- test$network
trend <- which(V(g)$role == "trendsetter")
conf <- which(V(g)$role == "conformist")

# distances: rows = sources (trendsetters), cols = all vertices
D <- distances(g, v = trend, to = V(g), mode = "all")
Dc <- D[, conf]  #keep just conformists

# calculate (a) the distance to the nearest trendsetter and (b) the average distance to
# trendsetters, across conformists
dists <- data.frame(shortest = apply(Dc, 2, min), average = apply(Dc, 2, function(x) mean(x[is.finite(x)])))
dists_long <- pivot_longer(dists, cols = c(shortest, average), names_to = "type", values_to = "distance_to_seed")

# plot the distribution
ggplot(dists_long, aes(x = distance_to_seed, fill = type)) + geom_histogram(bins = 20, alpha = 0.5, position = "identity") +
    theme_minimal()

Let’s take this as our experimental network configuration.

test$stats

#> $run
#> [1] 5
#> 
#> $seed
#> [1] 2253286
#> 
#> $num_nodes
#> [1] 50
#> 
#> $num_edges
#> [1] 113
#> 
#> $avg_degree
#> [1] 4.52
#> 
#> $sd_degree
#> [1] 4.15142
#> 
#> $net_density
#> [1] 0.0922449
#> 
#> $net_diameter
#> [1] 6
#> 
#> $avg_path_len
#> [1] 2.669388
#> 
#> $clust_coeff
#> [1] 0.2378049
#> 
#> $assort_deg
#> [1] -0.1081219
#> 
#> $deg_trait_cor
#> [1] 0.4120337
#> 
#> $components
#> [1] 1

table(degree(test$network))

#> 
#>  2  3  4  5  6  7 10 11 14 15 21 
#> 25  3  7  4  2  3  1  1  1  2  1

# identify trendsetters
ids <- which(V(test$network)$role == "trendsetter")

# check their centrality
sort(as.numeric(cbind(degree(test$network), V(test$network)$role)[ids]))

#> [1]  5  7 10 11 15

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

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

2.2 2. random network

with same number of ties

# ?erdos.renyi.game g <- erdos.renyi.game(n=50, p=0.05) ?sample_gnm()
set.seed(124124)
network <- sample_gnm(n = 50, m = 113)
V(network)$role <- sample(c(rep("trendsetter", 5), rep("conformist", 45)))

fplot_graph(network)

test <- fabm(network = network, max_rounds = 50, mi_threshold = 0.49, choice_rule = "probabilistic",
    plot = TRUE, histories = TRUE)

test$plot

# let's check the distribution of distances to 'trendsetters', among 'conformists'
g <- network
trend <- which(V(g)$role == "trendsetter")
conf <- which(V(g)$role == "conformist")

# distances: rows = sources (trendsetters), cols = all vertices
D <- distances(g, v = trend, to = V(g), mode = "all")
Dc <- D[, conf]  #keep just conformists

# calculate (a) the distance to the nearest trendsetter and (b) the average distance to
# trendsetters, across conformists
dists <- data.frame(shortest = apply(Dc, 2, min), average = apply(Dc, 2, function(x) mean(x[is.finite(x)])))
dists_long <- pivot_longer(dists, cols = c(shortest, average), names_to = "type", values_to = "distance_to_seed")

# plot the distribution
ggplot(dists_long, aes(x = distance_to_seed, fill = type)) + geom_histogram(bins = 20, alpha = 0.5, position = "identity") +
    theme_minimal()

# --- stats ---
stats <- list(num_nodes = vcount(g), num_edges = ecount(g), avg_degree = mean(degree(g)), sd_degree = sd(degree(g)),
    net_density = edge_density(g), net_diameter = diameter(g, directed = FALSE, unconnected = TRUE),
    avg_path_len = average.path.length(g, directed = FALSE), clust_coeff = transitivity(g, type = "global"),
    assort_deg = assortativity_degree(g), deg_trait_cor = fdegtraitcor(g)$cor, components = components(g)$no)
stats

#> $num_nodes
#> [1] 50
#> 
#> $num_edges
#> [1] 113
#> 
#> $avg_degree
#> [1] 4.52
#> 
#> $sd_degree
#> [1] 1.65665
#> 
#> $net_density
#> [1] 0.0922449
#> 
#> $net_diameter
#> [1] 5
#> 
#> $avg_path_len
#> [1] 2.730612
#> 
#> $clust_coeff
#> [1] 0.08387097
#> 
#> $assort_deg
#> [1] 0.1327246
#> 
#> $deg_trait_cor
#> [1] -0.02439024
#> 
#> $components
#> [1] 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(network))

# get roles
role_vector <- ifelse(V(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_n50_random.json")

---

3 import data

test <- read.csv("./rawdata/all_apps_wide_2026-02-02.csv")
times <- read.csv("./rawdata/PageTimes-2026-02-02.csv")

---

4 pilot results

On 2-2-2026, I recruited 80 Prolific participants, to populate a network of N=50 (with a 10% minority group):

4.1 arrival from Prolific

# subset experimental session
test <- test[test$session.code == "hwkqjhm0", ]
times <- times[times$session_code == "hwkqjhm0", ]

test <- test %>%
    transmute(participant_id = participant.code, participant_label = participant.label, id_in_session = participant.id_in_session,
        consent_given = consent.1.player.consent, consent_timestamp = consent.1.player.consent_timestamp,
        role = participant.role, is_dropout = participant.is_dropout, dropout_app = participant._current_app_name,
        comprehension_retries = comprehension.1.player.comprehension_retries, passed_comprehension = !participant._current_app_name %in%
            c("consent", "comprehension"), choice = unpop.1.player.choice)
test <- test[!is.na(test$consent_given), ]
test$bot <- ifelse(test$participant_label == "", 1, 0)

test <- test %>%
    mutate(consent_timestamp = ymd_hms(consent_timestamp), bot = factor(bot, levels = c(0, 1), labels = c("Prolific participant",
        "Bot"))) %>%
    arrange(consent_timestamp)

test <- test %>%
    mutate(arrival_order = row_number())

ggplot(test, aes(x = consent_timestamp, y = arrival_order)) + geom_point(aes(color = role, shape = bot),
    size = 3, alpha = 0.5) + scale_shape_manual(values = c(16, 2)) + scale_color_manual(values = c("blue",
    "red")) + labs(x = "Arrival time", y = "Arrival order", color = "Role", shape = "Type", title = "Participant arrivals from Prolific over time") +
    theme_minimal()

4.2 progression through the experiment

times <- times %>%
    mutate(timestamp = as_datetime(epoch_time_completed))

arrival_times <- times %>%
    group_by(participant_id_in_session) %>%
    summarize(arrival_time = min(timestamp), .groups = "drop")

times <- times %>%
    left_join(arrival_times, by = "participant_id_in_session") %>%
    mutate(participant_ordered = factor(participant_id_in_session, levels = arrival_times %>%
        arrange(arrival_time) %>%
        pull(participant_id_in_session)))

times_roles <- times %>%
    left_join(test %>%
        select(id_in_session, role), by = c(participant_id_in_session = "id_in_session"))

page_levels <- unique(times$page_name)

times_roles <- times_roles %>%
    mutate(page_name = factor(page_name, levels = page_levels))

custom_colors <- c(InitializeParticipant = "#c6dbef", ConsentPage = "#9ecae1", IntroductionPage = "#6baed6",
    ComprehensionPage = "#3182bd", NetworkFormationWaitPage = "#ffcc99", DecisionPage = "#ff9966", ResultsWaitPage = "#ff6666",
    ResultsPage = "#cc0033", FinalGameResults = "#660000")

# colored y-axis labels based on role
y_labels_colored <- times_roles %>%
    select(participant_ordered, role) %>%
    distinct() %>%
    arrange(participant_ordered) %>%
    mutate(label_colored = case_when(role == "Red" ~ paste0("<span style='color:red'>", participant_ordered,
        "</span>"), role == "Blue" ~ paste0("<span style='color:blue'>", participant_ordered, "</span>"),
        TRUE ~ paste0("<span style='color:darkgrey'>", participant_ordered, "</span>")))

# Create a named vector for scale_y_discrete labels
y_labels_vector <- y_labels_colored$label_colored
names(y_labels_vector) <- y_labels_colored$participant_ordered

# times2 <- times[times$participant_id_in_session %in% c(34,80),]

ggplot(times_roles, aes(x = timestamp, y = participant_ordered, color = page_name)) + geom_line(aes(group = participant_id_in_session),
    size = 1) + geom_point(size = 2) + scale_color_manual(values = custom_colors) + scale_y_discrete(labels = y_labels_vector) +
    labs(x = "Time", y = "Participant (ordered by arrival)", color = "Stage/Page", title = "Participant progression through experiment stages (by arrival)") +
    theme_minimal() + theme(axis.text.y = element_markdown(size = 6))

# by role
times_roles <- times %>%
    left_join(test %>%
        select(id_in_session, role), by = c(participant_id_in_session = "id_in_session"))

times_roles <- times_roles %>%
    mutate(page_name = factor(page_name, levels = page_levels))

colored_labels <- sapply(page_levels, function(stage) {
    paste0("<span style='color:", custom_colors[stage], "'>", stage, "</span>")
}, USE.NAMES = FALSE)

times_roles <- times_roles[!is.na(times_roles$role), ]

# Gantt-style plot faceted by participant role
ggplot(times_roles, aes(x = timestamp, y = page_name, group = participant_ordered, color = page_name)) +
    geom_step(direction = "hv", size = 0.5) + geom_point(size = 2) + scale_color_manual(values = custom_colors) +
    scale_y_discrete(labels = colored_labels) + labs(x = "Time", y = "Experiment Stage", color = "Stage/Page",
    title = "Participant progression through experiment stages by role") + facet_wrap(~role, scales = "free_y",
    ncol = 1) + theme_minimal() + theme(axis.text.y = element_markdown(size = 10), axis.text.x = element_text(angle = 45,
    hjust = 1), legend.position = "none")

4.3 comprehension

Retries

table(test$comprehension_retries, test$role)

#>     
#>      Blue Red
#>   0     5  38
#>   1    10  16
#>   2     3   8
#>   3     1   4
#>   4     1   3
#>   5     0   2
#>   7     0   1
#>   11    0   1
#>   12    0   1
#>   17    0   2
#>   21    0   1

# time per comprehension task
times2 <- times %>%
    filter(page_name == "ComprehensionPage") %>%
    mutate(time_spent = as.numeric(timestamp - arrival_time))

# exclude bots
bots <- test$id_in_session[test$bot == "Bot"]
times3 <- times2[!times2$participant_id_in_session %in% bots, ]

comprehension_data <- times3 %>%
    left_join(test %>%
        select(id_in_session, comprehension_retries, role), by = c(participant_id_in_session = "id_in_session"))

ggplot(comprehension_data, aes(x = time_spent, y = factor(comprehension_retries), color = role)) + geom_jitter(width = 0,
    height = 0.2, size = 2, alpha = 0.7) + scale_color_manual(values = c("blue", "red")) + labs(x = "Time spent on comprehension (seconds)",
    y = "Number of retries", color = "Role", title = "Time in seconds") + theme_minimal(base_size = 14)

4.4 choice behavior

# exclude bots and missing choices
test2 <- test %>%
    filter(bot != "Bot") %>%
    select(id_in_session, role, choice, comprehension_retries) %>%
    filter(!is.na(choice))

# scatter plot: choice vs role, size = retries
ggplot(test2, aes(x = role, y = as.numeric(choice), size = comprehension_retries)) + geom_jitter(width = 0.2,
    height = 0.05, alpha = 0.7, color = "steelblue") + scale_size_continuous(range = c(2, 8)) + labs(x = "Participant role",
    y = "Choice (1=blue, 0=red)", size = "Number of retries", title = "Participant choice by role and comprehension question retries") +
    theme_minimal(base_size = 14)

---

Copyright © Rob Franken