Inbreeding aversion analyses

Pre-registered hypothesis:

  1. Single women spend less time with male biological relatives during the FP (for paired women, we did not include questions about relative sex at this level of detail for reasons of length, but we would not make a different prediction).

  2. All women spend less time with relatives during the FP (this question was asked for all participating women).

load packages

source("0_helpers.R")
library(knitr) # um diese dokumente zu generieren
library(formr) # formr paket, kann daten laden
opts_chunk$set(warning = F, message = F, error = T, fig.width = 7, fig.height = 6)
library(lubridate) # mit datum-uhrzeit rechnen
library(pander) # schoen formatierte tabellen
library(ggplot2)
library(tidyr) # daten von lang nach breit und zurueck
library(dplyr) # daten zusammenfassen etc. immer als letztes laden
diary <- readRDS("diary_persons.rds")
library(brms)

time family (singles)

library(brms)
library(bayesplot)
model2_family <- brm(time_family ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person), data = diary %>% filter(reasons_for_exclusion == "" | reasons_for_exclusion == "not_single, "), cores = 4, backend = "cmdstanr",file = "brms_fits/model2_family3")
hypothesis(model2_family, "fertile_fab < 0", alpha = 0.01)
## Hypothesis Tests for class b:
##          Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
## 1 (fertile_fab) < 0    -0.01      0.06    -0.15     0.13       1.26      0.56     
## ---
## 'CI': 98%-CI for one-sided and 99%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 99%;
## for two-sided hypotheses, the value tested against lies outside the 99%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
hypothesis(model2_family, "fertile_fab:hormonal_contraceptionTRUE > 0", alpha = 0.01)
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
## 1 (fertile_fab:horm... > 0     0.07      0.12     -0.2     0.34       2.67      0.73     
## ---
## 'CI': 98%-CI for one-sided and 99%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 99%;
## for two-sided hypotheses, the value tested against lies outside the 99%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
model2_family_curve <- brm(time_family ~ s(RCD_fab,  by = hormonal_contraception, bs = "cc") + (1 | person), data =  diary %>% filter(reasons_for_exclusion == "" | reasons_for_exclusion == "not_single, "), backend = "cmdstanr", cores = 4, file ="brms_fits/model3_family_curve")

curve_plot(model2_family_curve, "Time spent with family")
library(brms)
outcome_var <- model$formula[[1]][[2]]

smooths30 <- conditional_effects(model, 
      effects = "RCD_fab:hormonal_contraception",
      spaghetti = TRUE, resolution = 60, 
      nsamples = 100, dpars = "mu",
      re_formula = NA)

library(tidyverse)

smo_data <- attributes(smooths30$`RCD_fab:hormonal_contraception`)$spaghetti %>% 
  mutate(nohc = if_else(hormonal_contraception == TRUE, "0", "1")) %>% 
  arrange(nohc)

smo_diff <- full_join(
  smo_data %>% 
    filter(nohc == "1") %>%
    mutate(sample__ = str_split_fixed(sample__, "_", n = 2)[,1]),
  smo_data %>% 
    filter(nohc == "0") %>% 
    mutate(sample__ = str_split_fixed(sample__, "_", n = 2)[,1]),
  by = c("RCD_fab", "sample__"), 
  suffix = c("_nohc", "_hc")
  ) %>% 
    mutate(estimate__ = estimate___nohc - estimate___hc,
           nohc = "diff")

smo_data <- bind_rows(smo_data, smo_diff)

if(model$family$family %in% c("gaussian", "bernoulli")) {
  sd_fun <- sd
} else if (model$family$family == "poisson") {
  sd_fun = mean
} 
y_axes <- bind_rows(
  model$data %>% 
    mutate(nohc = if_else(hormonal_contraception == TRUE, "0", "1")) %>% 
    group_by(nohc) %>% 
    summarise(
      n_women = n_distinct(person),
      outcomemeanplus1sd = min(mean(!!outcome_var, na.rm = T) + sd_fun(!!outcome_var, na.rm = T), max(!!outcome_var, na.rm = T)),
      outcomemeanminus1sd = max(mean(!!outcome_var, na.rm = T) - sd_fun(!!outcome_var, na.rm = T), min(!!outcome_var, na.rm = T))),
  model$data %>% 
    summarise(
      outcomemeanplus1sd = 0 + sd_fun(!!outcome_var, na.rm = T),
      outcomemeanminus1sd = 0 - sd_fun(!!outcome_var, na.rm = T)) %>% 
    mutate(nohc = "diff")
)

