# load libraries
library(tidyverse)
library(tidybayes)
library(bayesplot)
library(ggplot2)
library(brms)

1. Data Preparataion

First, we load the preprocessed data frame d, add a variable NetworkAge, and rename variable Condition as NetworkType. The categorical variable NetworkType is effect-coded such that the contrasts of Full (or Rich) and Star (or Sparse) network conditions are set to 0.5 and -0.5, respectively. The data frame d contains, as the measure of disagreement on gestural convention, the entropy averaged across quad members for each unique combination of Quad, Item, and Round.

# load data
# data frame d will be loaded.
# see: preprocess_data_entropy.R
load('data-entropy.RData')

# add NetworkAge
d <-
  d %>%
  mutate(NetworkAge = ifelse(Condition == "F", Round * 2, Round)) %>%
  rename(NetworkType = Condition)

# Use a non-default contrast (effect-coding)
# Full (Rich): 0.5, Sparse: -0.5
contrasts(d$NetworkType) <- c(0.5, -0.5)

We define the community-level amount of disagreement as the entropy values averaged across Item and call it as divergence (Divergence1). (c.f., We also analyzed an alternative divergence measure Divergence2, which was defined as the median entropy values per Round per Community/Quad.) This is because (1) we are interested in the convergence of convention at the language level (a whole set of utterances across items), not at the individual word level; (2) convergences on the same utterances for different items may not be independent from each other, so treating Item as a random factor may not be ideal. In the below, we compute Divergence1 (or Divergence2) and log-transform (with base 2) it to create an outcome variable of our interest logDivergence1 (or logDivergence2). The result patterns are similar so we focus on the measure logDivergence1.

d.quad <-
  d %>%
  group_by(NetworkType, Quad, Round, NetworkAge) %>%
  summarise(Divergence1 = mean(Entropy),
            Nfeatures1 = mean(Nfeatures),
            Divergence2 = median(Entropy),
            Nfeatures2 = median(Nfeatures)) %>%
  ungroup() %>%
  mutate(logDivergence1 = log2(Divergence1),
         logDivergence2 = log2(Divergence2))

d.quad %>% head()

Column Nfeatures represents the average number of gestures (features) used by a quad as a whole per item. It can be interpreted as the average utterance length.

Figures in the below suggest that logDivergence1 decreased approximately linearly with Round. In other words, Divergence1 decreases multiplicatively.

d.quad %>%
  ggplot(aes(x = Round, y = logDivergence1, group = Quad, color = NetworkType)) + 
  geom_line() + 
  scale_x_continuous(breaks = 1:4, labels = 1:4) + 
  labs(x = "Round", 
       y = "log Divergence 1") + 
  scale_color_discrete(labels = c('Rich','Sparse')) +
  theme_bw(12) + 
  theme(legend.position = c(0.85, 0.8),
        panel.grid.minor = element_blank())

Given these trends, it seems reasonable to model logDivergence1 as a linear function of Round. Because we are interested in whether the convergence rate is modulated by NetworkType, we fitted data to a Bayesian linear mixed-effects model with Round and Round-by-NetworkType as fixed-effect terms and by-Quad random intercepts and random slopes of Round with their random correlations. Because the contrasts of two levels of NetworkType were set to 0.5 and -0.5, the coefficient of the interaction term will be the estimate of slope difference between two network types. Note that we did not include the main effect of NetworkType because, under random assignments of participants to different communities of different network types, the expected value of logDivergence1 must be same at time point of 0 (not observed) when no interaction occurred between community members.

2. Model

We assume that logDivergence of community (or Quad) \(c\) from NetworkType of \(n\) at Round \(r\) is normally distributed with mean \(\mu_{n,c,r}\) and variance \(\sigma^2\).

Round-specific within-community log-divergence \(\mu_{n,c,r}\) is modeled as a linear function of Round \(R\) and NetworkType \(NT\) as follows:

where \(\beta_0\) is the intercept (i.e., log divergence when \(R_r = 0\) [not observed]), \(\beta_1\) is the slope of Round, and \(\beta_2\) is the slope difference between two network types: \({\beta_1}_{F} - {\beta_1}_{S}\). \(b_{0,c}\) is the by-Quad random intercept and \(b_{1,c}\) is the by-Quad random slope of Round. Because NetworkType is a between-Quad factor, we do not include the random slope of the interaction term.

The priors for those paramters are presented below.

The prior on \(\beta_0\) expresses the prior belief that divergence of community-specific language at Round 0 (before interaction) is (soft-)bounded within the interval \([0.002, 512] (= [2^{-9}, 2^{9}])\) bits. Recall that we defined divergence as entropy values averaged across items. Divergence of 512 bits is expected in the situation where each community member chooses each of 512 gestural components randomly (with probability of 0.5), which is unrealistic.

The prior on \(\beta_1\) expresses the prior belief that the rate of change in divergence with respect to Round is soft-bounded within the interval \([1/8, 8] (= [2^{-3}, 2^{3}])\). This means that the prior expresses the belief that divergence at Round 4 could reduce up to \(0.00024 (= 1/ 2^{3 \cdot 4})\) times the initial divergence or increase up to \(4096\) times the initial divergence. While this prior belief restricts the effect sizes of Round, it still allows highly unlikely values.

The prior on \(\beta_2\) expresses the prior belief that the ratio of the rate of divergence between Full and Sparse Network communities is soft-bounded within the interval \([1/8, 8] (= [2^{-3}, 2^{3}])\).

2.1. Investigation of Prior

prior1 <- prior(normal(0, 3), class = Intercept) +
  prior(normal(0, 1), class = b, coef = "Round") +
  prior(normal(0, 1), class = b, coef = "Round:NetworkType1") +
  prior(student_t(3, 0, 5), class = sd) +
  prior(student_t(3, 0, 5), class = sigma) + 
  prior(lkj(2), class = cor)
# d.quad.logdiv1.bm1.prior <- brm(
#   data = d.quad,
#   logDivergence1 ~ 1 + Round + NetworkType:Round +
#     (Round | Quad),
#   control = list(adapt_delta = 0.95),
#   sample_prior = "only",       # do not use likelihood
#   prior = prior1,
#   iter = 5000, warmup = 2000,
#   chains = 4, cores = 4,
#   seed = 1000)
# 
# saveRDS(d.quad.logdiv1.bm1.prior, file = "d-quad-logdiv1-bm1-prior.Rds")
d.quad.logdiv1.bm1.prior <- readRDS("d-quad-logdiv1-bm1-prior.Rds")

The figure in the below presents possible patterns (expected log divergence values under the current priors on model parameters), suggesting that the prior covers a wide range of patterns while excluding unrealistic log divergence values.

p <- plot(marginal_effects(d.quad.logdiv1.bm1.prior, 
                           effects = "Round:NetworkType",
                           method = "fitted", 
                           spaghetti=TRUE, nsamples=100), plot = FALSE)

p[[1]] + 
  theme_bw(15)