Analyses

Set-up

Packages

library(brms)
library(ggplot2)
library(kableExtra)
library(lme4)
library(lmerTest)
library(rmarkdown)
library(naniar)
library(performance)
library(see)
library(sjmisc)
library(tidyverse)

options(
  digits = 3
)
set.seed(170819)

Custom functions

# function to silence brms output
hush <- 
  function(
    code
    ){
    sink("/dev/null")
    tmp = code
    sink()
    return(tmp)
    }

Data

d <- read_csv("data/data.csv")

# same as above; but original file name:
# d <- read_csv("data/DataAggregated_T1T2_costsbenefits.csv")

# load image for work in IDE
# load("data/image.RData")

d <- d |> 
  rename(
    group = roles,
    op_expr = n_OE,
    gender = DE01_T1,
    age = DE02_01_T1,
    pol_stance = DE06_01_T1
  ) |> 
  mutate(
    female = as.logical(2 - gender),
    gender = factor(gender, labels = c("female", "male"))
  )

# recode to make as sum coding
d$anonymity_dev <- factor(d$anonymity)
contrasts(d$anonymity_dev) <- contr.sum(2)
d$persistence_dev <- factor(d$persistence)
contrasts(d$persistence_dev) <- contr.sum(2)

Missing data

Let’s first inspect how much data is missing

# Filter variables used for analyses
d_filt <- d |> 
  select(age, gender, pol_stance, topic, persistence, persistence_dev, anonymity, anonymity_dev, op_expr)

# Summary
miss_var_summary(d_filt)       # % missing per variable

# A tibble: 9 × 3
  variable        n_miss pct_miss
  <chr>            <int>    <num>
1 age                  0        0
2 gender               0        0
3 pol_stance           0        0
4 topic                0        0
5 persistence          0        0
6 persistence_dev      0        0
7 anonymity            0        0
8 anonymity_dev        0        0
9 op_expr              0        0

miss_case_summary(d_filt)      # % missing per row

# A tibble: 960 × 3
    case n_miss pct_miss
   <int>  <int>    <dbl>
 1     1      0        0
 2     2      0        0
 3     3      0        0
 4     4      0        0
 5     5      0        0
 6     6      0        0
 7     7      0        0
 8     8      0        0
 9     9      0        0
10    10      0        0
# ℹ 950 more rows

# Visualization
vis_miss(d_filt)               # Heatmap of missingness

gg_miss_var(d_filt)            # Bar plot: missing per variable

gg_miss_case(d_filt)           # Bar plot: missing per case

Shows that there are nore missing data for the variables we’ll analyse. Hence, no imputing necessary.

Descriptives

Let’s inspect distribution of opinion expressions.

ggplot(d, aes(op_expr)) +
  geom_histogram(binwidth = 1)

Looks like a zero-inflated poisson distribution. Confirms our preregistered approach to analyze data using zero-inflated Poisson approach.

nrow(d[d$op_expr == 0, ]) / nrow(d)

[1] 0.214

Overall, 21% of participants without any opinion expressions.

Let’s look at distribution of experimental groups.

d |> 
  select(persistence, anonymity) |> 
  table()

           anonymity
persistence high low
       high  240 240
       low   240 240

Distribution among groups perfect.

d |> 
  select(topic) |> 
  table()

topic
  climate    gender migration 
      320       320       320 

Distribution of topics also perfect.

Let’s check if groups are nested in topics:

is_nested(d$topic, d$group)

'f2' is nested within 'f1'

[1] TRUE

Indeed the case.

We first look at the experimental group’s descriptives

d |> 
  group_by(persistence) |> 
  summarize(op_expr_m = mean(op_expr)) |> 
  as.data.frame() |> 
  kable()

persistence	op_expr_m
high	9.26
low	9.29

Looking at persistence, we see there’s virtually no difference among groups.

d |> 
  group_by(anonymity) |> 
  summarize(op_expr_m = mean(op_expr)) |> 
  as.data.frame() |> 
  kable()

anonymity	op_expr_m
high	9.03
low	9.52

People who with less anonymity communicated more. But the difference isn’t particularly large.

d |> 
  group_by(persistence, anonymity) |> 
  summarize(op_expr_m = mean(op_expr)) |> 
  as.data.frame() |> 
  kable()

`summarise()` has grouped output by 'persistence'. You can override using the
`.groups` argument.

persistence	anonymity	op_expr_m
high	high	9.27
high	low	9.25
low	high	8.79
low	low	9.78

Looking at both groups combined, we see that low anonymity and low persistence created highest participation. But differences among groups aren’t large.

d |> 
  group_by(group) |> 
  summarize(
    anonymity = anonymity[1],
    persistence = persistence[1],
    topic = topic[1],
    op_expr_m = mean(op_expr)
    ) |> 
  rmarkdown::paged_table()

Looking at the various individual groups, we do see some difference. Generally, this shows that communication varied across groups.

d |> 
  group_by(topic) |> 
  summarize(op_expr_m = mean(op_expr)) |> 
  as.data.frame() |> 
  kable()

topic	op_expr_m
climate	9.63
gender	9.96
migration	8.22

Looking at topics specifically, we also see that there’s some variance.

Opinion expressions

Let’s look at the distribution of communicated number of words.

ggplot(d, aes(op_expr)) +
  geom_histogram(bins = 50)

Let’s also look at how many did not express any opinions at all.

Number:

length(which(d$op_expr == 0))

[1] 205

Percentage:

length(which(d$op_expr == 0)) / length(d$op_expr)

[1] 0.214

Stats:

summary(d$op_expr)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0     1.0     6.0     9.3    12.0   186.0 

Communicated words

Let’s look at the distribution of communicated number of words.

ggplot(d, aes(n_Words)) +
  geom_histogram(bins = 50)

Warning: Removed 198 rows containing non-finite outside the scale range
(`stat_bin()`).

Let’s also look at how many did not write any words at all.

length(which(d$n_Words == 0))

[1] 0

Everyone wrote at least one word.

Stats:

summary(d$n_Words)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1     149     336     496     602    8843     198 

Manipulation Check

Let’s see if respondents indeed perceived the experimental manipulations. Let’s first look at descriptives.

d |> 
  group_by(persistence) |> 
  summarize(
    "Perceived persistence" = mean(per_persistence, na.rm = TRUE)
    ) |> 
  as.data.frame() |> 
  kable()

persistence	Perceived persistence
high	1.78
low	4.19

The experimental manipulation worked.

model_pers <- lm(
  per_persistence ~ persistence_dev,
  d
)

