Load data

vcs <- rio::import("data_complete_2021_zscored.rds")

var_label(vcs$dominance) <- "Dominance"
var_label(vcs$neuro) <- "Neuroticism"
var_label(vcs$agree) <- "Agreeableness"
var_label(vcs$extra) <- "Extraversion"
var_label(vcs$openn) <- "Openness"
var_label(vcs$consc) <- "Conscientiousness"
var_label(vcs$soir_full) <- "Unrestricted sociosexuality"
var_label(vcs$f0) <- "Voice pitch"
var_label(vcs$pf) <- "Formants"
vcs$sex_c <- vcs$sex
contrasts(vcs$sex) <- contr.helmert(2)
var_label(vcs$sex) <- "Sex"
set.seed(1)
var_label(vcs$age) <- "Age"
vcs <- vcs %>% 
  mutate(
    age = if_else(dataset == 9, 20, age)/10,
    age_se = if_else(dataset == 9, 3, 0.5)/10
  )
vcs <- vcs %>% filter(!is.na(f0), !is.na(pf), !is.na(age))

warmup <- 2000
iter <- warmup + 2000
chains <- 4
control <- list(adapt_delta = 0.99)
priors <- c(
  prior(normal(0, 3), class = b)
)
library(brms)
options(mc.cores = parallel::detectCores(), brms.backend = "cmdstanr")
rstan::rstan_options(auto_write = TRUE)

variable_labels <- c("Intercept"= "Intercept",  "sex1" = "Sex [male]", "bsp_meageage_se" = "Age±SE", "f0" = "Voice pitch (f0)", "pf" = " Formant position (Pf)", "age" = "Age", "f0:sex1" = "Voice pitch (f0) × Sex")
effect_labels <- c("b_Intercept"= "Intercept",  "b_sex1" = "Sex [male]", "bsp_meageage_se" = "Age±SE", "b_f0" = "Voice pitch (f0)", "b_pf" = " Formant position (Pf)", "b_age" = "Age", "b_f0.sex1" = "Voice pitch (f0) × Sex")

Hypothesis 1: Dominance

Participants with lower voice pitch will have a more dominant personality.

Interpretation: The data are are consistent with substantial linear negative relationships between f0 and dominance. The 89% HDI for f0 falls entirely outside the ROPE, but the HDI for pf falls almost entirely within it. This evidence is consistent with a non-negligible association, where people with deeper voices are more dominant (after adjusting for age and gender). Further research is needed to verify if the association with pf is truly negligible in size.

Visually diagnose non-linearity

