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

Missingness patterns

codebook::md_pattern(vcs %>% select(dominance, extra, behavior, sex, age, f0, pf, dataset))
## # A tibble: 6 x 6
##   description                     behavior extra dominance var_miss n_miss
##   <chr>                              <dbl> <dbl>     <dbl>    <dbl>  <dbl>
## 1 Missing values per variable          234   797      1245     2276   2276
## 2 Missing values in 0 variables          1     1         1        0    969
## 3 Missing values in 2 variables          1     0         0        2    760
## 4 Missing values in 1 variables          1     1         0        1    267
## 5 Missing values in 2 variables          0     1         0        2    196
## 6 3 other, less frequent patterns        0     1         2        6     38
vcs <- vcs %>% filter(!is.na(f0), !is.na(pf), !is.na(age))

Neuroticism

Visually diagnose non-linearity

ggplot(vcs, aes(f0, neuro)) + 
  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, neuro)) + 
  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, neuro)) + 
  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:

Models

h0 <- brm(neuro ~ 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/neuro/h0") %>% 
  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_simple <- brm(neuro ~ 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/neuro/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.
compare_models <- loo_compare(h0, h1_simple)
compare_models
##           elpd_diff se_diff
## h1_simple  0.0       0.0   
## h0        -0.4       2.4
sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
                  pred.labels = variable_labels)
  neuro
Predictors Estimates CI (95%)
Intercept 0.04 -0.32;0.37
Voice pitch (f0) 0.16 0.02;0.31
Formant position (Pf) -0.05 -0.15;0.05
Sex [male] -0.15 -0.32;0.03
Age -0.01 -0.12;0.10
Random Effects
σ2 0.03
Ï„00 0.97
ICC 0.03
N dataset 7
Observations 1434
Marginal R2 / Conditional R2 0.111 / 0.111
Model details 
summary(h1_simple)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: neuro ~ f0 + pf + sex + age + (1 | dataset) 
##    Data: vcs (Number of observations: 1434) 
## 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.26      0.12     0.11     0.56 1.00     1962     3227
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.04      0.18    -0.32     0.37 1.00     3845     4355
## f0            0.16      0.07     0.02     0.31 1.00     5247     5331
## pf           -0.05      0.05    -0.15     0.05 1.00     5938     5342
## sex1         -0.15      0.09    -0.32     0.03 1.00     4447     5117
## age          -0.01      0.06    -0.12     0.10 1.00     6680     5021
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.94      0.02     0.91     0.98 1.00     7332     5501
## 
## 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)
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 49% [-0.23; 0.33]
Voice pitch (f0) Undecided 15% [ 0.05; 0.29]
Formant position (Pf) Undecided 86% [-0.14; 0.02]
Sex [male] Undecided 27% [-0.29;-0.00]
Age Accepted 100% [-0.10; 0.08] 
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).

Openness

Visually diagnose non-linearity

ggplot(vcs, aes(f0, openn)) + 
  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, openn)) + 
  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, openn)) + 
  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:

Models

h0 <- brm(openn ~ 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/openn/h0") %>% 
  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_simple <- brm(openn ~ 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/openn/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.
compare_models <- loo_compare(h0, h1_simple)
compare_models
##           elpd_diff se_diff
## h1_simple  0.0       0.0   
## h0        -0.6       2.1
sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
                  pred.labels = variable_labels)
  openn
Predictors Estimates CI (95%)
Intercept -0.12 -0.41;0.19
Voice pitch (f0) -0.17 -0.32;-0.01
Formant position (Pf) 0.00 -0.10;0.11
Sex [male] -0.22 -0.42;-0.04
Age 0.03 -0.08;0.14
Random Effects
σ2 0.01
Ï„00 1.00
ICC 0.01
N dataset 7
Observations 1434
Marginal R2 / Conditional R2 0.017 / 0.017
Model details 
summary(h1_simple)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: openn ~ f0 + pf + sex + age + (1 | dataset) 
##    Data: vcs (Number of observations: 1434) 
## 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.15      0.09     0.03     0.38 1.00     1634     2148
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.12      0.15    -0.41     0.19 1.00     5041     4951
## f0           -0.17      0.08    -0.32    -0.01 1.00     5115     5349
## pf            0.00      0.05    -0.10     0.11 1.00     4402     4960
## sex1         -0.22      0.10    -0.42    -0.04 1.00     3290     4546
## age           0.03      0.05    -0.08     0.14 1.00     6400     5901
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.99      0.02     0.96     1.03 1.00     7234     5030
## 
## 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)
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 41% [-0.36; 0.13]
Voice pitch (f0) Undecided 14% [-0.30;-0.05]
Formant position (Pf) Accepted 100% [-0.09; 0.09]
Sex [male] Undecided 5% [-0.38;-0.07]
Age Undecided 95% [-0.06; 0.12] 
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).