summary(model_pers)

Call:
lm(formula = per_persistence ~ persistence_dev, data = d)

Residuals:
   Min     1Q Median     3Q    Max 
-3.194 -0.777 -0.110  0.806  3.223 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)        2.9854     0.0378    79.1   <2e-16 ***
persistence_dev1  -1.2085     0.0378   -32.0   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1 on 703 degrees of freedom
  (255 observations deleted due to missingness)
Multiple R-squared:  0.593, Adjusted R-squared:  0.592 
F-statistic: 1.02e+03 on 1 and 703 DF,  p-value: <2e-16

The difference was statistically significant. Let’s now look at anonymity.

d |> 
  group_by(anonymity) |> 
  summarize(
    "Perceived anonymity" = mean(per_anonymity, na.rm = TRUE)
    ) |> 
  as.data.frame() |> 
  kable()

anonymity	Perceived anonymity
high	1.22
low	4.42

The experimental manipulation worked.

model_anon <- lm(
  per_anonymity ~ anonymity_dev,
  d
)
summary(model_anon)

Call:
lm(formula = per_anonymity ~ anonymity_dev, data = d)

Residuals:
   Min     1Q Median     3Q    Max 
-3.425 -0.216 -0.216  0.575  3.784 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)      2.8204     0.0348    81.0   <2e-16 ***
anonymity_dev1  -1.6042     0.0348   -46.1   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.924 on 703 degrees of freedom
  (255 observations deleted due to missingness)
Multiple R-squared:  0.751, Adjusted R-squared:  0.751 
F-statistic: 2.12e+03 on 1 and 703 DF,  p-value: <2e-16

The experimental manipulation worked.

Cross contamination

Let’s look at cross contamination, i.e. if anonymity affected perceived persistence and if persistence affects perceived anonymity.

d |> 
  group_by(anonymity) |> 
  summarize(
    "Perceived persistence" = mean(per_persistence, na.rm = TRUE)
    ) |> 
  as.data.frame() |> 
  kable()

anonymity	Perceived persistence
high	2.89
low	3.00

There was next to no difference among groups.

model_pers <- lm(
  per_persistence ~ anonymity_dev,
  d
)

summary(model_pers)

Call:
lm(formula = per_persistence ~ anonymity_dev, data = d)

Residuals:
   Min     1Q Median     3Q    Max 
-2.005 -1.671  0.108  1.662  2.108 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)      2.9485     0.0591   49.86   <2e-16 ***
anonymity_dev1  -0.0561     0.0591   -0.95     0.34    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.57 on 703 degrees of freedom
  (255 observations deleted due to missingness)
Multiple R-squared:  0.00128,   Adjusted R-squared:  -0.000141 
F-statistic: 0.901 on 1 and 703 DF,  p-value: 0.343

The difference was not statistically significant. No cross contamination re. anonymity.

d |> 
  group_by(persistence) |> 
  summarize(
    "Perceived anonymity" = mean(per_anonymity, na.rm = TRUE)
    ) |> 
  as.data.frame() |> 
  kable()

persistence	Perceived anonymity
high	2.76
low	2.93

There was next to no difference in the groups’ means.

model_anon <- lm(
  per_anonymity ~ persistence_dev,
  d
)
summary(model_anon)

Call:
lm(formula = per_anonymity ~ persistence_dev, data = d)

Residuals:
   Min     1Q Median     3Q    Max 
-1.933 -1.763 -0.763  2.067  2.237 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)        2.8479     0.0698   40.83   <2e-16 ***
persistence_dev1  -0.0848     0.0698   -1.22     0.22    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.85 on 703 degrees of freedom
  (255 observations deleted due to missingness)
Multiple R-squared:  0.0021,    Adjusted R-squared:  0.00068 
F-statistic: 1.48 on 1 and 703 DF,  p-value: 0.224

Again, no significant difference. In conclusion, no cross-contamination.

Random Allocation

Let’s check if random allocation across experimental conditions worked alright. Let’s first look at descriptives.

d |> 
  group_by(persistence, anonymity) |> 
  summarize(
    female_m = mean(female)
    , age_m = mean(age)
    , pol_stance_m = mean(pol_stance)
    ) |> 
  as.data.frame() |> 
  kable()

`summarise()` has grouped output by 'persistence'. You can override using the
`.groups` argument.

persistence	anonymity	female_m	age_m	pol_stance_m
high	high	0.613	43.2	5.57
high	low	0.662	39.6	5.44
low	high	0.558	41.7	5.76
low	low	0.683	40.8	5.15

There seem to be some differences. Let’s inspect manually.

model_persis <- lm(
  as.integer(persistence_dev) ~ age + gender + pol_stance,
  d
)
summary(model_persis)

Call:
lm(formula = as.integer(persistence_dev) ~ age + gender + pol_stance, 
    data = d)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.5352 -0.4994 -0.0006  0.4994  0.5385 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.536577   0.072667   21.15   <2e-16 ***
age         -0.000632   0.001402   -0.45     0.65    
gendermale   0.022885   0.034672    0.66     0.51    
pol_stance  -0.003460   0.008089   -0.43     0.67    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.501 on 956 degrees of freedom
Multiple R-squared:  0.000702,  Adjusted R-squared:  -0.00243 
F-statistic: 0.224 on 3 and 956 DF,  p-value: 0.88

Allocation of gender, age, and political stance on the two persistence groups was successful.

model_anon <- lm(
  as.integer(anonymity_dev) ~ age + gender + pol_stance,
  d
)
summary(model_anon)

Call:
lm(formula = as.integer(anonymity_dev) ~ age + gender + pol_stance, 
    data = d)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.6702 -0.4902  0.0125  0.4874  0.7263 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.77558    0.07192   24.69   <2e-16 ***
age         -0.00323    0.00139   -2.33   0.0201 *  
gendermale  -0.06686    0.03432   -1.95   0.0517 .  
pol_stance  -0.02142    0.00801   -2.68   0.0076 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.496 on 956 degrees of freedom
Multiple R-squared:  0.0211,    Adjusted R-squared:  0.018 
F-statistic: 6.87 on 3 and 956 DF,  p-value: 0.00014

However, allocation across anonymity groups wasn’t successful. In the subsequent analyses, let’s hence control for these sociodemographic variables.

Bayesian mixed effects modeling

We analyze the data using Bayesian modelling.

We use deviation/sum contrast coding (-.1, .1). Meaning, contrasts measure main effects of independent variables.

Fixed effects