ggplot(vcs, aes(f0, dominance)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 

ggplot(vcs, aes(pf, dominance)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 

ggplot(vcs, aes(age*10, dominance)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  scale_x_continuous("Age") +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

Decisions upon visual diagnosis: Visual diagnosis shows no clear sign of nonlinearities. We will fit linear effects for f0, pf, and age. It is not necessary to model measurement error in age, as dataset 9 (where age was not noted individually) is not included in these analyses. Visually, there could be an interaction between sex and f0.

Models

h1_simple <- brm(dominance ~ f0 + pf + sex + age + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, file_refit = "on_change", save_mevars = TRUE, 
          file = "models/dominance/h1_simple") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Warning: Rows containing NAs were excluded from the model.
## Warning: Argument 'save_mevars' is deprecated. Please use argument 'latent' in
## function 'save_pars()' instead.
h1_after_diagnosis <- brm(dominance ~ f0*sex + pf + sex + age + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, cores = chains, chains = chains,
          prior = priors, 
          control = control, file_refit = "on_change", save_mevars = TRUE,
          file = "models/dominance/h1_nonlinear") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Warning: Rows containing NAs were excluded from the model.

## Warning: Argument 'save_mevars' is deprecated. Please use argument 'latent' in
## function 'save_pars()' instead.
# conditional_smooths(h1_after_diagnosis)

compare_models <- loo_compare(h1_simple, h1_after_diagnosis)
compare_models
##                    elpd_diff se_diff
## h1_simple           0.0       0.0   
## h1_after_diagnosis -1.1       0.3
h1_sex_mod_dataset_varying <- brm(dominance ~ f0 + sex + age + 
                                    (1 + f0 + sex || dataset), data = vcs,
          prior = priors,
          iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
          control = control, file_refit = "on_change", save_mevars = TRUE,
          file = "models/dominance/h1_sex_mod_dataset_varying2") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Warning: Rows containing NAs were excluded from the model.

## Warning: Argument 'save_mevars' is deprecated. Please use argument 'latent' in
## function 'save_pars()' instead.
h1_after_diagnosis_melsm <- brm(
  bf(dominance ~ f0 + sex + age + (1 | dataset),
     sigma ~ f0 + sex + age + (1 | dataset)), data = vcs,
          prior = priors,
          iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
          control = control, file_refit = "on_change", save_mevars = TRUE,
          file = "models/dominance/h1_after_diagnosis_melsm") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Warning: Rows containing NAs were excluded from the model.

## Warning: Argument 'save_mevars' is deprecated. Please use argument 'latent' in
## function 'save_pars()' instead.
compare_models <- loo_compare(h1_simple, h1_after_diagnosis, h1_sex_mod_dataset_varying,
                      h1_after_diagnosis_melsm)
compare_models
##                            elpd_diff se_diff
## h1_after_diagnosis_melsm    0.0       0.0   
## h1_simple                  -1.7       3.2   
## h1_sex_mod_dataset_varying -2.1       3.3   
## h1_after_diagnosis         -2.8       3.2
r2s <- list(h1_simple=h1_simple, h1_after_diagnosis=h1_after_diagnosis, h1_sex_mod_dataset_varying=h1_sex_mod_dataset_varying,
                      h1_after_diagnosis_melsm=h1_after_diagnosis_melsm) %>% 
  map(~ { x <- as_tibble(bayes_R2(.))
          x$loo_R2 <- loo_R2(.)
          x
          }) %>%
  bind_rows(.id="model") %>% 
  arrange(desc(loo_R2))

sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
                  pred.labels = variable_labels)
  dominance
Predictors Estimates CI (95%)
Intercept -0.07 -0.80;0.68
Voice pitch (f0) -0.27 -0.45;-0.08
Formant position (Pf) -0.00 -0.13;0.12
Sex [male] -0.30 -0.52;-0.08
Age 0.04 -0.08;0.17
Random Effects
σ2 0.05
Ï„00 0.96
ICC 0.05
N dataset 4
Observations 985
Marginal R2 / Conditional R2 0.067 / 0.067 
equiv_test <- equivalence_test(h1_simple, ci = 0.89, range = c(-0.1, 0.1))
equiv_test_table <- equiv_test %>% 
  mutate(Parameter = recode(Parameter, !!!effect_labels)) %>% 
  mutate(
    ROPE_Percentage = sprintf("%2.f%%", 100*ROPE_Percentage),
    HDI = sprintf("[% .2f;% .2f]", HDI_low, HDI_high)) %>% 
  select(Parameter, H0 = ROPE_Equivalence, `inside ROPE` = ROPE_Percentage, `89% HDI` = HDI)
kable(equiv_test_table, caption = "Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).")
Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).
Parameter H0 inside ROPE 89% HDI
Intercept Undecided 31% [-0.61; 0.45]
Voice pitch (f0) Rejected 0% [-0.42;-0.12]
Formant position (Pf) Undecided 98% [-0.10; 0.11]
Sex [male] Rejected 0% [-0.47;-0.12]
Age Undecided 85% [-0.05; 0.14] 
plot(equiv_test) + scale_y_discrete("Effect", breaks = names(effect_labels), labels = effect_labels)
## Warning: Removed 3516 rows containing non-finite values (stat_density_ridges).

Model details 
summary(h1_simple)
## Warning: There were 3 divergent transitions after warmup. Increasing adapt_delta
## above may help. See http://mc-stan.org/misc/warnings.html#divergent-transitions-
## after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: dominance ~ f0 + pf + sex + age + (1 | dataset) 
##    Data: vcs (Number of observations: 985) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.53      0.41     0.16     1.68 1.00     1993     2686
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.07      0.37    -0.80     0.68 1.00     2616     2533
## f0           -0.27      0.09    -0.45    -0.08 1.00     5357     5343
## pf           -0.00      0.07    -0.13     0.12 1.00     6495     5685
## sex1         -0.30      0.11    -0.52    -0.08 1.00     4465     4603
## age           0.04      0.06    -0.08     0.17 1.00     7482     5631
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.97      0.02     0.93     1.02 1.00     8271     5405
## 
## Samples were drawn using sample(hmc). 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).
summary(h1_after_diagnosis)
## Warning: There were 4 divergent transitions after warmup. Increasing adapt_delta
## above may help. See http://mc-stan.org/misc/warnings.html#divergent-transitions-
## after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: dominance ~ f0 * sex + pf + sex + age + (1 | dataset) 
##    Data: vcs (Number of observations: 985) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.50      0.36     0.16     1.49 1.00     2341     3510
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.06      0.35    -0.80     0.60 1.00     3154     3304
## f0           -0.27      0.10    -0.46    -0.07 1.00     4870     4934
## sex1         -0.29      0.12    -0.52    -0.06 1.00     4091     4809
## pf           -0.00      0.06    -0.13     0.12 1.00     6343     5184
## age           0.04      0.06    -0.08     0.16 1.00     6846     5390
## f0:sex1       0.03      0.10    -0.17     0.22 1.00     7135     5738
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.97      0.02     0.93     1.02 1.00     8237     4988
## 
## Samples were drawn using sample(hmc). 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).
conditional_effects(h1_after_diagnosis_melsm)

# conditional_smooths(h1_after_diagnosis)
summary(h1_sex_mod_dataset_varying)
## Warning: There were 3 divergent transitions after warmup. Increasing adapt_delta
## above may help. See http://mc-stan.org/misc/warnings.html#divergent-transitions-
## after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: dominance ~ f0 + sex + age + (1 + f0 + sex || dataset) 
##    Data: vcs (Number of observations: 985) 
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup samples = 12000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.61      0.49     0.17     1.94 1.00     4104     5242
## sd(f0)            0.19      0.27     0.00     0.84 1.00     3248     3164
## sd(sex1)          0.32      0.42     0.01     1.48 1.00     3078     4575
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.11      0.42    -0.91     0.70 1.00     4814     4447
## f0           -0.27      0.19    -0.59     0.07 1.00     5027     3851
## sex1         -0.29      0.28    -0.84     0.30 1.00     4763     3834
## age           0.04      0.06    -0.08     0.16 1.00    14892     8393
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.97      0.02     0.93     1.01 1.00    15665     8678
## 
## Samples were drawn using sample(hmc). 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).
summary(h1_after_diagnosis_melsm)
## Warning: There were 5 divergent transitions after warmup. Increasing adapt_delta
## above may help. See http://mc-stan.org/misc/warnings.html#divergent-transitions-
## after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: dominance ~ f0 + sex + age + (1 | dataset) 
##          sigma ~ f0 + sex + age + (1 | dataset)
##    Data: vcs (Number of observations: 985) 
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup samples = 12000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 4) 
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)           0.54      0.42     0.17     1.68 1.00     3531     3988
## sd(sigma_Intercept)     0.11      0.13     0.00     0.43 1.00     2675     3580
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept          -0.12      0.35    -0.81     0.54 1.00     4919     3908
## sigma_Intercept     0.06      0.15    -0.23     0.33 1.00     5784     6023
## f0                 -0.30      0.10    -0.49    -0.11 1.00     8218     8432
## sex1               -0.33      0.10    -0.52    -0.14 1.00     8427     7752
## age                 0.05      0.06    -0.07     0.17 1.00    11668     8593
## sigma_f0            0.11      0.07    -0.03     0.24 1.00     8502     8442
## sigma_sex1          0.11      0.07    -0.02     0.25 1.00     8474     8383
## sigma_age          -0.02      0.05    -0.11     0.07 1.00     7440     8750
## 
## Samples were drawn using sample(hmc). 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).
kable(r2s, caption = "Explained variance, regular and leave-one-out.")
## Warning in `[<-.data.frame`(`*tmp*`, , isn, value = structure(list(Estimate =
## structure(c("0.0682743", : provided 8 variables to replace 5 variables
Explained variance, regular and leave-one-out.
model Estimate Est.Error Q2.5 Q97.5 loo_R2
h1_after_diagnosis_melsm 0.0682743 0.0146551 0.0412222 0.0987603 0.05358670
h1_simple 0.0678976 0.0145653 0.0409913 0.0978596 0.05217783
h1_sex_mod_dataset_varying 0.0696528 0.0147511 0.0427012 0.0998353 0.05090080
h1_after_diagnosis 0.0688751 0.0147137 0.0418394 0.0993437 0.05040002 
kable(coef(h1_sex_mod_dataset_varying)$dataset[, ,"f0"], caption = "Estimated association between f0 and outcome by dataset")
Estimated association between f0 and outcome by dataset
Estimate Est.Error Q2.5 Q97.5
2 -0.3044166 0.1171668 -0.5477592 -0.0850778
3 -0.2324883 0.1275848 -0.4676571 0.0382604
7 -0.2415646 0.1575541 -0.5270493 0.1070431
8 -0.2991568 0.1479890 -0.6182300 -0.0178871

Hypothesis 2: Extraversion

Participants with lower voice pitch will score higher on extraversion

Interpretation: The data are are consistent with a substantial linear negative relationship between f0 and extraversion. The 89% HDI for f0 excludes and falls entirely outside the ROPE. The 89% of the HDI for pf overlaps zero and the ROPE boundaries. This evidence is consistent with an association, where people with deeper voices have higher extraversion (after adjusting for age and gender), but further research is needed to test whether the association with pf is negligible in size.

Visually diagnose non-linearity

