# 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)

3. Bayesian Analysis of logDivergence1

# d.quad.logdiv1.bm1 <- brm(
#   data = d.quad,
#   logDivergence1 ~ 1 + Round + NetworkType:Round +
#     (Round | Quad),
#   control = list(adapt_delta = 0.99),
#   sample_prior = TRUE,
#   save_all_pars = TRUE,
#   prior = prior1,
#   iter = 5000, warmup = 2000,
#   chains = 4, cores = 4,
#   seed = 1000)
# 
# saveRDS(d.quad.logdiv1.bm1, file = "d-quad-logdiv1-bm1.Rds")
d.quad.logdiv1.bm1 <- readRDS("d-quad-logdiv1-bm1.Rds")
plot(d.quad.logdiv1.bm1, ask = FALSE, newpage = FALSE, N = 4)

The trace plots on the right column suggest the convergence of sampling distributions.

pp_check(d.quad.logdiv1.bm1, nsamples=100)

The thick black curve presents the distribution of empirical log divergence values (aggregated across all Quad, NetworkType, and Round values). The thin blue lines present the distributions of simulated data under the posterior joint distribution over model parameters.

print(summary(d.quad.logdiv1.bm1), digits = 3)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: logDivergence1 ~ 1 + Round + NetworkType:Round + (Round | Quad) 
##    Data: d.quad (Number of observations: 64) 
## Samples: 4 chains, each with iter = 5000; warmup = 2000; thin = 1;
##          total post-warmup samples = 12000
## 
## Group-Level Effects: 
## ~Quad (Number of levels: 16) 
##                      Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sd(Intercept)           0.264     0.076    0.132    0.434       5085 1.000
## sd(Round)               0.102     0.030    0.052    0.171       3582 1.001
## cor(Intercept,Round)   -0.103     0.317   -0.633    0.592       3822 1.001
## 
## Population-Level Effects: 
##                    Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## Intercept             2.445     0.080    2.284    2.606       6010 1.000
## Round                -0.411     0.031   -0.472   -0.350       5462 1.000
## Round:NetworkType1   -0.140     0.056   -0.250   -0.027       6045 1.001
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sigma    0.134     0.018    0.103    0.175       3686 1.002
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).
lab <- c('Slope of Round',
         'Slope difference\nbtw Network Structures')
temp <- c('b_Round',
          'b_Round:NetworkType1')

(p <- stanplot(d.quad.logdiv1.bm1, par = "^b_R",
           point_est = "median",
           prob = 0.5, 
           prob_outer = 0.95) +
  geom_vline(xintercept = 0, linetype = 3)  + 
  scale_y_discrete(labels=rev(lab), limits=rev(temp)) +
  xlim(-0.6, 0.005) +
  theme_bw(15))
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

p <- plot(marginal_effects(d.quad.logdiv1.bm1,
                           effects = "Round:NetworkType"), plot = FALSE)
p[[1]] + 
  geom_line(data = d.quad, 
            aes(x = as.numeric(Round), y = logDivergence1, group = Quad, 
                color = NetworkType,
                ymin = NULL, ymax = NULL, fill = NULL)) + 
  scale_color_discrete(labels = c('Rich','Sparse')) +
  scale_fill_discrete(labels = c('Rich', 'Sparse')) + 
  labs(y = "log divergence") +
  theme_bw(15)

The shaded regions present 95% credible intervals on conditional means of log divergence estimated from the present model. The spaghetti lines repsent observed data.

3.2. Bayes Factor

# prior1.base <- prior(normal(0, 3), class = Intercept) +
#   prior(normal(0, 1), class = b, coef = "Round") +
#   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.base <- brm(
#   data = d.quad,
#   logDivergence1 ~ 1 + Round + NetworkType +
#     (Round | Quad),
#   control = list(adapt_delta = 0.95),
#   save_all_pars = TRUE,
#   prior = prior1.base,
#   iter = 5000, warmup = 2000,
#   chains = 4, cores = 4,
#   seed = 1000)
# 
# saveRDS(d.quad.logdiv1.bm1.base, file = "d-quad-logdiv1-bm1-base.Rds")
d.quad.logdiv1.bm1.base <- readRDS("d-quad-logdiv1-bm1-base.Rds")
# set.seed(1005)
# BF10 <- c()
# for (i in 1:10) {
#   BF10 <- c(BF10, bayes_factor(d.quad.logdiv1.bm1, d.quad.logdiv1.bm1.base)$bf)
# }
# saveRDS(BF10, file = "d-quad-logdiv1-bm1-bf.Rds")
d.quad.logdiv1.bm1.BF10 <- readRDS(file = "d-quad-logdiv1-bm1-bf.Rds")
summary(d.quad.logdiv1.bm1.BF10)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.523   2.601   2.680   2.730   2.866   3.013