We preregistered to analyze fixed effects.

fit_fe_1 <- 
  hush(
    brm(
      op_expr ~ 
        1 + persistence_dev * anonymity_dev + age + female + pol_stance +
        (1 | topic/group)
      , data = d
      , chains = 4
      , cores = 4
      , iter = 6000
      , warmup = 2000
      , family = zero_inflated_poisson()
      , control = list(
        adapt_delta = .95
        , max_treedepth = 12
        )
      , save_pars = save_pars(all = TRUE)
      , silent = 2
      )
  )

Warning: There were 214 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Shows some convergence warnings. Let’s inspect model.

plot(fit_fe_1, ask = FALSE)

Trace-plots look alright.

Let’s look at results.

summary(fit_fe_1)

Warning: There were 214 divergent transitions after warmup. Increasing
adapt_delta above 0.95 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

 Family: zero_inflated_poisson 
  Links: mu = log; zi = identity 
Formula: op_expr ~ 1 + persistence_dev * anonymity_dev + age + female + pol_stance + (1 | topic/group) 
   Data: d (Number of observations: 960) 
  Draws: 4 chains, each with iter = 6000; warmup = 2000; thin = 1;
         total post-warmup draws = 16000

Multilevel Hyperparameters:
~topic (Number of levels: 3) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.27      0.31     0.01     1.18 1.00     1935     4004

~topic:group (Number of levels: 48) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.40      0.05     0.32     0.50 1.00     2863     4835

Regression Coefficients:
                                Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept                           1.99      0.20     1.56     2.43 1.00
persistence_dev1                   -0.01      0.06    -0.12     0.11 1.00
anonymity_dev1                     -0.01      0.06    -0.13     0.10 1.00
age                                 0.01      0.00     0.01     0.01 1.00
femaleTRUE                         -0.00      0.02    -0.05     0.05 1.00
pol_stance                         -0.02      0.01    -0.03    -0.01 1.00
persistence_dev1:anonymity_dev1     0.02      0.06    -0.09     0.14 1.00
                                Bulk_ESS Tail_ESS
Intercept                           3751     2124
persistence_dev1                    2818     4655
anonymity_dev1                      2685     5035
age                                22907    13685
femaleTRUE                          9989    10431
pol_stance                         11560    10620
persistence_dev1:anonymity_dev1     2773     4109

Further Distributional Parameters:
   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
zi     0.21      0.01     0.19     0.24 1.00    10401     9549

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

No significant effect emerged.

Let’s inspect ICC

var_ratio_fe <- performance::variance_decomposition(
  fit_fe_1
  , by_group = TRUE)
var_ratio_fe

# Random Effect Variances and ICC

Conditioned on: all random effects

## Variance Ratio (comparable to ICC)
Ratio: 0.42  CI 95%: [-0.24 0.72]

## Variances of Posterior Predicted Distribution
Conditioned on fixed effects: 31.94  CI 95%: [15.38 68.22]
Conditioned on rand. effects: 55.01  CI 95%: [49.75 60.91]

## Difference in Variances
Difference: 22.97  CI 95%: [-13.03 40.39]

41.886 percent of variance in opinion expressions explained by both topics and groups.

Let’s visualize results to see what they exactly mean.

p <- plot(
  conditional_effects(
    fit_fe_1
    ), 
  ask = FALSE,
  plot = FALSE
  )

p_anon <- 
  p[["anonymity_dev"]] +
  xlab("Anonymity") +
  ylab("Opinion expression") +
  scale_x_discrete(
    limits = rev
     ) +
  scale_y_continuous(
    limits = c(5, 14)
    , breaks = c(6, 8, 10, 12, 14)
    )

p_pers <- 
  p[["persistence_dev"]] +
  xlab("Persistence") +
  ylab("Opinion expression") +
  scale_x_discrete(
    limits = rev
   ) +
  scale_y_continuous(
    limits = c(5, 14)
    , breaks = c(6, 8, 10, 12, 14)
    ) +
  theme(
    axis.title.y = element_blank()
    )

p_int <- 
  p[["persistence_dev:anonymity_dev"]] +
  xlab("Persistence") +
  scale_x_discrete(
    limits = rev
     ) +
  scale_color_discrete(
    labels = c("low", "high")
    ) +
  guides(
    fill = "none",
    color = guide_legend(
      title = "Anonymity"
      )
    ) +
  theme(
    axis.title.y = element_blank()
    ) +
  scale_y_continuous(
    limits = c(5, 14)
    , breaks = c(6, 8, 10, 12, 14)
    )

plot <- cowplot::plot_grid(
  p_anon, p_pers, p_int, 
  labels = c('A', 'B', "C"), 
  nrow = 1,
  rel_widths = c(2, 2, 3)
  )

plot

ggsave("figures/results.png", plot, width = 8, height = 4)

Shows that there are no main effects. There seems to be a (nonsignificant) interaction effect. In low persistence environment, anonymity is conducive to communication; in high it’s the opposite.

Let’s look at posteriors

p_1 <- 
  pp_check(fit_fe_1) + 
  labs(title = "Zero-inflated poisson")

Using 10 posterior draws for ppc type 'dens_overlay' by default.

p_1

The actual distribution cannot be precisely reproduced, but it’s also not too far off.

Random effects

We preregistered to explore and compare models with random effects. So let’s model how the experimental conditions affect the outcomes differently depending on topic.

fit_re_1 <- 
  hush(
    brm(
      op_expr ~ 
        1 + persistence_dev * anonymity_dev + age + female + pol_stance +
        (1 + persistence_dev * anonymity_dev | topic) + 
        (1 | topic:group)
      , data = d
      , chains = 4
      , cores = 4
      , iter = 6000
      , warmup = 2000
      , family = zero_inflated_poisson()
      , control = list(
        adapt_delta = .95
        , max_treedepth = 15
        )
      , save_pars = save_pars(all = TRUE)
    )
  )

Warning: There were 757 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess

Shows some convergence warnings.

Let’s inspect model.

plot(fit_re_1, ask = FALSE)

Traceplots look alright.

Let’s look at results.

summary(fit_re_1)

Warning: There were 757 divergent transitions after warmup. Increasing
adapt_delta above 0.95 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

 Family: zero_inflated_poisson 
  Links: mu = log; zi = identity 
Formula: op_expr ~ 1 + persistence_dev * anonymity_dev + age + female + pol_stance + (1 + persistence_dev * anonymity_dev | topic) + (1 | topic:group) 
   Data: d (Number of observations: 960) 
  Draws: 4 chains, each with iter = 6000; warmup = 2000; thin = 1;
         total post-warmup draws = 16000