ggplot(vcs, aes(f0, extra)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour =  "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 

ggplot(vcs, aes(pf, extra)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 

ggplot(vcs, aes(age*10, extra)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  scale_x_continuous("Age") +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

Decisions upon visual diagnosis: Visual diagnostic shows signs of a potential interaction by sex coupled with nonlinearity for pf. We will fit separate splines by sex for pf. It is necessary to model measurement error in age, as dataset 9 (where age was not noted individually) is included in these analyses.

Models

h1_simple <- brm(extra ~ f0 + pf + sex + me(age, age_se) + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, file_refit = "on_change", save_mevars = TRUE, 
          file = "models/extra/h1_simple") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Warning: Rows containing NAs were excluded from the model.
## Warning: Argument 'save_mevars' is deprecated. Please use argument 'latent' in
## function 'save_pars()' instead.
h1_after_diagnosis <- brm(extra ~ f0 + s(pf, by = sex) + sex + me(age, age_se) + (1 | dataset), data = vcs,
          iter = iter + 1000, warmup = warmup + 1000, cores = chains, chains = chains,
          prior = priors, 
          control = control, file_refit = "on_change", save_mevars = TRUE,
          file = "models/extra/h1_nonlinear2") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Warning: Rows containing NAs were excluded from the model.

## Warning: Argument 'save_mevars' is deprecated. Please use argument 'latent' in
## function 'save_pars()' instead.
conditional_smooths(h1_after_diagnosis)

compare_models <- loo_compare(h1_simple, h1_after_diagnosis)
compare_models
##                    elpd_diff se_diff
## h1_simple           0.0       0.0   
## h1_after_diagnosis -2.1       0.8
h1_sex_mod_dataset_varying <- brm(extra ~ f0 + sex + me(age, age_se) + 
                                    (1 + f0 + sex || dataset), data = vcs,
          prior = priors,
          iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
          control = control, file_refit = "on_change", save_mevars = TRUE,
          file = "models/extra/h1_sex_mod_dataset_varying2") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Warning: Rows containing NAs were excluded from the model.

## Warning: Argument 'save_mevars' is deprecated. Please use argument 'latent' in
## function 'save_pars()' instead.
h1_after_diagnosis_melsm <- brm(
  bf(extra ~ f0 + sex + me(age, age_se) + (1 | dataset),
     sigma ~ f0 + sex + me(age, age_se) + (1 | dataset)), data = vcs,
          prior = priors,
          iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
          control = control, file_refit = "on_change", save_mevars = TRUE,
          file = "models/extra/h1_after_diagnosis_melsm") %>%
  add_criterion("loo") %>%
  add_criterion("bayes_R2") %>%
  add_criterion("loo_R2")
## Warning: Rows containing NAs were excluded from the model.

## Warning: Argument 'save_mevars' is deprecated. Please use argument 'latent' in
## function 'save_pars()' instead.
compare_models <- loo_compare(h1_simple, h1_after_diagnosis, h1_sex_mod_dataset_varying ,h1_after_diagnosis_melsm
                      )
compare_models
##                            elpd_diff se_diff
## h1_after_diagnosis_melsm     0.0       0.0  
## h1_simple                  -13.8       5.4  
## h1_sex_mod_dataset_varying -14.7       5.4  
## h1_after_diagnosis         -15.8       5.4
r2s <- list(h1_simple=h1_simple, h1_after_diagnosis=h1_after_diagnosis, h1_sex_mod_dataset_varying=h1_sex_mod_dataset_varying,h1_after_diagnosis_melsm=h1_after_diagnosis_melsm) %>% 
  map(~ { x <- as_tibble(bayes_R2(.))
          x$loo_R2 <- loo_R2(.)
          x
          }) %>%
  bind_rows(.id="model") %>% 
  arrange(desc(loo_R2))

sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
                  pred.labels = variable_labels)
  extra
Predictors Estimates CI (95%)
Intercept 0.19 -0.21;0.58
Voice pitch (f0) -0.23 -0.38;-0.08
Formant position (Pf) 0.08 -0.03;0.18
Sex [male] -0.28 -0.46;-0.10
Age±SE -0.11 -0.22;0.01
Random Effects
σ2 0.03
Ï„00 0.97
ICC 0.03
N dataset 7
Observations 1433
Marginal R2 / Conditional R2 0.056 / 0.056 
equiv_test <- equivalence_test(h1_simple, ci = 0.89, range = c(-0.1, 0.1))
equiv_test_table <- equiv_test %>% 
  mutate(Parameter = recode(Parameter, !!!effect_labels)) %>% 
  mutate(
    ROPE_Percentage = sprintf("%2.f%%", 100*ROPE_Percentage),
    HDI = sprintf("[% .2f;% .2f]", HDI_low, HDI_high)) %>% 
  select(Parameter, H0 = ROPE_Equivalence, `inside ROPE` = ROPE_Percentage, `89% HDI` = HDI)
kable(equiv_test_table, caption = "Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).")
Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).
Parameter H0 inside ROPE 89% HDI
Intercept Undecided 28% [-0.14; 0.50]
Voice pitch (f0) Rejected 0% [-0.35;-0.11]
Formant position (Pf) Undecided 68% [-0.00; 0.16]
Sex [male] Rejected 0% [-0.43;-0.14]
Age±SE Undecided 46% [-0.20;-0.01] 
plot(equiv_test) + scale_y_discrete("Effect", breaks = names(effect_labels), labels = effect_labels)
## Warning: Removed 3516 rows containing non-finite values (stat_density_ridges).

Model details 
summary(h1_simple)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: extra ~ f0 + pf + sex + me(age, age_se) + (1 | dataset) 
##    Data: vcs (Number of observations: 1433) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.35      0.15     0.16     0.73 1.00     2404     4372
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept       0.19      0.20    -0.21     0.58 1.00     3114     4163
## f0             -0.23      0.08    -0.38    -0.08 1.00     7221     6134
## pf              0.08      0.05    -0.03     0.18 1.00     8287     6673
## sex1           -0.28      0.09    -0.46    -0.10 1.00     6190     6003
## meageage_se    -0.11      0.06    -0.22     0.01 1.00     7134     5697
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.97      0.02     0.94     1.01 1.00    12826     5752
## 
## Samples were drawn using sample(hmc). 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).
summary(h1_after_diagnosis)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: extra ~ f0 + s(pf, by = sex) + sex + me(age, age_se) + (1 | dataset) 
##    Data: vcs (Number of observations: 1433) 
## Samples: 4 chains, each with iter = 5000; warmup = 3000; thin = 1;
##          total post-warmup samples = 8000
## 
## Smooth Terms: 
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(spfsexfemale_1)     0.68      0.65     0.02     2.37 1.00     5523     4613
## sds(spfsexmale_1)       0.84      0.76     0.02     2.86 1.00     5205     4347
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.35      0.15     0.16     0.74 1.00     3172     5300
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept           0.15      0.21    -0.25     0.56 1.00     3981     4779
## f0                 -0.23      0.08    -0.38    -0.08 1.00    10304     6090
## sex1               -0.30      0.10    -0.49    -0.11 1.00     9329     6016
## spf:sexfemale_1     1.04      1.51    -2.16     4.22 1.00     6769     5245
## spf:sexmale_1       0.50      1.66    -2.70     3.99 1.00     7882     5473
## meageage_se        -0.11      0.06    -0.22     0.01 1.00    10385     5509
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.97      0.02     0.94     1.01 1.00    17508     5732
## 
## Samples were drawn using sample(hmc). 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).
conditional_effects(h1_after_diagnosis_melsm)

# conditional_smooths(h1_after_diagnosis_melsm)
summary(h1_sex_mod_dataset_varying)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: extra ~ f0 + sex + me(age, age_se) + (1 + f0 + sex || dataset) 
##    Data: vcs (Number of observations: 1433) 
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup samples = 12000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.37      0.17     0.17     0.81 1.00     4274     5566
## sd(f0)            0.09      0.09     0.00     0.31 1.00     4374     4828
## sd(sex1)          0.12      0.11     0.00     0.41 1.00     3729     5006
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept       0.18      0.22    -0.23     0.60 1.00     4786     5527
## f0             -0.23      0.09    -0.41    -0.05 1.00     8654     7548
## sex1           -0.32      0.12    -0.53    -0.08 1.00     6527     6153
## meageage_se    -0.10      0.06    -0.22     0.01 1.00     9621     9231
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.97      0.02     0.94     1.01 1.00    17997     7853
## 
## Samples were drawn using sample(hmc). 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).
summary(h1_after_diagnosis_melsm)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: extra ~ f0 + sex + me(age, age_se) + (1 | dataset) 
##          sigma ~ f0 + sex + me(age, age_se) + (1 | dataset)
##    Data: vcs (Number of observations: 1433) 
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup samples = 12000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)           0.37      0.16     0.17     0.77 1.00     3166     5080
## sd(sigma_Intercept)     0.14      0.07     0.05     0.31 1.00     3043     5691
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept             0.18      0.20    -0.20     0.58 1.00     3126     4796
## sigma_Intercept       0.04      0.12    -0.22     0.27 1.00     4058     5424
## f0                   -0.20      0.08    -0.35    -0.05 1.00     6137     7650
## sex1                 -0.32      0.08    -0.48    -0.17 1.00     5940     7828
## sigma_f0              0.09      0.06    -0.02     0.20 1.00     6680     7842
## sigma_sex1            0.12      0.06     0.00     0.23 1.00     6574     7975
## meageage_se          -0.10      0.05    -0.20    -0.00 1.00     7974     8486
## sigma_meageage_se    -0.01      0.04    -0.10     0.08 1.00     5065     7343
## 
## Samples were drawn using sample(hmc). 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).
kable(r2s, caption = "Explained variance, regular and leave-one-out.")
## Warning in `[<-.data.frame`(`*tmp*`, , isn, value = structure(list(Estimate =
## structure(c("0.0560665", : provided 8 variables to replace 5 variables
Explained variance, regular and leave-one-out.
model Estimate Est.Error Q2.5 Q97.5 loo_R2
h1_simple 0.0560665 0.0115298 0.0347388 0.0802184 0.04263608
h1_after_diagnosis_melsm 0.0529681 0.0112895 0.0321689 0.0762337 0.04255709
h1_sex_mod_dataset_varying 0.0574352 0.0114818 0.0365357 0.0813209 0.04124206
h1_after_diagnosis 0.0579205 0.0115332 0.0364228 0.0818447 0.03981800 
kable(coef(h1_sex_mod_dataset_varying)$dataset[, ,"f0"], caption = "Estimated association between f0 and outcome by dataset")
Estimated association between f0 and outcome by dataset
Estimate Est.Error Q2.5 Q97.5
2 -0.2000055 0.0933028 -0.3742787 -0.0116282
3 -0.1965949 0.1002789 -0.3872111 0.0046707
4 -0.2098497 0.1084163 -0.4160660 0.0126459
7 -0.2143766 0.1160460 -0.4378958 0.0328257
8 -0.2703890 0.1278774 -0.5753656 -0.0530679
9 -0.2426635 0.1038801 -0.4564030 -0.0422624
10 -0.2544931 0.1234924 -0.5380544 -0.0361293 
equivalence_test(h1_after_diagnosis_melsm, ci = 0.89, range = c(-0.1, 0.1))
## # Test for Practical Equivalence
## 
##   ROPE: [-0.10 0.10]
## 
## Parameter         |        H0 | inside ROPE |       89% HDI
## -----------------------------------------------------------
## Intercept         | Undecided |     29.54 % | [-0.15  0.47]
## sigma_Intercept   | Undecided |     62.68 % | [-0.15  0.24]
## f0                | Undecided |      4.70 % | [-0.32 -0.08]
## sex1              |  Rejected |      0.00 % | [-0.45 -0.20]
## sigma_f0          | Undecided |     54.06 % | [ 0.00  0.18]
## sigma_sex1        | Undecided |     38.65 % | [ 0.02  0.21]
## meageage_se       | Undecided |     48.70 % | [-0.18 -0.02]
## sigma_meageage_se |  Accepted |    100.00 % | [-0.08  0.06]