Conscientiousness

Visually diagnose non-linearity

ggplot(vcs, aes(f0, consc)) + 
  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, consc)) + 
  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, consc)) + 
  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:

Models

h0 <- brm(consc ~ 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/consc/h0") %>% 
  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_simple <- brm(consc ~ 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/consc/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.
compare_models <- loo_compare(h0, h1_simple)
compare_models
##           elpd_diff se_diff
## h0         0.0       0.0   
## h1_simple -1.8       0.8
sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
                  pred.labels = variable_labels)
  consc
Predictors Estimates CI (95%)
Intercept -0.14 -0.50;0.25
Voice pitch (f0) -0.02 -0.17;0.13
Formant position (Pf) 0.04 -0.06;0.14
Sex [male] -0.09 -0.27;0.09
Age 0.01 -0.10;0.12
Random Effects
σ2 0.08
Ï„00 0.93
ICC 0.08
N dataset 7
Observations 1434
Marginal R2 / Conditional R2 0.091 / 0.091
Model details 
summary(h1_simple)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: consc ~ f0 + pf + sex + age + (1 | dataset) 
##    Data: vcs (Number of observations: 1434) 
## 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.32      0.14     0.15     0.67 1.00     2286     3266
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.14      0.19    -0.50     0.25 1.00     3366     4156
## f0           -0.02      0.08    -0.17     0.13 1.00     4931     4417
## pf            0.04      0.05    -0.06     0.14 1.00     5104     5500
## sex1         -0.09      0.09    -0.27     0.09 1.00     4116     4700
## age           0.01      0.06    -0.10     0.12 1.00     5986     5357
## 
## 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     7254     5090
## 
## 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)
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 34% [-0.42; 0.18]
Voice pitch (f0) Undecided 90% [-0.13; 0.11]
Formant position (Pf) Undecided 94% [-0.04; 0.12]
Sex [male] Undecided 53% [-0.24; 0.05]
Age Undecided 100% [-0.08; 0.10] 
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).

SOI-R Desire

Visually diagnose non-linearity

ggplot(vcs, aes(f0, desire)) + 
  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, desire)) + 
  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, desire)) + 
  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 show 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

h0 <- brm(desire ~ 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/desire/h0") %>% 
  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_simple <- brm(desire ~ 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/desire/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.
compare_models <- loo_compare(h0, h1_simple)
compare_models
##           elpd_diff se_diff
## h1_simple  0.0       0.0   
## h0        -5.1       3.6
sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
                  pred.labels = variable_labels)
  desire
Predictors Estimates CI (95%)
Intercept 0.31 -0.05;0.66
Voice pitch (f0) -0.22 -0.33;-0.10
Formant position (Pf) 0.01 -0.07;0.10
Sex [male] 0.23 0.09;0.37
Age -0.12 -0.21;-0.04
Random Effects
σ2 0.05
Ï„00 0.96
ICC 0.05
N dataset 9
Observations 2000
Marginal R2 / Conditional R2 0.217 / 0.217
Model details 
summary(h1_simple)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: desire ~ f0 + pf + sex + age + (1 | dataset) 
##    Data: vcs (Number of observations: 2000) 
## 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.42      0.14     0.24     0.75 1.00     1940     3460
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.31      0.18    -0.05     0.66 1.00     2793     3569
## f0           -0.22      0.06    -0.33    -0.10 1.00     5752     5038
## pf            0.01      0.04    -0.07     0.10 1.00     6768     5540
## sex1          0.23      0.07     0.09     0.37 1.00     5338     5154
## age          -0.12      0.05    -0.21    -0.04 1.00     7426     5354
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.89      0.01     0.86     0.92 1.00     7909     5360
## 
## 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)
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 7% [ 0.02; 0.59]
Voice pitch (f0) Rejected 0% [-0.31;-0.13]
Formant position (Pf) Accepted 100% [-0.06; 0.08]
Sex [male] Rejected 0% [ 0.11; 0.34]
Age Undecided 27% [-0.20;-0.05] 
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).

SOI-R Attitude

Visually diagnose non-linearity

ggplot(vcs, aes(f0, attitude)) + 
  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, attitude)) + 
  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, attitude)) + 
  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 show 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