4. Frequentist Analysis of logDivergence1

library(lme4)
library(lmerTest)
library(optimx)
d.quad.logdiv1.lmer1 <- lmer(
  logDivergence1 ~ 1 + Round + NetworkType:Round +
    (Round | Quad), 
  REML = TRUE,
  data = d.quad,
  control = lmerControl(
    optimizer ='optimx', optCtrl=list(method='L-BFGS-B'))
  )
print(summary(d.quad.logdiv1.lmer1, ddf = "Kenward-Roger"), digits = 3)
## Linear mixed model fit by REML. t-tests use Kenward-Roger's method [
## lmerModLmerTest]
## Formula: logDivergence1 ~ 1 + Round + NetworkType:Round + (Round | Quad)
##    Data: d.quad
## Control: 
## lmerControl(optimizer = "optimx", optCtrl = list(method = "L-BFGS-B"))
## 
## REML criterion at convergence: -6.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1097 -0.6169  0.0055  0.5553  1.5917 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  Quad     (Intercept) 0.06253  0.250         
##           Round       0.00921  0.096    -0.27
##  Residual             0.01548  0.124         
## Number of obs: 64, groups:  Quad, 16
## 
## Fixed effects:
##                    Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)          2.4458     0.0732 15.0000   33.41  1.7e-15 ***
## Round               -0.4110     0.0277 14.0342  -14.82  5.8e-10 ***
## Round:NetworkType1  -0.1369     0.0533 14.0000   -2.57    0.022 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Round 
## Round       -0.437       
## Rnd:NtwrkT1  0.000  0.000

5. Bayesian Analysis of logDivergence2

In the below, we report the analysis of an alternative measure of divergence.

# d.quad.logdiv2.bm1.base <- brm(
#   data = d.quad,
#   logDivergence2 ~ 1 + Round + NetworkType +
#     (Round | Quad),
#   control = list(adapt_delta = 0.99, 
#                  max_treedepth = 15),
#   save_all_pars = TRUE,
#   prior = prior1.base,
#   iter = 5000, warmup = 2000,
#   chains = 4, cores = 4,
#   seed = 1000)
# 
# saveRDS(d.quad.logdiv2.bm1.base, file = "d-quad-logdiv2-bm1-base.Rds")
d.quad.logdiv2.bm1.base <- readRDS("d-quad-logdiv2-bm1-base.Rds")
# d.quad.logdiv2.bm1 <- brm(
#   data = d.quad,
#   logDivergence2 ~ 1 + Round + NetworkType:Round +
#     (Round | Quad),
#   control = list(adapt_delta = 0.95),
#   save_all_pars = TRUE,
#   prior = prior1,
#   iter = 5000, warmup = 2000,
#   chains = 4, cores = 4,
#   seed = 1000)
# 
# saveRDS(d.quad.logdiv2.bm1, file = "d-quad-logdiv2-bm1.Rds")
d.quad.logdiv2.bm1 <- readRDS("d-quad-logdiv2-bm1.Rds")
plot(d.quad.logdiv2.bm1, ask = FALSE, newpage = FALSE, N = 4)

print(summary(d.quad.logdiv2.bm1), digits = 3)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: logDivergence2 ~ 1 + Round + NetworkType:Round + (Round | Quad) 
##    Data: d.quad (Number of observations: 64) 
## Samples: 4 chains, each with iter = 5000; warmup = 2000; thin = 1;
##          total post-warmup samples = 12000
## 
## Group-Level Effects: 
## ~Quad (Number of levels: 16) 
##                      Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sd(Intercept)           0.302     0.094    0.137    0.508       4454 1.001
## sd(Round)               0.087     0.037    0.016    0.165       2295 1.001
## cor(Intercept,Round)   -0.140     0.363   -0.717    0.657       4819 1.001
## 
## Population-Level Effects: 
##                    Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## Intercept             2.482     0.099    2.283    2.674       8675 1.000
## Round                -0.465     0.032   -0.530   -0.402       8820 1.000
## Round:NetworkType1   -0.202     0.056   -0.313   -0.093       7793 1.000
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sigma    0.186     0.025    0.144    0.243       3271 1.001
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).
pp_check(d.quad.logdiv2.bm1, nsamples=100)