Hypothesis 3: Agreeableness

Participants with lower voice pitch will score higher on agreeableness.

Interpretation: The credible intervals for f0 and pf include zero, but the 89% HDI does not fall entirely within the ROPE. This evidence is consistent with a negligible association between f0/pf and agreeableness. Further research is needed to verify if the associations are truly negligible in size.

Visually diagnose non-linearity

ggplot(vcs, aes(f0, agree)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour =  "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 

ggplot(vcs, aes(pf, agree)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 

ggplot(vcs, aes(age*10, agree)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  scale_x_continuous("Age") +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

Decisions upon visual diagnosis: Visual diagnostic shows signs of nonlinearity for f0 and pf, but no clear interactions by sex. We will fit splines for f0, and pf. It is necessary to model measurement error in age, as dataset 9 (where age was not noted individually) is included in these analyses.

Models

h1_simple <- brm(agree ~ f0 + pf + sex + me(age, age_se) + (1 | dataset), data = vcs,
          iter = iter + 1000, warmup = warmup + 1000, chains = chains, cores = chains,
          prior = priors,
          control = control, file_refit = "on_change", save_mevars = TRUE, 
          file = "models/agree/h1_simple2") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Warning: Rows containing NAs were excluded from the model.
## Warning: Argument 'save_mevars' is deprecated. Please use argument 'latent' in
## function 'save_pars()' instead.
h1_after_diagnosis <- brm(agree ~ f0 * sex + pf + sex + me(age, age_se) + (1 | dataset), data = vcs,
                          iter = iter + 1500, warmup = warmup + 1500, cores = chains, chains = chains,
                          prior = priors, 
                          control = control, file_refit = "on_change", save_mevars = TRUE,
                          file = "models/agree/h1_nonlinear4") %>% 
    add_criterion("loo") %>% 
    add_criterion("bayes_R2") %>% 
    add_criterion("loo_R2") 
## Warning: Rows containing NAs were excluded from the model.

## Warning: Argument 'save_mevars' is deprecated. Please use argument 'latent' in
## function 'save_pars()' instead.
compare_models <- loo_compare(h1_simple, h1_after_diagnosis)
compare_models
##                    elpd_diff se_diff
## h1_after_diagnosis  0.0       0.0   
## h1_simple          -2.8       2.5
h1_sex_mod_dataset_varying <- brm(agree ~ f0 * sex + sex + me(age, age_se) + 
                                    (1 + f0 + sex || dataset), data = vcs,
          prior = priors,
          iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
          control = control, file_refit = "on_change", save_mevars = TRUE,
          file = "models/agree/h1_sex_mod_dataset_varying2") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Warning: Rows containing NAs were excluded from the model.

## Warning: Argument 'save_mevars' is deprecated. Please use argument 'latent' in
## function 'save_pars()' instead.
h1_after_diagnosis_melsm <- brm(
  bf(agree ~ f0 * sex + sex + me(age, age_se) + (1 | dataset),
     sigma ~ f0 * sex + sex + me(age, age_se) + (1 | dataset)), data = vcs,
          prior = priors,
          iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
          control = control, file_refit = "on_change", save_mevars = TRUE,
          file = "models/agree/h1_after_diagnosis_melsm2") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Warning: Rows containing NAs were excluded from the model.

## Warning: Argument 'save_mevars' is deprecated. Please use argument 'latent' in
## function 'save_pars()' instead.
compare_models <- loo_compare(h1_simple, h1_after_diagnosis, h1_sex_mod_dataset_varying,
                      h1_after_diagnosis_melsm)
compare_models
##                            elpd_diff se_diff
## h1_sex_mod_dataset_varying  0.0       0.0   
## h1_after_diagnosis_melsm   -1.0       4.3   
## h1_after_diagnosis         -3.1       2.6   
## h1_simple                  -5.9       3.5
r2s <- list(h1_simple=h1_simple, h1_after_diagnosis=h1_after_diagnosis, h1_sex_mod_dataset_varying=h1_sex_mod_dataset_varying,
                      h1_after_diagnosis_melsm=h1_after_diagnosis_melsm) %>% 
  map(~ { x <- as_tibble(bayes_R2(.))
          x$loo_R2 <- loo_R2(.)
          x
          }) %>%
  bind_rows(.id="model") %>% 
  arrange(desc(loo_R2))

sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
                  pred.labels = variable_labels)
  agree
Predictors Estimates CI (95%)
Intercept 0.09 -0.31;0.50
Voice pitch (f0) 0.03 -0.12;0.18
Formant position (Pf) 0.03 -0.07;0.13
Sex [male] -0.00 -0.18;0.18
Age±SE -0.06 -0.17;0.05
Random Effects
σ2 0.08
Ï„00 0.92
ICC 0.08
N dataset 7
Observations 1434
Marginal R2 / Conditional R2 0.093 / 0.093 
equiv_test <- equivalence_test(h1_simple, ci = 0.89, range = c(-0.1, 0.1))
equiv_test_table <- equiv_test %>% 
  mutate(Parameter = recode(Parameter, !!!effect_labels)) %>% 
  mutate(
    ROPE_Percentage = sprintf("%2.f%%", 100*ROPE_Percentage),
    HDI = sprintf("[% .2f;% .2f]", HDI_low, HDI_high)) %>% 
  select(Parameter, H0 = ROPE_Equivalence, `inside ROPE` = ROPE_Percentage, `89% HDI` = HDI)
kable(equiv_test_table, caption = "Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).")
Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).
Parameter H0 inside ROPE 89% HDI
Intercept Undecided 40% [-0.23; 0.41]
Voice pitch (f0) Undecided 87% [-0.09; 0.15]
Formant position (Pf) Undecided 94% [-0.05; 0.12]
Sex [male] Undecided 81% [-0.14; 0.15]
Age±SE Undecided 81% [-0.15; 0.04] 
plot(equiv_test) + scale_y_discrete("Effect", breaks = names(effect_labels), labels = effect_labels)
## Warning: Removed 3516 rows containing non-finite values (stat_density_ridges).

Model details 
summary(h1_simple)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: agree ~ f0 + pf + sex + me(age, age_se) + (1 | dataset) 
##    Data: vcs (Number of observations: 1434) 
## Samples: 4 chains, each with iter = 5000; warmup = 3000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.37      0.16     0.19     0.78 1.00     2594     4030
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept       0.09      0.21    -0.31     0.50 1.00     3516     4128
## f0              0.03      0.08    -0.12     0.18 1.00     8110     5574
## pf              0.03      0.05    -0.07     0.13 1.00     9227     6191
## sex1           -0.00      0.09    -0.18     0.18 1.00     6400     5786
## meageage_se    -0.06      0.06    -0.17     0.05 1.00     8322     5632
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.96      0.02     0.92     0.99 1.00    14174     6030
## 
## Samples were drawn using sample(hmc). 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).
summary(h1_after_diagnosis)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: agree ~ f0 * sex + pf + sex + me(age, age_se) + (1 | dataset) 
##    Data: vcs (Number of observations: 1434) 
## Samples: 4 chains, each with iter = 5500; warmup = 3500; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.38      0.16     0.19     0.80 1.00     3595     4732
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept       0.33      0.23    -0.10     0.79 1.00     4759     4862
## f0              0.09      0.08    -0.07     0.24 1.00    11065     6589
## sex1            0.09      0.10    -0.11     0.28 1.00     9094     5941
## pf              0.03      0.05    -0.07     0.13 1.00    13323     6413
## f0:sex1         0.20      0.08     0.05     0.35 1.00    13795     6043
## meageage_se    -0.07      0.06    -0.18     0.04 1.00    11475     6074
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.95      0.02     0.92     0.99 1.00    17519     5558
## 
## Samples were drawn using sample(hmc). 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).
conditional_effects(h1_after_diagnosis_melsm)

# conditional_smooths(h1_after_diagnosis_melsm)
summary(h1_sex_mod_dataset_varying)
## Warning: There were 1 divergent transitions after warmup. Increasing adapt_delta
## above may help. See http://mc-stan.org/misc/warnings.html#divergent-transitions-
## after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: agree ~ f0 * sex + sex + me(age, age_se) + (1 + f0 + sex || dataset) 
##    Data: vcs (Number of observations: 1434) 
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup samples = 12000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.42      0.19     0.20     0.91 1.00     4088     5358
## sd(f0)            0.16      0.12     0.01     0.46 1.00     3418     4195
## sd(sex1)          0.13      0.13     0.00     0.47 1.00     3872     5220
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept       0.27      0.25    -0.21     0.77 1.00     5354     6301
## f0              0.06      0.11    -0.15     0.27 1.00     8556     6721
## sex1            0.04      0.13    -0.21     0.27 1.00     7186     5599
## f0:sex1         0.15      0.09    -0.03     0.31 1.00    11444     9722
## meageage_se    -0.06      0.06    -0.17     0.05 1.00    13844     9028
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.95      0.02     0.92     0.99 1.00    19640     8421
## 
## Samples were drawn using sample(hmc). 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).
summary(h1_after_diagnosis_melsm)
##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: agree ~ f0 * sex + sex + me(age, age_se) + (1 | dataset) 
##          sigma ~ f0 * sex + sex + me(age, age_se) + (1 | dataset)
##    Data: vcs (Number of observations: 1434) 
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup samples = 12000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)           0.38      0.16     0.19     0.78 1.00     3487     5587
## sd(sigma_Intercept)     0.04      0.04     0.00     0.13 1.00     4576     5330
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept             0.31      0.22    -0.11     0.75 1.00     3967     5651
## sigma_Intercept       0.18      0.12    -0.05     0.40 1.00     7996     7623
## f0                    0.09      0.08    -0.07     0.24 1.00     8839     7855
## sex1                  0.05      0.09    -0.12     0.22 1.00     7963     8072
## f0:sex1               0.18      0.07     0.04     0.33 1.00    13956     8619
## sigma_f0             -0.08      0.06    -0.20     0.04 1.00     7125     7759
## sigma_sex1           -0.08      0.07    -0.21     0.05 1.00     6958     7416
## sigma_f0:sex1        -0.02      0.06    -0.14     0.10 1.00    12290    10347
## meageage_se          -0.07      0.05    -0.18     0.03 1.00    10929     8780
## sigma_meageage_se    -0.10      0.04    -0.17    -0.02 1.00     6713     6605
## 
## Samples were drawn using sample(hmc). 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).
kable(r2s, caption = "Explained variance, regular and leave-one-out.")
## Warning in `[<-.data.frame`(`*tmp*`, , isn, value = structure(list(Estimate =
## structure(c("0.1048387", : provided 8 variables to replace 5 variables
Explained variance, regular and leave-one-out.
model Estimate Est.Error Q2.5 Q97.5 loo_R2
h1_sex_mod_dataset_varying 0.1048387 0.0144906 0.0773468 0.1343458 0.08802218
h1_after_diagnosis_melsm 0.0970211 0.0140251 0.0704472 0.1252379 0.08525821
h1_after_diagnosis 0.0982942 0.0142484 0.0715919 0.1272165 0.08421256
h1_simple 0.0931076 0.0138371 0.0664255 0.1209109 0.08042350 
kable(coef(h1_sex_mod_dataset_varying)$dataset[, ,"f0"], caption = "Estimated association between f0 and outcome by dataset")
Estimated association between f0 and outcome by dataset
Estimate Est.Error Q2.5 Q97.5
2 0.1677698 0.1086444 -0.0354107 0.3929548
3 0.0088634 0.1120270 -0.2185637 0.2179185
4 0.0233515 0.1389023 -0.2866741 0.2746139
7 0.0618898 0.1398835 -0.2259794 0.3402301
8 0.0004170 0.1578933 -0.3732197 0.2647875
9 0.0171983 0.1254549 -0.2450564 0.2500621
10 0.1680425 0.1742130 -0.1147971 0.5806874

