# load libraries
library(tidyverse)
library(tidybayes)
library(bayesplot)
library(ggplot2)
library(brms)
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.
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:
\(\mu_{n,c,r} = (\beta_0 + b_{0,c}) + ((\beta_1 + b_{1,c}) + \beta_2 NT_{n}) R_{r}\)
\(\sigma \sim \text{half-Student-t}(3, 0, 5)\)
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.
\(\beta_0 \sim N(0, 3^2)\)
\(\beta_1 \sim N(0, 1^2)\)
\(\beta_2 \sim N(0, 1^2)\)
\(\begin{align} \begin{bmatrix} b_{0,q} \\ b_{1,q} \end{bmatrix} & \sim N(\begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} \sigma_{0}^2 & \rho_{01} \sigma_0 \sigma_1\\ \rho_{01} \sigma_0 \sigma_1 & \sigma_{1}^2 \end{bmatrix}) \end{align}\)
\(\sigma_0 \sim \text{half-Student-t}(3, 0, 5)\)
\(\sigma_1 \sim \text{half-Student-t}(3, 0, 5)\)
\(\begin{align} \begin{bmatrix} 1 & \rho_{01}\\ \rho_{01} & 1 \end{bmatrix} \end{align} \sim \text{LKJ}(2)\)
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}])\).
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)
# 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.
# 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
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
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)
# 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
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
# 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).
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)
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