5.1. Bayes Factor

# set.seed(1005)
# BF10 <- c()
# for (i in 1:10) {
#   BF10 <- c(BF10, bayes_factor(d.quad.logdiv2.bm1, d.quad.logdiv2.bm1.base)$bf)
# }
# saveRDS(BF10, file = "d-quad-logdiv2-bm1-bf.Rds")
d.quad.logdiv2.bm1.BF10 <- readRDS(file = "d-quad-logdiv2-bm1-bf.Rds")
summary(d.quad.logdiv2.bm1.BF10)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   15.96   16.53   17.67   17.45   18.17   18.75

6. Frequentist Analysis of logDivergence2

d.quad.logdiv2.lmer1 <- lmer(
  logDivergence2 ~ 1 + Round + NetworkType:Round +
    (Round | Quad), 
  REML = TRUE,
  data = d.quad,
  control = lmerControl(
    optimizer ='optimx', optCtrl=list(method='L-BFGS-B'))
  )
print(summary(d.quad.logdiv2.lmer1, ddf = "Kenward-Roger"), digits = 3)
## Linear mixed model fit by REML. t-tests use Kenward-Roger's method [
## lmerModLmerTest]
## Formula: logDivergence2 ~ 1 + Round + NetworkType:Round + (Round | Quad)
##    Data: d.quad
## Control: 
## lmerControl(optimizer = "optimx", optCtrl = list(method = "L-BFGS-B"))
## 
## REML criterion at convergence: 19.4
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -2.029 -0.417  0.050  0.521  1.780 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  Quad     (Intercept) 0.09061  0.3010        
##           Round       0.00853  0.0923   -0.42
##  Residual             0.02934  0.1713        
## Number of obs: 64, groups:  Quad, 16
## 
## Fixed effects:
##                    Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)          2.4829     0.0917 15.0000   27.07  3.8e-14 ***
## Round               -0.4656     0.0300 14.1215  -15.52  2.9e-10 ***
## Round:NetworkType1  -0.1933     0.0513 14.0000   -3.77   0.0021 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Round 
## Round       -0.599       
## Rnd:NtwrkT1  0.000  0.000

7. Network Age

# see preprocess_data_entropy_networkage.R for preprocessing data
load('data-entropy-networkage.RData')

# add NetworkAge
d <-
  d %>%
  mutate(NetworkAge = ifelse(Window == 1, 2, 8)) %>%
  mutate(Window = as_factor(Window)) %>%
  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)
# Window 1: -0.5, Window 2: 0.5
contrasts(d$Window) <- c(-0.5, 0.5)
d.quad.w <-
  d %>%
  group_by(NetworkType, Quad, Window, NetworkAge) %>%
  summarise(Divergence1 = mean(Entropy),
            Nfeatures1 = mean(Nfeatures),
            Divergence2 = median(Entropy),
            Nfeatures2 = median(Nfeatures)) %>%
  ungroup() %>%
  mutate(logDivergence1 = log2(Divergence1),
         logDivergence2 = log2(Divergence2))
d.quad.w %>%
  ggplot(aes(x = NetworkAge, y = logDivergence1, group = Quad, color = NetworkType)) + 
  geom_line() + 
  scale_x_continuous(breaks = c(2,8), labels = c(2,8)) +
  labs(x = "Average Network Age", 
       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())

For each quad, we calculated community-level divergence measures using the gesures produced by the A-X dyads (where X can be B, C, or D) at network ages of 0, 2, 4 (average network age = 2, or Window 1) or at network ages of 6, 8, 10 (average network age = 8, or Window 2).

7.1. NetworkType x Window

prior10 <- prior(normal(0, 3), class = Intercept) +
  prior(normal(0, 1), class = b, coef = "Window1") +
  prior(normal(0, 1), class = b, coef = "NetworkType1") +
  prior(normal(0, 1), class = b, coef = "Window1:NetworkType1") +
  prior(student_t(3, 0, 5), class = sd) +
  prior(student_t(3, 0, 5), class = sigma)