Multilevel Hyperparameters:
~topic (Number of levels: 3) 
                                                      Estimate Est.Error
sd(Intercept)                                             0.35      0.45
sd(persistence_dev1)                                      0.25      0.38
sd(anonymity_dev1)                                        0.25      0.40
sd(persistence_dev1:anonymity_dev1)                       0.44      0.53
cor(Intercept,persistence_dev1)                          -0.01      0.46
cor(Intercept,anonymity_dev1)                            -0.01      0.46
cor(persistence_dev1,anonymity_dev1)                     -0.01      0.46
cor(Intercept,persistence_dev1:anonymity_dev1)            0.06      0.46
cor(persistence_dev1,persistence_dev1:anonymity_dev1)     0.01      0.46
cor(anonymity_dev1,persistence_dev1:anonymity_dev1)      -0.01      0.46
                                                      l-95% CI u-95% CI Rhat
sd(Intercept)                                             0.01     1.64 1.00
sd(persistence_dev1)                                      0.00     1.34 1.00
sd(anonymity_dev1)                                        0.00     1.36 1.00
sd(persistence_dev1:anonymity_dev1)                       0.02     2.00 1.00
cor(Intercept,persistence_dev1)                          -0.83     0.82 1.00
cor(Intercept,anonymity_dev1)                            -0.83     0.82 1.00
cor(persistence_dev1,anonymity_dev1)                     -0.83     0.84 1.00
cor(Intercept,persistence_dev1:anonymity_dev1)           -0.79     0.86 1.00
cor(persistence_dev1,persistence_dev1:anonymity_dev1)    -0.83     0.84 1.00
cor(anonymity_dev1,persistence_dev1:anonymity_dev1)      -0.84     0.84 1.00
                                                      Bulk_ESS Tail_ESS
sd(Intercept)                                             3999     7589
sd(persistence_dev1)                                      3022     1235
sd(anonymity_dev1)                                        2973     4156
sd(persistence_dev1:anonymity_dev1)                        812     1075
cor(Intercept,persistence_dev1)                          18914    11709
cor(Intercept,anonymity_dev1)                            11052    10770
cor(persistence_dev1,anonymity_dev1)                     14598    12328
cor(Intercept,persistence_dev1:anonymity_dev1)            4205     1403
cor(persistence_dev1,persistence_dev1:anonymity_dev1)    13374    12801
cor(anonymity_dev1,persistence_dev1:anonymity_dev1)       3222     1168

~topic:group (Number of levels: 48) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.40      0.05     0.32     0.50 1.00     1232     2167

Regression Coefficients:
                                Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept                           2.00      0.29     1.36     2.59 1.00
persistence_dev1                   -0.00      0.22    -0.43     0.47 1.00
anonymity_dev1                     -0.01      0.25    -0.45     0.45 1.00
age                                 0.01      0.00     0.01     0.01 1.00
femaleTRUE                         -0.00      0.02    -0.05     0.04 1.00
pol_stance                         -0.02      0.01    -0.03    -0.01 1.00
persistence_dev1:anonymity_dev1     0.05      0.39    -0.69     1.38 1.01
                                Bulk_ESS Tail_ESS
Intercept                           5679     5933
persistence_dev1                    7250     6100
anonymity_dev1                      4316     4669
age                                17356    12335
femaleTRUE                         25008     9823
pol_stance                         14780    10789
persistence_dev1:anonymity_dev1      681      124

Further Distributional Parameters:
   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
zi     0.21      0.01     0.19     0.24 1.00    16420    10765

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Again, no main or interaction effects.

Let’s see if the random effects model fits better

cores <- parallel::detectCores()

fit_fe_1 <- add_criterion(
  fit_fe_1
  , "kfold" 
  , K = 5
  , cores = cores
  )

Fitting model 1 out of 5

Fitting model 2 out of 5

Fitting model 3 out of 5

Fitting model 4 out of 5

Fitting model 5 out of 5

fit_re_1 <- add_criterion(
  fit_re_1
  , "kfold"
  , K = 5
  , cores = cores
  )

Fitting model 1 out of 5

Start sampling

Fitting model 2 out of 5

Start sampling

Fitting model 3 out of 5

Start sampling

Fitting model 4 out of 5

Start sampling

Fitting model 5 out of 5

Start sampling

comp_1 <- loo_compare(fit_fe_1, fit_re_1, criterion = "kfold")
comp_1

         elpd_diff se_diff
fit_re_1    0.0       0.0 
fit_fe_1 -149.4      69.3 