h0 <- brm(attitude ~ 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/attitude/h0") %>% 
  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_simple <- brm(attitude ~ 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/attitude/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.
compare_models <- loo_compare(h0, h1_simple)
compare_models
##           elpd_diff se_diff
## h1_simple  0.0       0.0   
## h0        -7.5       4.5
sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
                  pred.labels = variable_labels)
  attitude
Predictors Estimates CI (95%)
Intercept -0.02 -0.45;0.41
Voice pitch (f0) -0.24 -0.35;-0.12
Formant position (Pf) -0.05 -0.14;0.03
Sex [male] -0.07 -0.21;0.07
Age -0.00 -0.09;0.09
Random Effects
σ2 0.15
Ï„00 0.85
ICC 0.15
N dataset 9
Observations 1999
Marginal R2 / Conditional R2 0.199 / 0.199
Model details 
summary(h1_simple)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: attitude ~ f0 + pf + sex + age + (1 | dataset) 
##    Data: vcs (Number of observations: 1999) 
## 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.54      0.17     0.32     0.97 1.00     1999     3478
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.02      0.22    -0.45     0.41 1.00     2313     3203
## f0           -0.24      0.06    -0.35    -0.12 1.00     5468     5182
## pf           -0.05      0.04    -0.14     0.03 1.00     6146     5099
## sex1         -0.07      0.07    -0.21     0.07 1.00     4659     5209
## age          -0.00      0.04    -0.09     0.09 1.00     7365     5977
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.90      0.01     0.87     0.93 1.00     7885     5516
## 
## 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)
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 42% [-0.36; 0.33]
Voice pitch (f0) Rejected 0% [-0.33;-0.15]
Formant position (Pf) Undecided 91% [-0.12; 0.02]
Sex [male] Undecided 68% [-0.18; 0.05]
Age Accepted 100% [-0.07; 0.07] 
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).

Sociosexuality

Visually diagnose non-linearity

ggplot(vcs, aes(f0, soir_full)) + 
  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, soir_full)) + 
  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, soir_full)) + 
  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 show 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

h0 <- brm(soir_full ~ 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/soir_full/h0") %>% 
  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_simple <- brm(soir_full ~ 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/soir_full/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.
compare_models <- loo_compare(h0, h1_simple)
compare_models
##           elpd_diff se_diff
## h1_simple   0.0       0.0  
## h0        -10.7       5.2
sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
                  pred.labels = variable_labels)
  soir_full
Predictors Estimates CI (95%)
Intercept -0.23 -0.66;0.19
Voice pitch (f0) -0.28 -0.39;-0.17
Formant position (Pf) -0.03 -0.11;0.06
Sex [male] -0.04 -0.18;0.10
Age 0.08 -0.01;0.17
Random Effects
σ2 0.15
Ï„00 0.86
ICC 0.15
N dataset 9
Observations 1995
Marginal R2 / Conditional R2 0.228 / 0.228
Model details 
summary(h1_simple)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: soir_full ~ f0 + pf + sex + age + (1 | dataset) 
##    Data: vcs (Number of observations: 1995) 
## 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.54      0.17     0.32     0.98 1.00     1929     3158
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.23      0.22    -0.66     0.19 1.00     2151     3008
## f0           -0.28      0.06    -0.39    -0.17 1.00     5663     5491
## pf           -0.03      0.04    -0.11     0.06 1.00     6390     5260
## sex1         -0.04      0.07    -0.18     0.10 1.00     5236     4933
## age           0.08      0.04    -0.01     0.17 1.00     7140     5153
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.88      0.01     0.85     0.91 1.00     7794     5282
## 
## 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)
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 23% [-0.57; 0.11]
Voice pitch (f0) Rejected 0% [-0.38;-0.20]
Formant position (Pf) Accepted 100% [-0.09; 0.04]
Sex [male] Undecided 87% [-0.15; 0.08]
Age Undecided 68% [ 0.01; 0.15] 
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).

Summary of results

Models