smo_data <- smo_data %>% 
  left_join(y_axes)

smo_data <- smo_data %>% 
  mutate(nohc = fct_inorder(case_when(
    nohc == "1" ~ paste0("Cycling (n=", n_women, ")"),
    nohc == "0" ~ paste0("HC user (n=", n_women, ")"),
    nohc == "diff" ~ "Difference"
  )))
  
alphas <- diary %>% 
  ungroup() %>% 
  select(RCD_fab, prc_stirn_b_squished) %>% 
  mutate(RCD_fab = round(RCD_fab)) %>% 
  distinct()
ggplot(smo_data, aes(RCD_fab)) +
  # grow the Y axes to show mean±1SD
  geom_blank(aes(ymin = outcomemeanminus1sd, 
                 ymax = outcomemeanplus1sd)) +
  # show fertile window probabilities as BG color
  geom_rect(aes(xmin = RCD_fab, xmax = RCD_fab + 1,
                alpha = prc_stirn_b_squished), 
            ymin = -Inf, ymax = Inf, 
            fill = '#A45B9D', color = NA, 
            group = 1, data = alphas) +
  # vertical line at the est. day of ovulation
  geom_vline(xintercept = -14, linetype = 'dashed') +
  
  # draw samples from spline
  geom_line(aes(y = estimate__, colour = nohc, group = sample__), stat = 'identity', alpha = 0.08) +
  
  # color lines
  scale_color_manual("Contraception", 
                     values = c('black', 'red', "#4B7BB4"), 
                     guide = FALSE) +
  
  # make the alpha range fit the prc_stirn_b values
  scale_alpha(guide = F, range = c(0,0.58)) +
  
  # irrelevant when showing individual spline estimates
  # scale_fill_manual("Contraception", 
  #                   labels = c("0" = "HC user", "1" = "Cycling", "diff" = "Difference"), 
  #                   values = c("0" = 'black', '1' = 'red', "diff" = "#4B7BB4"), 
  #                   breaks = c("0", "1", "diff"), guide = FALSE) +
  
  # facet the plot by contraception 
  facet_wrap( ~ nohc, scales = "free_y", ncol = 3) +
  
  # fancy labelling for X axis
  scale_x_continuous("Days until next menstrual onset (standardized to 29-day cycle)", 
                     breaks = c(-29, -21, -14, -6,-5,-4,-3,-2, -1), 
                     labels = c("29", "22", "14\nest. ovulation", "6\n", "\n", "\n", "\n\npremens.", "\n", "1\n")) +
  
  # This variable comes from the function
  scale_y_continuous(outcome) +
  
  cowplot::theme_cowplot(font_size = 10, 
                         font_family = "Helvetica") +
  
  theme(strip.background = element_rect(fill = "#FFFFFF"),
        strip.text = element_text(size = 11))

ggsave(paste0("figures/",outcome, "_curve.png"), width = 6, height = 3)
ggsave(paste0("figures/",outcome, "_curve.pdf"), width = 6, height = 3)
sjPlot::tab_model(model2_family, ci.hyphen = ";", show.ci = 0.99, show.icc = FALSE, bpe = "mean")
  time_family