Although model comparisons showed that the model with random effects fitted better, the difference was not significant (Δ ELPD = -149.36, 95% CI [-285.246, -13.467]. Hence, for reasons of parsimony the model with fixed effects is preferred.

Null-Model

Let’s also inspect a model without random intercepts to see if including random intercepts is worthwhile.

fit_nm_1 <- 
  hush(
    brm(
      op_expr ~ 
        1 + persistence_dev * anonymity_dev + age + female + pol_stance
      , data = d
      , chains = 4
      , cores = 4
      , iter = 6000
      , warmup = 2000
      , family = zero_inflated_poisson()
      , control = list(
        adapt_delta = .95
        , max_treedepth = 12
        )
      , save_pars = save_pars(all = TRUE)
      , silent = 2
      )
  )

summary(fit_nm_1)

 Family: zero_inflated_poisson 
  Links: mu = log; zi = identity 
Formula: op_expr ~ 1 + persistence_dev * anonymity_dev + age + female + pol_stance 
   Data: d (Number of observations: 960) 
  Draws: 4 chains, each with iter = 6000; warmup = 2000; thin = 1;
         total post-warmup draws = 16000

Regression Coefficients:
                                Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept                           1.96      0.06     1.85     2.06 1.00
persistence_dev1                    0.00      0.01    -0.02     0.02 1.00
anonymity_dev1                     -0.04      0.01    -0.06    -0.02 1.00
age                                 0.01      0.00     0.01     0.01 1.00
femaleTRUE                          0.02      0.02    -0.02     0.07 1.00
pol_stance                         -0.01      0.01    -0.02     0.01 1.00
persistence_dev1:anonymity_dev1     0.02      0.01    -0.00     0.04 1.00
                                Bulk_ESS Tail_ESS
Intercept                          20193    13428
persistence_dev1                   20156    11001
anonymity_dev1                     20783    10932
age                                19306    11896
femaleTRUE                         20703    11806
pol_stance                         22715    12379
persistence_dev1:anonymity_dev1    20156    12002

Further Distributional Parameters:
   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
zi     0.21      0.01     0.19     0.24 1.00    20932    11802

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Shows significant effect of anonymity on opinion expression.

Let’s compare models.

fit_nm_1 <- add_criterion(
  fit_nm_1
  , "kfold"
  , K = 5
  , cores = cores
  )

Fitting model 1 out of 5

Fitting model 2 out of 5

Fitting model 3 out of 5

Fitting model 4 out of 5

Fitting model 5 out of 5

comp_2 <- loo_compare(fit_fe_1, fit_nm_1, criterion = "kfold")
comp_2

         elpd_diff se_diff
fit_fe_1    0.0       0.0 
fit_nm_1 -364.9     134.0 

The model comparisons showed that the model with random intercepts fitted significantly better than the null model with fixed intercepts (Δ ELPD = -364.85, 95% CI [-627.549, -102.152].

Hurdle

Let’s now estimate a fixed effects model with hurdles.

fit_hrdl_1 <- 
  hush(
    brm(
      bf(
        op_expr ~ 
          1 + persistence_dev * anonymity_dev + age + female + pol_stance +
          (1 | topic) + 
          (1 | topic:group),
        zi ~ 
          1 + persistence_dev * anonymity_dev + age + female + pol_stance +
          (1 | topic) + 
          (1 | topic:group)
      )
    , data = d
    , chains = 4
    , cores = 4
    , iter = 6000
    , warmup = 2000
    , family = zero_inflated_poisson()
    , control = list(
      adapt_delta = .95
      , max_treedepth = 15
      )
    )
  )

Warning: There were 422 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess

Again, some warnings.

Let’s inspect model.

plot(fit_hrdl_1, ask = FALSE)

Trace-plots look alright.

summary(fit_hrdl_1)

Warning: There were 422 divergent transitions after warmup. Increasing
adapt_delta above 0.95 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

 Family: zero_inflated_poisson 
  Links: mu = log; zi = logit 
Formula: op_expr ~ 1 + persistence_dev * anonymity_dev + age + female + pol_stance + (1 | topic) + (1 | topic:group) 
         zi ~ 1 + persistence_dev * anonymity_dev + age + female + pol_stance + (1 | topic) + (1 | topic:group)
   Data: d (Number of observations: 960) 
  Draws: 4 chains, each with iter = 6000; warmup = 2000; thin = 1;
         total post-warmup draws = 16000

Multilevel Hyperparameters:
~topic (Number of levels: 3) 
                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)        0.26      0.31     0.01     1.20 1.00     1220      292
sd(zi_Intercept)     0.28      0.42     0.01     1.52 1.00     3380     4173

~topic:group (Number of levels: 48) 
                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)        0.40      0.05     0.32     0.50 1.00     2278     1113
sd(zi_Intercept)     0.24      0.14     0.02     0.52 1.00     3893     7197

Regression Coefficients:
                                   Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept                              2.00      0.18     1.60     2.39 1.01
zi_Intercept                          -1.69      0.49    -2.61    -0.72 1.00
persistence_dev1                      -0.01      0.06    -0.13     0.11 1.00
anonymity_dev1                        -0.01      0.06    -0.13     0.11 1.00
age                                    0.01      0.00     0.01     0.01 1.00
femaleTRUE                            -0.00      0.02    -0.05     0.05 1.00
pol_stance                            -0.02      0.01    -0.03    -0.01 1.00
persistence_dev1:anonymity_dev1        0.02      0.06    -0.09     0.16 1.01
zi_persistence_dev1                    0.03      0.09    -0.15     0.20 1.00
zi_anonymity_dev1                      0.01      0.09    -0.17     0.18 1.00
zi_age                                 0.01      0.01    -0.00     0.03 1.00
zi_femaleTRUE                          0.19      0.17    -0.15     0.53 1.00
zi_pol_stance                         -0.05      0.04    -0.12     0.04 1.00
zi_persistence_dev1:anonymity_dev1    -0.01      0.09    -0.18     0.17 1.00
                                   Bulk_ESS Tail_ESS
Intercept                              4517     2523
zi_Intercept                           8977     5131
persistence_dev1                       3100     2940
anonymity_dev1                         2296     5882
age                                   17399    12354
femaleTRUE                            19167    11310
pol_stance                            21511    10355
persistence_dev1:anonymity_dev1        1098      208
zi_persistence_dev1                   12625     3777
zi_anonymity_dev1                     15515    11855
zi_age                                17492    11736
zi_femaleTRUE                         19399    11393
zi_pol_stance                         19998    11223
zi_persistence_dev1:anonymity_dev1    15922    10687

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Same results, no main effects, slightly larger but still non-significant interaction effect.

Exploratory Analyses

Frequentist

Look at results from a frequentist perspective.

Fixed effects

Estimate nested model.

fit_fe_1_frq <- 
  lmer(
    op_expr ~ 
      1 + 
      (1 | topic/group) + 
      persistence_dev * anonymity_dev + age + female + pol_stance
    , data = d
    )

summary(fit_fe_1_frq)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: op_expr ~ 1 + (1 | topic/group) + persistence_dev * anonymity_dev +  
    age + female + pol_stance
   Data: d

REML criterion at convergence: 7604

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.484 -0.555 -0.196  0.264 13.549 

Random effects:
 Groups      Name        Variance Std.Dev.
 group:topic (Intercept) 1.08e+01 3.28e+00
 topic       (Intercept) 8.92e-08 2.99e-04
 Residual                1.55e+02 1.25e+01
Number of obs: 960, groups:  group:topic, 48; topic, 3

Fixed effects:
                                Estimate Std. Error       df t value Pr(>|t|)
(Intercept)                       6.0450     2.1478 885.1185    2.81    0.005
persistence_dev1                 -0.0197     0.6214  43.8644   -0.03    0.975
anonymity_dev1                   -0.3437     0.6242  44.6401   -0.55    0.585
age                               0.0860     0.0357 943.4320    2.41    0.016
femaleTRUE                       -0.2434     0.8783 939.5946   -0.28    0.782
pol_stance                       -0.0319     0.2058 942.6445   -0.16    0.877
persistence_dev1:anonymity_dev1   0.1953     0.6226  44.1939    0.31    0.755
                                  