Hypothesis 4: SOI-R behavior

Participants with lower voice pitch will report having a more unrestricted sociosexual behavior.

Interpretation: The data are are consistent with substantial linear negative relationship between f0 and unrestricted sociosexuality. The 89% HDI for f0 falls entirely outside the ROPE, but the HDI for pf falls almost entirely within it. This evidence is consistent with a non-negligible association, where people with deeper voices have a less restricted sociosexuality (after adjusting for age and gender). Further research is needed to verify if the association with pf is truly negligible in size.

Visually diagnose non-linearity

ggplot(vcs, aes(f0, behavior)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 

ggplot(vcs, aes(pf, behavior)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 

ggplot(vcs, aes(age*10, behavior)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  scale_x_continuous("Age") +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

Decisions upon visual diagnosis: We will fit splines for f0, pf, and age. It seems likely that the nonlinearities in f0 and pf will not both be apparent after adjusting for one another. There is no evidence in the plots that there is an interaction with sex, except for age.

Models

h1_simple <- brm(behavior ~ f0 + pf + sex + age + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, file_refit = "on_change", save_mevars = TRUE, 
          file = "models/behavior/h1_simple") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Warning: Rows containing NAs were excluded from the model.
## Warning: Argument 'save_mevars' is deprecated. Please use argument 'latent' in
## function 'save_pars()' instead.
h1_after_diagnosis <- brm(behavior ~  s(f0) + s(pf) + sex + s(age, by = sex) + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, cores = chains, chains = chains,
          prior = priors, 
          control = control, file_refit = "on_change", save_mevars = TRUE,
          file = "models/behavior/h1_nonlinear") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Warning: Rows containing NAs were excluded from the model.

## Warning: Argument 'save_mevars' is deprecated. Please use argument 'latent' in
## function 'save_pars()' instead.
conditional_smooths(h1_after_diagnosis)

compare_models <- loo_compare(h1_simple, h1_after_diagnosis)
compare_models
##                    elpd_diff se_diff
## h1_after_diagnosis  0.0       0.0   
## h1_simple          -6.6       4.5
h1_sex_mod_dataset_varying <- brm(behavior ~ f0 + sex + s(age, by = sex) + 
                                    (1 + f0 + sex || dataset), data = vcs,
          prior = priors,
          iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
          control = control, file_refit = "on_change", save_mevars = TRUE,
          file = "models/behavior/h1_sex_mod_dataset_varying2") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Warning: Rows containing NAs were excluded from the model.

## Warning: Argument 'save_mevars' is deprecated. Please use argument 'latent' in
## function 'save_pars()' instead.
h1_after_diagnosis_melsm <- brm(
  bf(behavior ~  f0 + sex + s(age, by = sex) + (1 | dataset),
     sigma ~ f0 + sex + s(age, by = sex) + (1 | dataset)), data = vcs,
          prior = priors,
          iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
          control = control, file_refit = "on_change", save_mevars = TRUE,
          file = "models/behavior/h1_after_diagnosis_melsm") %>%
  add_criterion("loo") %>%
  add_criterion("bayes_R2") %>%
  add_criterion("loo_R2")
## Warning: Rows containing NAs were excluded from the model.

## Warning: Argument 'save_mevars' is deprecated. Please use argument 'latent' in
## function 'save_pars()' instead.
compare_models <- loo_compare(h1_simple, h1_after_diagnosis, h1_sex_mod_dataset_varying, h1_after_diagnosis_melsm)
compare_models
##                            elpd_diff se_diff
## h1_after_diagnosis_melsm      0.0       0.0 
## h1_sex_mod_dataset_varying -102.8      15.6 
## h1_after_diagnosis         -103.5      15.6 
## h1_simple                  -110.1      16.3
r2s <- list(h1_simple=h1_simple, h1_after_diagnosis=h1_after_diagnosis, h1_sex_mod_dataset_varying=h1_sex_mod_dataset_varying, h1_after_diagnosis_melsm=h1_after_diagnosis_melsm) %>% 
  map(~ { x <- as_tibble(bayes_R2(.))
          x$loo_R2 <- loo_R2(.)
          x
          }) %>%
  bind_rows(.id="model") %>% 
  arrange(desc(loo_R2))

sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
                  pred.labels = variable_labels)
  behavior
Predictors Estimates CI (95%)
Intercept -0.85 -1.20;-0.50
Voice pitch (f0) -0.21 -0.33;-0.09
Formant position (Pf) -0.02 -0.10;0.07
Sex [male] -0.24 -0.38;-0.09
Age 0.32 0.23;0.41
Random Effects
σ2 0.12
Ï„00 0.88
ICC 0.12
N dataset 9
Observations 1996
Marginal R2 / Conditional R2 0.173 / 0.173 
equiv_test <- equivalence_test(h1_simple, ci = 0.89, range = c(-0.1, 0.1))
equiv_test_table <- equiv_test %>% 
  mutate(Parameter = recode(Parameter, !!!effect_labels)) %>% 
  mutate(
    ROPE_Percentage = sprintf("%2.f%%", 100*ROPE_Percentage),
    HDI = sprintf("[% .2f;% .2f]", HDI_low, HDI_high)) %>% 
  select(Parameter, H0 = ROPE_Equivalence, `inside ROPE` = ROPE_Percentage, `89% HDI` = HDI)
kable(equiv_test_table, caption = "Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).")
Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).
Parameter H0 inside ROPE 89% HDI
Intercept Rejected 0% [-1.12;-0.56]
Voice pitch (f0) Rejected 0% [-0.31;-0.12]
Formant position (Pf) Accepted 100% [-0.09; 0.05]
Sex [male] Rejected 0% [-0.35;-0.12]
Age Rejected 0% [ 0.24; 0.39] 
plot(equiv_test) + scale_y_discrete("Effect", breaks = names(effect_labels), labels = effect_labels)
## Warning: Removed 3516 rows containing non-finite values (stat_density_ridges).

Model details 
summary(h1_simple)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: behavior ~ f0 + pf + sex + age + (1 | dataset) 
##    Data: vcs (Number of observations: 1996) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 9) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.41      0.14     0.23     0.75 1.00     1925     3771
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.85      0.18    -1.20    -0.50 1.00     2520     4240
## f0           -0.21      0.06    -0.33    -0.09 1.00     5719     5077
## pf           -0.02      0.04    -0.10     0.07 1.00     6518     5000
## sex1         -0.24      0.07    -0.38    -0.09 1.00     4768     5122
## age           0.32      0.05     0.23     0.41 1.00     6975     5146
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.91      0.01     0.89     0.94 1.00     7722     4831
## 
## Samples were drawn using sample(hmc). 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).
summary(h1_after_diagnosis)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: behavior ~ s(f0) + s(pf) + sex + s(age, by = sex) + (1 | dataset) 
##    Data: vcs (Number of observations: 1996) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Smooth Terms: 
##                      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sds(sf0_1)               0.74      0.73     0.03     2.68 1.00     2220
## sds(spf_1)               0.47      0.45     0.02     1.69 1.00     3481
## sds(sagesexfemale_1)     1.60      0.87     0.36     3.69 1.00     2776
## sds(sagesexmale_1)       1.20      0.80     0.21     3.21 1.00     2213
##                      Tail_ESS
## sds(sf0_1)               4085
## sds(spf_1)               3697
## sds(sagesexfemale_1)     2573
## sds(sagesexmale_1)       3376
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 9) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.40      0.14     0.22     0.77 1.00     1869     2986
## 
## Population-Level Effects: 
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept           -0.07      0.15    -0.37     0.22 1.01     1283     2173
## sex1                -0.24      0.09    -0.42    -0.05 1.00     6940     5539
## sf0_1               -1.49      1.37    -4.05     1.72 1.00     4147     3485
## spf_1               -0.15      1.09    -2.38     2.25 1.00     4332     3916
## sage:sexfemale_1     1.70      2.18    -3.11     5.64 1.00     4452     5120
## sage:sexmale_1       0.64      1.65    -2.99     3.65 1.00     4192     5371
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.91      0.01     0.88     0.94 1.00    11489     5753
## 
## Samples were drawn using sample(hmc). 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).
conditional_effects(h1_simple)