# d.quad.logdiv1.w.bm1 <- brm(
#   data = d.quad.w,
#   logDivergence1 ~ 1 + Window + NetworkType + Window:NetworkType +
#     (1 | Quad),
#   control = list(adapt_delta = 0.95),
#   sample_prior = TRUE,
#   save_all_pars = TRUE,
#   prior = prior10,
#   iter = 5000, warmup = 2000,
#   chains = 4, cores = 4,
#   seed = 1000)
# 
# saveRDS(d.quad.logdiv1.w.bm1, file = "d-quad-logdiv1-w-bm1.Rds")
d.quad.logdiv1.w.bm1 <- readRDS("d-quad-logdiv1-w-bm1.Rds")
plot(d.quad.logdiv1.w.bm1, ask = FALSE, newpage = FALSE, N = 3)

print(summary(d.quad.logdiv1.w.bm1), digits=3)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: logDivergence1 ~ 1 + Window + NetworkType + Window:NetworkType + (1 | Quad) 
##    Data: d.quad.w (Number of observations: 32) 
## Samples: 4 chains, each with iter = 5000; warmup = 2000; thin = 1;
##          total post-warmup samples = 12000
## 
## Group-Level Effects: 
## ~Quad (Number of levels: 16) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sd(Intercept)    0.284     0.082    0.143    0.466       2638 1.002
## 
## Population-Level Effects: 
##                      Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## Intercept               1.607     0.081    1.443    1.771       3579 1.000
## Window1                -0.668     0.072   -0.812   -0.524      12011 1.000
## NetworkType1           -0.080     0.165   -0.404    0.251       3559 1.001
## Window1:NetworkType1    0.150     0.147   -0.140    0.437      11060 1.000
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sigma    0.202     0.044    0.138    0.311       2510 1.004
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).
pp_check(d.quad.logdiv1.w.bm1, nsamples = 100)

7.2. NetworkType x NetworkAge

Here we use NetworkAge (a continuous variable) instead of Window (a categorical variable) as a predictor. Similar to the main analysis, we introduce a constraint that log divergence should be same in both NetworkType conditions at NetworkAge of 0 (not observed).

prior11 <- prior(normal(0, 3), class = Intercept) +
  prior(normal(0, 1), class = b, coef = "NetworkAge") +
  prior(normal(0, 1), class = b, coef = "NetworkAge:NetworkType1") +
  prior(student_t(3, 0, 5), class = sd) +
  prior(student_t(3, 0, 5), class = sigma)
# d.quad.logdiv1.w.bm2 <- brm(
#   data = d.quad.w,
#   logDivergence1 ~ 1 + NetworkAge + NetworkAge:NetworkType +
#     (1 | Quad),
#   control = list(adapt_delta = 0.95),
#   sample_prior = TRUE,
#   save_all_pars = TRUE,
#   prior = prior11,
#   iter = 5000, warmup = 2000,
#   chains = 4, cores = 4,
#   seed = 1000)
# 
# saveRDS(d.quad.logdiv1.w.bm2, file = "d-quad-logdiv1-w-bm2.Rds")
d.quad.logdiv1.w.bm2 <- readRDS("d-quad-logdiv1-w-bm2.Rds")
plot(d.quad.logdiv1.w.bm2, ask = FALSE, newpage = FALSE, N = 3)

print(summary(d.quad.logdiv1.w.bm2), digits=3)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: logDivergence1 ~ 1 + NetworkAge + NetworkAge:NetworkType + (1 | Quad) 
##    Data: d.quad.w (Number of observations: 32) 
## Samples: 4 chains, each with iter = 5000; warmup = 2000; thin = 1;
##          total post-warmup samples = 12000
## 
## Group-Level Effects: 
## ~Quad (Number of levels: 16) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sd(Intercept)    0.283     0.081    0.135    0.464       2371 1.001
## 
## Population-Level Effects: 
##                         Estimate Est.Error l-95% CI u-95% CI Eff.Sample
## Intercept                  2.163     0.102    1.958    2.368       4754
## NetworkAge                -0.112     0.012   -0.136   -0.087      13612
## NetworkAge:NetworkType1    0.010     0.020   -0.030    0.048       5600
##                          Rhat
## Intercept               1.001
## NetworkAge              1.000
## NetworkAge:NetworkType1 1.001
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sigma    0.204     0.045    0.138    0.315       2576 1.002
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).
pp_check(d.quad.logdiv1.w.bm2, nsamples = 100)