(Intercept)                     **
persistence_dev1                  
anonymity_dev1                    
age                             * 
femaleTRUE                        
pol_stance                        
persistence_dev1:anonymity_dev1   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) prss_1 anny_1 age    fmTRUE pl_stn
prsstnc_dv1  0.015                                   
annymty_dv1  0.053  0.000                            
age         -0.750 -0.010 -0.049                     
femaleTRUE  -0.463 -0.014  0.041  0.252              
pol_stance  -0.538 -0.009 -0.057 -0.003  0.062       
prsstn_1:_1  0.019  0.000 -0.002 -0.045 -0.034  0.038
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

Quite weird that topic doesn’t get any variance at all. Perhaps due to small cluster size? With Bayesian estimation, it worked alright. Also, again no significant effects.

Estimate without nesting.

fit_fe_2_frq <- 
  lmer(
    op_expr ~ 
      1 + 
      (1 | group) +
      persistence_dev * anonymity_dev + age + female + pol_stance + topic
    , data = d
    )

summary(fit_fe_2_frq)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: op_expr ~ 1 + (1 | group) + persistence_dev * anonymity_dev +  
    age + female + pol_stance + topic
   Data: d

REML criterion at convergence: 7598

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.513 -0.548 -0.188  0.265 13.577 

Random effects:
 Groups   Name        Variance Std.Dev.
 group    (Intercept)  11       3.32   
 Residual             155      12.46   
Number of obs: 960, groups:  group, 48

Fixed effects:
                                Estimate Std. Error       df t value Pr(>|t|)
(Intercept)                       6.3760     2.3261 511.8117    2.74   0.0063
persistence_dev1                 -0.0195     0.6257  41.8742   -0.03   0.9753
anonymity_dev1                   -0.3439     0.6284  42.6053   -0.55   0.5871
age                               0.0854     0.0357 942.0312    2.39   0.0169
femaleTRUE                       -0.2641     0.8788 938.1103   -0.30   0.7638
pol_stance                       -0.0324     0.2058 941.3972   -0.16   0.8751
topicgender                       0.4373     1.5334  41.9702    0.29   0.7769
topicmigration                   -1.3078     1.5329  41.9180   -0.85   0.3984
persistence_dev1:anonymity_dev1   0.1960     0.6268  42.1850    0.31   0.7561
                                  
(Intercept)                     **
persistence_dev1                  
anonymity_dev1                    
age                             * 
femaleTRUE                        
pol_stance                        
topicgender                       
topicmigration                    
persistence_dev1:anonymity_dev1   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) prss_1 anny_1 age    fmTRUE pl_stn tpcgnd tpcmgr
prsstnc_dv1  0.014                                                 
annymty_dv1  0.049  0.000                                          
age         -0.701 -0.010 -0.049                                   
femaleTRUE  -0.423 -0.014  0.041  0.252                            
pol_stance  -0.489 -0.009 -0.057 -0.004  0.063                     
topicgender -0.325  0.000 -0.001  0.017 -0.024 -0.020              
topicmigrtn -0.337  0.000 -0.001  0.024  0.000 -0.016  0.500       
prsstn_1:_1  0.018  0.000 -0.001 -0.045 -0.033  0.038 -0.001 -0.002

Also shows no significant effects.

For curiosity, estimate also without hierarchical structure.

fit_fe_3_frq <- 
  lm(
    op_expr ~ 
      1 + 
      persistence_dev * anonymity_dev + topic + age + female + pol_stance
    , data = d
    )

summary(fit_fe_3_frq)

Call:
lm(formula = op_expr ~ 1 + persistence_dev * anonymity_dev + 
    topic + age + female + pol_stance, data = d)

Residuals:
   Min     1Q Median     3Q    Max 
-12.95  -7.46  -2.92   3.00 177.59 

Coefficients:
                                Estimate Std. Error t value Pr(>|t|)   
(Intercept)                       5.7503     2.2095    2.60   0.0094 **
persistence_dev1                 -0.0223     0.4148   -0.05   0.9571   
anonymity_dev1                   -0.3598     0.4191   -0.86   0.3908   
topicgender                       0.4291     1.0174    0.42   0.6733   
topicmigration                   -1.3114     1.0166   -1.29   0.1974   
age                               0.0901     0.0362    2.49   0.0129 * 
femaleTRUE                       -0.2010     0.8932   -0.23   0.8220   
pol_stance                        0.0398     0.2087    0.19   0.8489   
persistence_dev1:anonymity_dev1   0.2004     0.4166    0.48   0.6306   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.8 on 951 degrees of freedom
Multiple R-squared:  0.0115,    Adjusted R-squared:  0.00315 
F-statistic: 1.38 on 8 and 951 DF,  p-value: 0.202

Also here, no significant effects.

Gender

As preregistered, let’s see if effects differ across genders.

fit_fe_gen <- 
  hush(
    brm(
      op_expr ~ 
        1 + persistence_dev * anonymity_dev * gender + age + pol_stance +
        (1 | topic/group)
      , data = d
      , chains = 4
      , cores = 4
      , iter = 8000
      , warmup = 2000
      , family = zero_inflated_poisson()
      , control = list(
        adapt_delta = .95
        , max_treedepth = 12
        )
      )
  )

Warning: There were 603 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Again, some warnings.

Let’s inspect model.

plot(fit_fe_gen, ask = FALSE)

Traceplots look alright.

Let’s look at results.

summary(fit_fe_gen)

Warning: There were 603 divergent transitions after warmup. Increasing
adapt_delta above 0.95 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

 Family: zero_inflated_poisson 
  Links: mu = log; zi = identity 
Formula: op_expr ~ 1 + persistence_dev * anonymity_dev * gender + age + pol_stance + (1 | topic/group) 
   Data: d (Number of observations: 960) 
  Draws: 4 chains, each with iter = 8000; warmup = 2000; thin = 1;
         total post-warmup draws = 24000

Multilevel Hyperparameters:
~topic (Number of levels: 3) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.21      0.21     0.01     0.84 1.00     3091     4051

~topic:group (Number of levels: 48) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.40      0.05     0.32     0.50 1.00     3874     8425

Regression Coefficients:
                                           Estimate Est.Error l-95% CI u-95% CI
Intercept                                      1.99      0.15     1.67     2.32
persistence_dev1                               0.02      0.06    -0.10     0.14
anonymity_dev1                                -0.04      0.06    -0.16     0.08
gendermale                                    -0.00      0.02    -0.05     0.04
age                                            0.01      0.00     0.01     0.01
pol_stance                                    -0.02      0.01    -0.03    -0.01
persistence_dev1:anonymity_dev1                0.01      0.06    -0.11     0.12
persistence_dev1:gendermale                   -0.08      0.02    -0.12    -0.03
anonymity_dev1:gendermale                      0.07      0.02     0.03     0.11
persistence_dev1:anonymity_dev1:gendermale     0.05      0.02     0.00     0.09
                                           Rhat Bulk_ESS Tail_ESS