conditional_smooths(h1_after_diagnosis)

summary(h1_sex_mod_dataset_varying)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: behavior ~ f0 + sex + s(age, by = sex) + (1 + f0 + sex || dataset) 
##    Data: vcs (Number of observations: 1996) 
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup samples = 12000
## 
## Smooth Terms: 
##                      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sds(sagesexfemale_1)     1.56      0.86     0.35     3.62 1.00     5000
## sds(sagesexmale_1)       1.25      0.82     0.19     3.28 1.00     4133
##                      Tail_ESS
## sds(sagesexfemale_1)     4793
## sds(sagesexmale_1)       4422
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 9) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.40      0.14     0.20     0.74 1.00     5187     6496
## sd(f0)            0.09      0.08     0.00     0.29 1.00     3942     6609
## sd(sex1)          0.13      0.13     0.00     0.46 1.00     2993     5700
## 
## Population-Level Effects: 
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept           -0.07      0.15    -0.38     0.23 1.00     4620     6460
## f0                  -0.20      0.07    -0.35    -0.05 1.00    11443     7794
## sex1                -0.20      0.10    -0.39     0.01 1.00     8423     7107
## sage:sexfemale_1     1.71      2.13    -2.85     5.56 1.00     8680     8629
## sage:sexmale_1       0.65      1.69    -3.05     3.77 1.00     9322     8428
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.91      0.01     0.88     0.94 1.00    21421     7927
## 
## Samples were drawn using sample(hmc). 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).
summary(h1_after_diagnosis_melsm)
##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: behavior ~ f0 + sex + s(age, by = sex) + (1 | dataset) 
##          sigma ~ f0 + sex + s(age, by = sex) + (1 | dataset)
##    Data: vcs (Number of observations: 1996) 
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup samples = 12000
## 
## Smooth Terms: 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sds(sagesexfemale_1)           1.59      0.85     0.37     3.62 1.00     5139
## sds(sagesexmale_1)             1.04      0.70     0.17     2.82 1.00     4535
## sds(sigma_sagesexfemale_1)     0.57      0.52     0.02     1.90 1.00     3743
## sds(sigma_sagesexmale_1)       0.78      0.55     0.10     2.22 1.00     4437
##                            Tail_ESS
## sds(sagesexfemale_1)           5000
## sds(sagesexmale_1)             4483
## sds(sigma_sagesexfemale_1)     5544
## sds(sigma_sagesexmale_1)       4582
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 9) 
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)           0.41      0.14     0.23     0.75 1.00     3749     5127
## sd(sigma_Intercept)     0.30      0.10     0.17     0.56 1.00     3670     6236
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                 -0.05      0.15    -0.34     0.25 1.00     2878
## sigma_Intercept           -0.19      0.11    -0.41     0.02 1.00     2977
## f0                        -0.17      0.05    -0.27    -0.07 1.00    10941
## sex1                      -0.17      0.06    -0.28    -0.06 1.00    10214
## sigma_f0                  -0.01      0.04    -0.09     0.08 1.00     9978
## sigma_sex1                 0.06      0.05    -0.04     0.15 1.00     9138
## sage:sexfemale_1           1.73      2.14    -2.93     5.67 1.00     7697
## sage:sexmale_1             0.58      1.59    -2.91     3.40 1.00     6873
## sigma_sage:sexfemale_1     0.19      1.55    -3.27     3.29 1.00     5882
## sigma_sage:sexmale_1       0.47      1.52    -2.41     3.87 1.00     6153
##                        Tail_ESS
## Intercept                  3876
## sigma_Intercept            4384
## f0                         9219
## sex1                       8282
## sigma_f0                   8645
## sigma_sex1                 8671
## sage:sexfemale_1           8618
## sage:sexmale_1             8213
## sigma_sage:sexfemale_1     4636
## sigma_sage:sexmale_1       6857
## 
## Samples were drawn using sample(hmc). 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).
kable(r2s, caption = "Explained variance, regular and leave-one-out.")
## Warning in `[<-.data.frame`(`*tmp*`, , isn, value = structure(list(Estimate =
## structure(c("0.1785773", : provided 8 variables to replace 5 variables
Explained variance, regular and leave-one-out.
model Estimate Est.Error Q2.5 Q97.5 loo_R2
h1_after_diagnosis_melsm 0.1785773 0.0139789 0.1513007 0.2063147 0.17122084
h1_sex_mod_dataset_varying 0.1841389 0.0142630 0.1559592 0.2118939 0.16994368
h1_after_diagnosis 0.1840598 0.0142224 0.1569562 0.2124536 0.16927791
h1_simple 0.1729142 0.0138575 0.1460109 0.2003496 0.16441387 
kable(coef(h1_sex_mod_dataset_varying)$dataset[, ,"f0"], caption = "Estimated association between f0 and outcome by dataset")
Estimated association between f0 and outcome by dataset
Estimate Est.Error Q2.5 Q97.5
1 -0.2102720 0.0800788 -0.3710549 -0.0538281
2 -0.2581127 0.0983075 -0.4915334 -0.0959797
3 -0.1569973 0.0976920 -0.3212690 0.0642548
4 -0.1488008 0.1106081 -0.3282048 0.1130664
5 -0.2076882 0.0936743 -0.3958332 -0.0188280
6 -0.2002083 0.1028056 -0.4058073 0.0125026
7 -0.2441958 0.1133408 -0.5125175 -0.0492377
8 -0.1814103 0.1089961 -0.3833508 0.0665820
11 -0.2281538 0.1284119 -0.5205214 0.0176658 
equivalence_test(h1_after_diagnosis_melsm, ci = 0.89, range = c(-0.1, 0.1))
## # Test for Practical Equivalence
## 
##   ROPE: [-0.10 0.10]
## 
## Parameter              |        H0 | inside ROPE |       89% HDI
## ----------------------------------------------------------------
## Intercept              | Undecided |     59.26 % | [-0.26  0.18]
## sigma_Intercept        | Undecided |     12.12 % | [-0.36 -0.04]
## f0                     | Undecided |      4.57 % | [-0.25 -0.08]
## sex1                   | Undecided |      5.61 % | [-0.26 -0.08]
## sigma_f0               |  Accepted |    100.00 % | [-0.08  0.06]
## sigma_sex1             | Undecided |     85.67 % | [-0.02  0.13]
## sage.sexfemale_1       | Undecided |      2.52 % | [-1.69  5.07]
## sage.sexmale_1         | Undecided |      4.57 % | [-1.96  2.99]
## sigma_sage.sexfemale_1 | Undecided |      7.97 % | [-2.10  2.43]
## sigma_sage.sexmale_1   | Undecided |      8.00 % | [-1.81  2.88]