sessionInfo()
## R version 3.6.2 (2019-12-12)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] optimx_2018-7.10 lmerTest_3.1-0   lme4_1.1-21      Matrix_1.2-18   
##  [5] brms_2.8.0       Rcpp_1.0.1       bayesplot_1.6.0  tidybayes_1.0.4 
##  [9] forcats_0.4.0    stringr_1.4.0    dplyr_0.8.0.1    purrr_0.3.2     
## [13] readr_1.3.1      tidyr_0.8.3      tibble_2.1.1     ggplot2_3.1.1   
## [17] tidyverse_1.2.1 
## 
## loaded via a namespace (and not attached):
##   [1] minqa_1.2.4               colorspace_1.4-1         
##   [3] ellipsis_0.2.0.1          ggridges_0.5.1           
##   [5] rsconnect_0.8.13          ggstance_0.3.1           
##   [7] markdown_0.9              base64enc_0.1-3          
##   [9] rstudioapi_0.10           rstan_2.18.2             
##  [11] svUnit_0.7-12             DT_0.5                   
##  [13] mvtnorm_1.0-10            lubridate_1.7.4          
##  [15] xml2_1.2.0                bridgesampling_0.6-0     
##  [17] splines_3.6.2             knitr_1.22               
##  [19] shinythemes_1.1.2         jsonlite_1.6             
##  [21] nloptr_1.2.1              pbkrtest_0.4-7           
##  [23] broom_0.5.2               shiny_1.3.2              
##  [25] compiler_3.6.2            httr_1.4.0               
##  [27] backports_1.1.4           assertthat_0.2.1         
##  [29] lazyeval_0.2.2            cli_1.1.0                
##  [31] later_0.8.0               htmltools_0.3.6          
##  [33] prettyunits_1.0.2         tools_3.6.2              
##  [35] igraph_1.2.4.1            coda_0.19-2              
##  [37] gtable_0.3.0              glue_1.3.1               
##  [39] reshape2_1.4.3            cellranger_1.1.0         
##  [41] nlme_3.1-143              crosstalk_1.0.0          
##  [43] xfun_0.6                  ps_1.3.0                 
##  [45] rvest_0.3.3               mime_0.6                 
##  [47] miniUI_0.1.1.1            gtools_3.8.1             
##  [49] MASS_7.3-51.4             zoo_1.8-5                
##  [51] scales_1.0.0              colourpicker_1.0         
##  [53] hms_0.4.2                 promises_1.0.1           
##  [55] Brobdingnag_1.2-6         parallel_3.6.2           
##  [57] inline_0.3.15             shinystan_2.5.0          
##  [59] yaml_2.2.0                gridExtra_2.3            
##  [61] loo_2.1.0                 StanHeaders_2.18.1       
##  [63] stringi_1.4.3             dygraphs_1.1.1.6         
##  [65] boot_1.3-24               pkgbuild_1.0.3           
##  [67] rlang_0.4.0               pkgconfig_2.0.2          
##  [69] matrixStats_0.54.0        evaluate_0.14            
##  [71] lattice_0.20-38           rstantools_1.5.1         
##  [73] htmlwidgets_1.3           labeling_0.3             
##  [75] tidyselect_0.2.5          processx_3.3.0           
##  [77] plyr_1.8.4                magrittr_1.5             
##  [79] R6_2.4.0                  generics_0.0.2           
##  [81] pillar_1.4.2              haven_2.1.0              
##  [83] withr_2.1.2               xts_0.11-2               
##  [85] abind_1.4-5               modelr_0.1.4             
##  [87] crayon_1.3.4              arrayhelpers_1.0-20160527
##  [89] rmarkdown_1.12            grid_3.6.2               
##  [91] readxl_1.3.1              callr_3.2.0              
##  [93] threejs_0.3.1             digest_0.6.20            
##  [95] xtable_1.8-4              numDeriv_2016.8-1        
##  [97] httpuv_1.5.1              stats4_3.6.2             
##  [99] munsell_0.5.0             shinyjs_1.0