Intercept                                  1.00     5113     5226
persistence_dev1                           1.00     3885     6451
anonymity_dev1                             1.00     4466     6784
gendermale                                 1.00    12495     8754
age                                        1.00    16372     7194
pol_stance                                 1.00    15721    15008
persistence_dev1:anonymity_dev1            1.00     4269     7151
persistence_dev1:gendermale                1.00    13393    12807
anonymity_dev1:gendermale                  1.00    15059    13599
persistence_dev1:anonymity_dev1:gendermale 1.00    14762    13666

Further Distributional Parameters:
   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
zi     0.21      0.01     0.19     0.24 1.00    14835    13287

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Indeed, several gender effects.

• For females, the effect of persistence is larger, that is more positive.
• For females, the effect of anonymity is smaller, that is more negative.
• For females, the interaction effect is also a bit smaller, that is more negative.

Let’s visualize results.

p_gen <- plot(
  conditional_effects(
    fit_fe_gen
    ), 
  ask = FALSE,
  plot = FALSE
  )

p_gen_pers <- 
  p_gen[["persistence_dev:gender"]] +
  xlab("Persistence") +
  ylab("Opinion expression") +
  scale_y_continuous(
    limits = c(4, 15),
    breaks = c(5, 7.5, 10, 12.5, 15)
  ) +
  scale_x_discrete(
    limits = rev
  ) +
  guides(
    fill = "none"
    , color = "none"
    )

p_gen_anon <- 
  p_gen[["anonymity_dev:gender"]] +
  xlab("Anonymity") +
  ylab("Opinion expression") +
  scale_y_continuous(
    limits = c(3.5, 15),
    breaks = c(5, 7.5, 10, 12.5, 15)
  ) +
  theme(
    axis.title.y = element_blank()
    ) +
  guides(
    fill = "none"
    ) + 
  scale_x_discrete(
    limits = rev
  ) +
  scale_color_discrete(
    name = "Gender"
    )

plot_gen <- cowplot::plot_grid(
  p_gen_pers, p_gen_anon, 
  labels = c('A', 'B'), 
  nrow = 1,
  rel_widths = c(4, 5)
  )

plot_gen

ggsave("figures/results_gen.png", plot_gen, width = 8, height = 4)

Benefits

Let’s see if benefits differ across experimental groups.

We first look at the experimental group’s descriptives

d |> 
  group_by(persistence) |> 
  summarize(benefits_m = mean(benefits, na.rm = TRUE)) |> 
  as.data.frame() |> 
  kable()

persistence	benefits_m
high	3.12
low	3.23

Looking at persistence, we see people with lower persistence reporting slightly higher benefits.

d |> 
  group_by(anonymity) |> 
  summarize(benefits_m = mean(benefits, na.rm = TRUE)) |> 
  as.data.frame() |> 
  kable()

anonymity	benefits_m
high	3.15
low	3.20

Looking at anonymity, we see people with low anonymity reporting marginally higher benefits.

d |> 
  group_by(persistence, anonymity) |> 
  summarize(benefits_m = mean(benefits, na.rm = T)) |> 
  as.data.frame() |> 
  kable()

`summarise()` has grouped output by 'persistence'. You can override using the
`.groups` argument.

persistence	anonymity	benefits_m
high	high	3.07
high	low	3.18
low	high	3.22
low	low	3.23

Looking at both groups combined, we see that low anonymity and low persistence yielded highest benefits.

Let’s look if effects are significant.

fit_fe_ben_1 <- 
  hush(
    brm(
      benefits ~ 
        1 + persistence_dev * anonymity_dev  + age + female + pol_stance +
        (1 | topic/group)
      , data = d
      , chains = 4
      , cores = 4
      , iter = 6000
      , warmup = 2000
      , control = list(
        adapt_delta = .95
        , max_treedepth = 12
        )
      )
  )

Warning: Rows containing NAs were excluded from the model.

Warning: There were 160 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Let’s inspect model.

plot(fit_fe_ben_1, ask = FALSE)

Traceplots look alright.

Let’s look at results.

summary(fit_fe_ben_1)

Warning: There were 160 divergent transitions after warmup. Increasing
adapt_delta above 0.95 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: benefits ~ 1 + persistence_dev * anonymity_dev + age + female + pol_stance + (1 | topic/group) 
   Data: d (Number of observations: 705) 
  Draws: 4 chains, each with iter = 6000; warmup = 2000; thin = 1;
         total post-warmup draws = 16000

Multilevel Hyperparameters:
~topic (Number of levels: 3) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.13      0.17     0.00     0.63 1.00     2141     2071

~topic:group (Number of levels: 48) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.08      0.05     0.00     0.17 1.00     3801     6162

Regression Coefficients:
                                Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept                           3.21      0.17     2.85     3.55 1.00
persistence_dev1                   -0.05      0.03    -0.11     0.01 1.00
anonymity_dev1                     -0.03      0.03    -0.09     0.03 1.00
age                                -0.00      0.00    -0.01     0.00 1.00
femaleTRUE                         -0.07      0.06    -0.19     0.04 1.00
pol_stance                          0.02      0.01    -0.01     0.04 1.00
persistence_dev1:anonymity_dev1    -0.02      0.03    -0.08     0.04 1.00
                                Bulk_ESS Tail_ESS
Intercept                           2614      903
persistence_dev1                   10311     8192
anonymity_dev1                      9308     9903
age                                18813    12042
femaleTRUE                         13006     9174
pol_stance                         10421    10613
persistence_dev1:anonymity_dev1    10541     7780

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.74      0.02     0.70     0.78 1.00    12561     9668

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

No significant effects. But note that effect of persistence on perceived benefits only marginally not significant.

Costs

Let’s see if perceived differed across experimental groups.

We first look at the experimental group’s descriptives

d |> 
  group_by(persistence) |> 
  summarize(costs = mean(costs, na.rm = TRUE)) |> 
  as.data.frame() |> 
  kable()

persistence	costs
high	1.99
low	1.99

Looking at persistence, we see both groups report equal costs.

d |> 
  group_by(anonymity) |> 
  summarize(costs = mean(costs, na.rm = TRUE)) |> 
  as.data.frame() |> 
  kable()