Summary of results

Models

dominance <- readRDS("models/dominance/h1_simple.rds")
agree <- readRDS("models/agree/h1_simple2.rds")
extra <- readRDS("models/extra/h1_simple.rds")
soi <- readRDS("models/behavior/h1_simple.rds")
theme_set(theme_cowplot())

modify <- function(outcome) {
  list(scale_y_discrete(outcome, breaks = names(effect_labels), labels = effect_labels),
  guides(fill = "none"),
  theme(panel.spacing = unit(c(0), "cm"),
        plot.margin = unit(c(0, 0.5, 0, 0.5), "cm"))
)
}
remove_x_axis <- list(
  scale_x_continuous("", limits = c(-0.5, 0.4)),
  theme(axis.ticks.x = element_blank(),
        axis.text.x = element_blank(),
        axis.line.x = element_blank())
)

p_dominance <- plot(equivalence_test(dominance, parameters = c("b_f0", "b_pf", "b_sex1"), ci = 0.89)) + modify("Dominance") + remove_x_axis
p_agree <- plot(equivalence_test(agree, parameters = c("b_f0", "b_pf", "b_sex1"), ci = 0.89)) +
  modify("Agreeableness") + remove_x_axis
p_extra <- plot(equivalence_test(extra, parameters = c("b_f0", "b_pf", "b_sex1"), ci = 0.89))  +
  modify("Extraversion") + remove_x_axis
p_soi <- plot(equivalence_test(soi, parameters = c("b_f0", "b_pf", "b_sex1"), ci = 0.89)) +
  modify("SOI-R behavior") + xlim(c(-0.5, 0.4))

soir_full <- readRDS("models/soir_full/h1_simple.rds")
p_soir_full <- plot(equivalence_test(soir_full, parameters = c("b_f0", "b_pf", "b_sex1"), ci = 0.89)) +
  modify("Sociosexuality") + xlim(c(-0.5, 0.4))