Predictors Estimates CI (99%)
Intercept 1.71 1.64;1.78
premenstrual_phase_fab:
TRUE 		-0.01 -0.08;0.05
menstruation 0.08 0.01;0.16
fertile_fab -0.01 -0.17;0.15
hormonal_contraception:
TRUE 		0.06 -0.07;0.20
premenstrual_phase_fabTRUE.hormonal_contraceptionTRUE 0.06 -0.07;0.18
menstruation.hormonal_contraceptionTRUE 0.00 -0.15;0.16
fertile_fab.hormonal_contraceptionTRUE 0.07 -0.24;0.37
Random Effects
σ2 1.27
τ00 person 0.24
τ11 person.fertile_fab 0.48
ρ01  
ρ01  
N person 793
Observations 24450
Marginal R2 / Conditional R2 0.002 / 0.158 
summary(model2_family)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: time_family ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person) 
##    Data: diary %>% filter(reasons_for_exclusion == "" | rea (Number of observations: 24450) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Group-Level Effects: 
## ~person (Number of levels: 793) 
##                            Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)                  0.49      0.02     0.45     0.54 1.00     2011     2789
## sd(fertile_fab)                0.69      0.06     0.54     0.84 1.00     1204     2107
## cor(Intercept,fertile_fab)    -0.31      0.07    -0.47    -0.13 1.00     3488     3232
## 
## Population-Level Effects: 
##                                                       Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept                                                 1.71      0.03     1.64     1.78 1.00     2155
## premenstrual_phase_fabTRUE                               -0.01      0.03    -0.08     0.05 1.00     5308
## menstruation                                              0.08      0.03     0.01     0.16 1.00     6070
## fertile_fab                                              -0.01      0.06    -0.17     0.15 1.00     4801
## hormonal_contraceptionTRUE                                0.06      0.05    -0.07     0.20 1.00     2318
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE     0.06      0.05    -0.07     0.18 1.00     5131
## menstruation:hormonal_contraceptionTRUE                   0.00      0.06    -0.15     0.16 1.00     5411
## fertile_fab:hormonal_contraceptionTRUE                    0.07      0.12    -0.24     0.37 1.00     4271
##                                                       Tail_ESS
## Intercept                                                 3115
## premenstrual_phase_fabTRUE                                3246
## menstruation                                              3577
## fertile_fab                                               3520
## hormonal_contraceptionTRUE                                3016
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE     3348
## menstruation:hormonal_contraceptionTRUE                   3506
## fertile_fab:hormonal_contraceptionTRUE                    3671
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.13      0.01     1.11     1.14 1.00     5864     2978
## 
## 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).
marginal_effects(model2_family)

custom_forest(model2_family, pars = "fertile_fab") +
    theme(axis.text.y = element_blank()) + 
    geom_hline(yintercept = 0, linetype = 'dashed', size = 1, color = "black") +
    theme(axis.text.y = element_blank(),
          axis.ticks.y = element_blank(),
          legend.position = c(0,1),
          legend.justification = c(0,1)) +
    facet_null() +
    ylab("") +
    scale_y_continuous("Est. fertile phase effect + 80% CI", breaks = -2:2)

unrelated men

model2_unrelated_men <- brm(person_is_unrelated_man ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person), data =  diary %>% filter(reasons_for_exclusion == ""), backend = "cmdstanr", family = poisson(), cores = 4, file ="brms_fits/unrelated_men")

sjPlot::tab_model(model2_unrelated_men, ci.hyphen = ";", show.ci = 0.99, show.icc = FALSE, bpe = "mean")
  person_is_unrelated_man
Predictors Incidence Rate Ratios CI (99%)
Intercept 0.20 0.13;0.27
premenstrual_phase_fab:
TRUE 		0.93 0.82;1.05
Est. menstruation 1.07 0.93;1.21
Est. fertile window prob.
(BC+i) 		0.97 0.66;1.37
hormonal_contraception:
TRUE 		0.87 0.39;1.66
premenstrual_phase_fabTRUE.hormonal_contraceptionTRUE 0.97 0.73;1.27
menstruation.hormonal_contraceptionTRUE 0.96 0.68;1.30
fertile_fab.hormonal_contraceptionTRUE 1.20 0.67;2.08
Random Effects
σ2 0.42
τ00 0.19
N person 263
Observations 12251
Marginal R2 / Conditional R2 0.001 / 0.363 
summary(model2_unrelated_men)
##  Family: poisson 
##   Links: mu = log 
## Formula: person_is_unrelated_man ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person) 
##    Data: diary %>% filter(reasons_for_exclusion == "") (Number of observations: 12251) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Group-Level Effects: 
## ~person (Number of levels: 263) 
##                            Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)                  1.72      0.11     1.47     2.02 1.01      677     1712
## sd(fertile_fab)                0.26      0.16     0.00     0.69 1.01      501     1380
## cor(Intercept,fertile_fab)     0.07      0.50    -0.97     0.99 1.00     4001     2084
## 
## Population-Level Effects: 
##                                                       Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept                                                -1.64      0.14    -2.03    -1.30 1.03      190
## premenstrual_phase_fabTRUE                               -0.08      0.05    -0.19     0.04 1.00     4423
## menstruation                                              0.06      0.05    -0.07     0.19 1.00     5082
## fertile_fab                                              -0.04      0.14    -0.42     0.32 1.00     3856
## hormonal_contraceptionTRUE                               -0.18      0.29    -0.93     0.51 1.01      291
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE    -0.04      0.11    -0.31     0.24 1.00     4182
## menstruation:hormonal_contraceptionTRUE                  -0.05      0.12    -0.39     0.26 1.00     3892
## fertile_fab:hormonal_contraceptionTRUE                    0.16      0.23    -0.41     0.73 1.00     4001
##                                                       Tail_ESS
## Intercept                                                  362
## premenstrual_phase_fabTRUE                                3061
## menstruation                                              3272
## fertile_fab                                               3128
## hormonal_contraceptionTRUE                                 570
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE     3100
## menstruation:hormonal_contraceptionTRUE                   3183
## fertile_fab:hormonal_contraceptionTRUE                    2923
## 
## 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).
marginal_effects(model2_unrelated_men)