anonymity	costs
high	1.89
low	2.09

Looking at anonymity, we see people with low anonymity report slightly higher costs.

d |> 
  group_by(persistence, anonymity) |> 
  summarize(costs = mean(costs, na.rm = TRUE)) |> 
  as.data.frame() |> 
  kable()

`summarise()` has grouped output by 'persistence'. You can override using the
`.groups` argument.

persistence	anonymity	costs
high	high	1.90
high	low	2.07
low	high	1.87
low	low	2.11

Looking at both groups combined, we see that highest costs were reported by group with low anonymity and low persistence.

Let’s look if effects are significant.

fit_fe_costs_1 <- 
  hush(
    brm(
      costs ~ 
        1 + persistence_dev * anonymity_dev + age + female + pol_stance +
        (1 | topic/group)
      , data = d
      , chains = 4
      , cores = 4
      , iter = 8000
      , warmup = 2000
      , control = list(
        adapt_delta = .95
        , max_treedepth = 12
        )
      )
  )

Warning: Rows containing NAs were excluded from the model.

Warning: There were 186 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Let’s inspect model.

plot(fit_fe_costs_1, ask = FALSE)

Traceplots look alright.

Let’s look at results.

summary(fit_fe_costs_1)

Warning: There were 186 divergent transitions after warmup. Increasing
adapt_delta above 0.95 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: costs ~ 1 + persistence_dev * anonymity_dev + age + female + pol_stance + (1 | topic/group) 
   Data: d (Number of observations: 705) 
  Draws: 4 chains, each with iter = 8000; warmup = 2000; thin = 1;
         total post-warmup draws = 24000

Multilevel Hyperparameters:
~topic (Number of levels: 3) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.13      0.19     0.00     0.66 1.00     3819     2147

~topic:group (Number of levels: 48) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.07      0.05     0.00     0.17 1.00     7146     9973

Regression Coefficients:
                                Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept                           2.48      0.19     2.11     2.86 1.00
persistence_dev1                    0.00      0.03    -0.07     0.07 1.00
anonymity_dev1                     -0.08      0.03    -0.15    -0.02 1.00
age                                -0.01      0.00    -0.02    -0.01 1.00
femaleTRUE                          0.01      0.07    -0.12     0.14 1.00
pol_stance                         -0.00      0.02    -0.04     0.03 1.00
persistence_dev1:anonymity_dev1     0.02      0.03    -0.05     0.09 1.00
                                Bulk_ESS Tail_ESS
Intercept                           5196     2440
persistence_dev1                   26569    16662
anonymity_dev1                     25889    17303
age                                34343    18126
femaleTRUE                         25883    16648
pol_stance                         34809    17089
persistence_dev1:anonymity_dev1    24338    13993

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.84      0.02     0.80     0.89 1.00    30590    16079

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

We find that anonymity does reduce costs.

Mediation

Let’s see if perceived benefits and costs were associated with increased opinion expressions.

fit_fe_med <- 
  hush(
    brm(
      op_expr ~ 
        1 + persistence_dev * anonymity_dev + benefits + costs  + age + female + pol_stance + 
        (1 | topic/group)
      , data = d
      , chains = 4
      , cores = 4
      , iter = 6000
      , warmup = 2000
      , family = zero_inflated_poisson()
      , control = list(
        adapt_delta = .95
        , max_treedepth = 12
        )
      )
  )

Warning: Rows containing NAs were excluded from the model.

Warning: There were 195 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Let’s look at results.

summary(fit_fe_med)

Warning: There were 195 divergent transitions after warmup. Increasing
adapt_delta above 0.95 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

 Family: zero_inflated_poisson 
  Links: mu = log; zi = identity 
Formula: op_expr ~ 1 + persistence_dev * anonymity_dev + benefits + costs + age + female + pol_stance + (1 | topic/group) 
   Data: d (Number of observations: 705) 
  Draws: 4 chains, each with iter = 6000; warmup = 2000; thin = 1;
         total post-warmup draws = 16000

Multilevel Hyperparameters:
~topic (Number of levels: 3) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.19      0.23     0.01     0.86 1.00     2218     4466

~topic:group (Number of levels: 48) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.41      0.05     0.32     0.51 1.00     3002     4917

Regression Coefficients:
                                Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept                           1.96      0.16     1.63     2.30 1.00
persistence_dev1                   -0.03      0.06    -0.14     0.09 1.00
anonymity_dev1                      0.00      0.06    -0.12     0.12 1.00
benefits                            0.11      0.02     0.08     0.14 1.00
costs                              -0.09      0.01    -0.12    -0.07 1.00
age                                 0.01      0.00     0.01     0.01 1.00
femaleTRUE                          0.00      0.02    -0.05     0.05 1.00
pol_stance                         -0.02      0.01    -0.03    -0.00 1.00
persistence_dev1:anonymity_dev1     0.02      0.06    -0.10     0.14 1.00
                                Bulk_ESS Tail_ESS
Intercept                           4690     4632
persistence_dev1                    2325     3690
anonymity_dev1                      2431     3821
benefits                           12168    11167
costs                              11977    10166
age                                21747    12851
femaleTRUE                         12484    10183
pol_stance                         11991     9162
persistence_dev1:anonymity_dev1     2662     4144

Further Distributional Parameters:
   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
zi     0.08      0.01     0.06     0.10 1.00    12415    10719

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

We find that increased perceived costs are associated with decreased opinion expressions. Increased benefits are associated with increased opinion expressions. Let’s check if overall effect is significant.

anon_costs_a_b <- fixef(fit_fe_costs_1)["anonymity_dev1", "Estimate"]
anon_costs_a_se <- fixef(fit_fe_costs_1)["anonymity_dev1", "Est.Error"]
anon_costs_a_dis <- rnorm(10000, anon_costs_a_b, anon_costs_a_se)

anon_costs_b_b <- fixef(fit_fe_med)["benefits", "Estimate"]
anon_costs_b_se <- fixef(fit_fe_med)["benefits", "Est.Error"]
anon_costs_b_dis <- rnorm(10000, anon_costs_b_b, anon_costs_b_se)

anon_costs_ab_dis <- anon_costs_a_dis * anon_costs_b_dis
anon_costs_ab_m <- median(anon_costs_ab_dis)
anon_costs_ab_ll <- quantile(anon_costs_ab_dis, .025)
anon_costs_ab_ul <- quantile(anon_costs_ab_dis, .975)

The effect is significant (b = -0.01, 95% MC CI [-0.02, 0]).

Save

save.image("data/image.RData")