legend <- get_legend(
  p_soir_full + 
    guides(fill = guide_legend(nrow = 1)) +
    theme(legend.position = "bottom", legend.justification = "center")
)
## Warning: Removed 2637 rows containing non-finite values (stat_density_ridges).
neg_spacer <- -0.17
aligned_plots <- align_plots(p_dominance, p_extra, p_agree, p_soi, align="hv")
## Warning: Removed 2637 rows containing non-finite values (stat_density_ridges).

## Warning: Removed 2637 rows containing non-finite values (stat_density_ridges).

## Warning: Removed 2637 rows containing non-finite values (stat_density_ridges).

## Warning: Removed 2637 rows containing non-finite values (stat_density_ridges).
plot_grid(ncol = 1,  rel_heights = c(1, neg_spacer, 0.8, neg_spacer, 1, neg_spacer, 1.1, .15),
          aligned_plots[[1]], NULL, aligned_plots[[2]], NULL, aligned_plots[[3]],NULL, aligned_plots[[4]], legend)

ggsave("preregistered_rope.pdf", width = 8, height = 10)

Tabular

get_coefs <- function(model) {
model_summary <- summary(model)
nonvarying <- model_summary$fixed %>% 
  as.data.frame() %>% 
  rownames_to_column("term") %>% 
  mutate(group = "non-varying",
         outcome = as.character(model_summary$formula[[1]][[2]]))
varying <- model_summary$random %>% 
  map(~ rownames_to_column(as.data.frame(.), "term")) %>% 
  bind_rows(.id = "group") %>% 
  left_join(model_summary$ngrps %>% 
    as_tibble() %>% 
    gather(group, n), by = "group") %>% 
  mutate(#group = paste0(group, " (n=", n, ")"),
         outcome = as.character(model_summary$formula[[1]][[2]]))

two_digits <- function(x) { sprintf("%.2f", x) }
coefs <- bind_rows(nonvarying, 
                   varying)
coefs <- coefs %>% select(group, outcome, term, Estimate, `l-95% CI`, `u-95% CI`)
coefs$Estimate = two_digits(coefs$Estimate)
coefs$`l-95% CI` = two_digits(coefs$`l-95% CI`)
coefs$`u-95% CI` = two_digits(coefs$`u-95% CI`)
coefs
}

coefs <- bind_rows(get_coefs(dominance), 
                   get_coefs(extra), get_coefs(agree), get_coefs(soi))
## Warning: There were 3 divergent transitions after warmup. Increasing adapt_delta
## above may help. See http://mc-stan.org/misc/warnings.html#divergent-transitions-
## after-warmup
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.

## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
wide_coefs <- coefs %>% 
  mutate(Estimate = paste0(Estimate, " [", `l-95% CI`, ";", `u-95% CI`, "]")) %>% 
    select(group, outcome, term, Estimate) %>% 
  mutate(term = recode(term, "meageage_se" = "Age(±SE)", "age" = "Age(±SE)"),
         outcome = fct_inorder(outcome),
         term = fct_inorder(term)) %>% 
  spread(outcome, Estimate) %>% 
  arrange(desc(group), term) %>% 
  mutate(group = fct_inorder(group))

outcome_names <- c("Extraversion" = "extra", "Dominance" = "dominance",
                   "Agreeableness" = "agree", "SOI-R behavior" = "behavior")

library(kableExtra)
wide_coefs %>% 
  mutate(term = str_replace_all(term, variable_labels)) %>% 
  rename(!!!syms(outcome_names), Term = term) %>% 
  select(-group) %>% 
  kable(caption = "Exploratory analysis results") %>%
  pack_rows(index = table(wide_coefs$group)) %>% 
  column_spec(column = 1, width = "3.3cm") %>% 
  column_spec(column = 2:(ncol(wide_coefs)-1), width = "3cm") %>% 
  add_header_above(c(" " = 1, "Estimated effect on each outcome [95% CI]" = ncol(wide_coefs)-2)) %>% 
  footnote("Estimated associations between outcomes, for which we had made preregistered predictions.")
Exploratory analysis results
Estimated effect on each outcome [95% CI]
Term Dominance Extraversion Agreeableness SOI-R behavior
non-varying
Intercept -0.07 [-0.80;0.68] 0.19 [-0.21;0.58] 0.09 [-0.31;0.50] -0.85 [-1.20;-0.50]
Voice pitch (f0) -0.27 [-0.45;-0.08] -0.23 [-0.38;-0.08] 0.03 [-0.12;0.18] -0.21 [-0.33;-0.09]
Formant position (Pf) -0.00 [-0.13;0.12] 0.08 [-0.03;0.18] 0.03 [-0.07;0.13] -0.02 [-0.10;0.07]
Sex [male] -0.30 [-0.52;-0.08] -0.28 [-0.46;-0.10] -0.00 [-0.18;0.18] -0.24 [-0.38;-0.09]
Age(±SE) 0.04 [-0.08;0.17] -0.11 [-0.22;0.01] -0.06 [-0.17;0.05] 0.32 [0.23;0.41]
dataset
sd(Intercept) 0.53 [0.16;1.68] 0.35 [0.16;0.73] 0.37 [0.19;0.78] 0.41 [0.23;0.75]
Note:
Estimated associations between outcomes, for which we had made preregistered predictions.

Scatterplots

dominance <- ggplot(vcs, aes(f0, dominance)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(aes(colour = sex, fill = sex), method = "lm") +
  scale_color_manual("Sex", values = c("#40A5BF", "black")) +
  scale_fill_manual("Sex", values = c("#40A5BF", "black"))  +
  guides(color = F, fill = F) +
  xlab("") +
  ylab("Dominance") +
  theme(panel.spacing = unit(c(0), "cm"),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))

extra <- ggplot(vcs, aes(f0, extra)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(aes(colour = sex, fill = sex), method = "lm") +
  scale_color_manual("Sex", values = c("#40A5BF", "black")) +
  scale_fill_manual("Sex", values = c("#40A5BF", "black"))  +
  guides(color = F, fill = F) +
  xlab("") +
  ylab("Extraversion") +
  theme(panel.spacing = unit(c(0), "cm"),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))

agree <- ggplot(vcs, aes(f0, agree)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(aes(colour = sex, fill = sex), method = "lm") +
  scale_color_manual("Sex", values = c("#40A5BF", "black")) +
  scale_fill_manual("Sex", values = c("#40A5BF", "black"))  +
  guides(color = F, fill = F) +
  xlab("Voice pitch (f0)") +
  ylab("Agreeableness") +
  theme(panel.spacing = unit(c(0), "cm"),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))


behavior <- ggplot(vcs, aes(f0, behavior)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(aes(colour = sex, fill = sex), method = "lm") +
  scale_color_manual("Sex", values = c("#40A5BF", "black")) +
  scale_fill_manual("Sex", values = c("#40A5BF", "black"))  +
  guides(color = F, fill = F) +
  xlab("Voice pitch (f0)") +
  ylab("SOI-R Behavior") +
  theme(panel.spacing = unit(c(0), "cm"),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))


plot_grid(ncol = 2,  dominance, extra, agree, behavior, rel_heights = c(1, 1, 1, 1),
          align = "hv")
## Warning: Removed 1245 rows containing non-finite values (stat_smooth).
## Warning: Removed 1245 rows containing missing values (geom_point).
## Warning: Removed 797 rows containing non-finite values (stat_smooth).
## Warning: Removed 797 rows containing missing values (geom_point).
## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing missing values (geom_point).
## Warning: Removed 234 rows containing non-finite values (stat_smooth).
## Warning: Removed 234 rows containing missing values (geom_point).

ggsave("scatterplots_prereg.pdf", width = 8, height = 8)