unrelated women

model2_unrelated_women <- brm(person_is_unrelated_woman ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person), data = diary, backend = "cmdstanr", family = poisson(), cores = 4, file ="brms_fits/unrelated_woman")

sjPlot::tab_model(model2_unrelated_women, ci.hyphen = ";", show.ci = 0.99, show.icc = FALSE, bpe = "mean")
  person_is_unrelated_woman
Predictors Incidence Rate Ratios CI (99%)
Intercept 0.21 0.14;0.29
premenstrual_phase_fab:
TRUE 		0.94 0.85;1.03
Est. menstruation 0.99 0.88;1.11
Est. fertile window prob.
(BC+i) 		0.89 0.60;1.29
hormonal_contraception:
TRUE 		0.83 0.36;1.64
premenstrual_phase_fabTRUE.hormonal_contraceptionTRUE 1.13 0.88;1.47
menstruation.hormonal_contraceptionTRUE 1.36 1.04;1.79
fertile_fab.hormonal_contraceptionTRUE 1.27 0.69;2.24
Random Effects
σ2 0.71
τ00 0.20
N person 263
Observations 12251
Marginal R2 / Conditional R2 0.000 / 0.400 
summary(model2_unrelated_women)
##  Family: poisson 
##   Links: mu = log 
## Formula: person_is_unrelated_woman ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person) 
##    Data: diary (Number of observations: 12251) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Group-Level Effects: 
## ~person (Number of levels: 263) 
##                            Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)                  1.76      0.11     1.50     2.09 1.00      567     1046
## sd(fertile_fab)                0.34      0.16     0.01     0.70 1.01      520      940
## cor(Intercept,fertile_fab)     0.10      0.40    -0.89     0.97 1.00     4361     1969
## 
## Population-Level Effects: 
##                                                       Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept                                                -1.56      0.14    -1.95    -1.23 1.01      215
## premenstrual_phase_fabTRUE                               -0.07      0.04    -0.16     0.03 1.00     5170
## menstruation                                             -0.01      0.04    -0.13     0.10 1.00     5543
## fertile_fab                                              -0.12      0.15    -0.51     0.25 1.00     4331
## hormonal_contraceptionTRUE                               -0.22      0.29    -1.01     0.49 1.01      289
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE     0.12      0.10    -0.12     0.39 1.00     4002
## menstruation:hormonal_contraceptionTRUE                   0.30      0.11     0.04     0.58 1.00     4295
## fertile_fab:hormonal_contraceptionTRUE                    0.21      0.23    -0.37     0.81 1.00     4140
##                                                       Tail_ESS
## Intercept                                                  659
## premenstrual_phase_fabTRUE                                3125
## menstruation                                              3338
## fertile_fab                                               3324
## hormonal_contraceptionTRUE                                 894
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE     3415
## menstruation:hormonal_contraceptionTRUE                   3646
## fertile_fab:hormonal_contraceptionTRUE                    3442
## 
## 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).
marginal_effects(model2_unrelated_women)