openn <- readRDS("models/openn/h1_simple.rds")
neuro <- readRDS("models/neuro/h1_simple.rds")
consc <- readRDS("models/consc/h1_simple.rds")
desire <- readRDS("models/desire/h1_simple.rds")
attitude <- readRDS("models/attitude/h1_simple.rds")
soir_full <- readRDS("models/soir_full/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_extra <- plot(equivalence_test(neuro, parameters = c("b_f0", "b_pf", "b_sex1"), ci = 0.89)) + modify("Neuroticism") + remove_x_axis
p_openn <- plot(equivalence_test(openn, parameters = c("b_f0", "b_pf", "b_sex1"), ci = 0.89)) +
  modify("Openness") + remove_x_axis
p_consc <- plot(equivalence_test(consc, parameters = c("b_f0", "b_pf", "b_sex1"), ci = 0.89))  +
  modify("Conscientiousness") + remove_x_axis
p_attitude <- plot(equivalence_test(attitude, parameters = c("b_f0", "b_pf", "b_sex1"), ci = 0.89)) +
  modify("SOI-R attitude") + remove_x_axis
p_desire <- plot(equivalence_test(desire, parameters = c("b_f0", "b_pf", "b_sex1"), ci = 0.89)) +
  modify("SOI-R desire") + remove_x_axis
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_extra + 
    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.155
plot_grid(ncol = 1,  rel_heights = c(1, neg_spacer, 1, neg_spacer, 1,neg_spacer, 1,neg_spacer, 1, neg_spacer, 1.1, .15), align = "v", axis = "b",
          p_extra, NULL, p_openn,NULL, p_consc,NULL, p_attitude, NULL, p_desire, NULL, p_soir_full, legend)
## 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).

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

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

ggsave("exploratory_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(neuro), 
                   get_coefs(openn), get_coefs(consc), get_coefs(desire),
                   get_coefs(attitude),get_coefs(soir_full))

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("Neuroticism" = "neuro", "Conscientiousness" = "consc",
                   "Openness" = "openn", "Sociosexuality" = "soir_full",
                   "SOI-R desire" = "desire", "SOI-R attitude" = "attitude")

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 no predictions.")
Exploratory analysis results
Estimated effect on each outcome [95% CI]
Term Neuroticism Openness Conscientiousness SOI-R desire SOI-R attitude Sociosexuality
non-varying
Intercept 0.04 [-0.32;0.37] -0.12 [-0.41;0.19] -0.14 [-0.50;0.25] 0.31 [-0.05;0.66] -0.02 [-0.45;0.41] -0.23 [-0.66;0.19]
Voice pitch (f0) 0.16 [0.02;0.31] -0.17 [-0.32;-0.01] -0.02 [-0.17;0.13] -0.22 [-0.33;-0.10] -0.24 [-0.35;-0.12] -0.28 [-0.39;-0.17]
Formant position (Pf) -0.05 [-0.15;0.05] 0.00 [-0.10;0.11] 0.04 [-0.06;0.14] 0.01 [-0.07;0.10] -0.05 [-0.14;0.03] -0.03 [-0.11;0.06]
Sex [male] -0.15 [-0.32;0.03] -0.22 [-0.42;-0.04] -0.09 [-0.27;0.09] 0.23 [0.09;0.37] -0.07 [-0.21;0.07] -0.04 [-0.18;0.10]
Age(±SE) -0.01 [-0.12;0.10] 0.03 [-0.08;0.14] 0.01 [-0.10;0.12] -0.12 [-0.21;-0.04] -0.00 [-0.09;0.09] 0.08 [-0.01;0.17]
dataset
sd(Intercept) 0.26 [0.11;0.56] 0.15 [0.03;0.38] 0.32 [0.15;0.67] 0.42 [0.24;0.75] 0.54 [0.32;0.97] 0.54 [0.32;0.98]
Note:
Estimated associations between outcomes, for which we had made no predictions.

Scatterplots

outcome_names <- c("Neuroticism" = "neuro", "Conscientiousness" = "consc",
                   "Openness" = "openn", "Sociosexuality" = "soir_full",
                   "SOI-R desire" = "desire", "SOI-R attitude" = "attitude")

neuro <- ggplot(vcs, aes(f0, neuro)) + 
  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("Neuroticism") +
  theme(panel.spacing = unit(c(0), "cm"),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))

consc <- ggplot(vcs, aes(f0, consc)) + 
  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("Conscientiousness") +
  theme(panel.spacing = unit(c(0), "cm"),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))

openn <- ggplot(vcs, aes(f0, openn)) + 
  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("Openness") +
  theme(panel.spacing = unit(c(0), "cm"),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))


soir_full <- ggplot(vcs, aes(f0, soir_full)) + 
  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("Sociosexuality") +
  theme(panel.spacing = unit(c(0), "cm"),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))


desire <- ggplot(vcs, aes(f0, desire)) + 
  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 Desire") +
  theme(panel.spacing = unit(c(0), "cm"),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))


attitude <- ggplot(vcs, aes(f0, attitude)) + 
  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 Attitude") +
  theme(panel.spacing = unit(c(0), "cm"),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))

plot_grid(ncol = 2,  openn, neuro, consc, soir_full, attitude, desire,
          align = "hv")
## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 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 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing missing values (geom_point).
## Warning: Removed 235 rows containing non-finite values (stat_smooth).
## Warning: Removed 235 rows containing missing values (geom_point).
## Warning: Removed 231 rows containing non-finite values (stat_smooth).
## Warning: Removed 231 rows containing missing values (geom_point).
## Warning: Removed 230 rows containing non-finite values (stat_smooth).
## Warning: Removed 230 rows containing missing values (geom_point).

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