Robustness analyses
Julie Driebe & Ruben Arslan
28 February 2019
Robustness analyses
knitr::opts_chunk$set(error = TRUE)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 = 10, fig.height = 10)
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)
options(brms.backend = "cmdstanr",
mc.cores = 4)time spent with family
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, file = "brms_fits/model2_family3")
do_model(model2_family, diary %>% mutate(reasons_for_exclusion = str_replace_all(reasons_for_exclusion, "not_single, ", "")), model_prefix = "family_")Robustness checks
Basic model
library(effects); library(lme4); library(lmerTest); library(knitr); library(brms)
summary(model)## 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).outcome <- model$formula[[1]][[2]]Full random slopes
add_mod_formula = update.formula(
formula(model),
new = as.formula(. ~ . - (1 + fertile_fab | person) + (1 + premenstrual_phase_fab + menstruation + fertile_fab | person ))) # reorder so that the triptych looks nice
full_random_slopes = update(model, recompile = TRUE, iter = 6000, warmup = 4000, cores = 4, formula = add_mod_formula, newdata = diary %>% filter(reasons_for_exclusion == ""), file= paste0("brms_fits/", model_prefix, "rob_check_full_random_slopes"))
summary(full_random_slopes)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: time_family ~ premenstrual_phase_fab + menstruation + fertile_fab + hormonal_contraception + (1 + premenstrual_phase_fab + menstruation + fertile_fab | person) + premenstrual_phase_fab:hormonal_contraception + menstruation:hormonal_contraception + fertile_fab:hormonal_contraception
## Data: diary %>% filter(reasons_for_exclusion == "") (Number of observations: 24450)
## Samples: 4 chains, each with iter = 6000; warmup = 4000; thin = 1;
## total post-warmup samples = 8000
##
## 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.51 0.02 0.45 0.57 1.00 2910 4345
## sd(premenstrual_phase_fabTRUE) 0.32 0.04 0.22 0.41 1.00 1162 2178
## sd(menstruation) 0.40 0.04 0.28 0.51 1.00 1313 2325
## sd(fertile_fab) 0.76 0.08 0.56 0.95 1.00 975 2601
## cor(Intercept,premenstrual_phase_fabTRUE) -0.23 0.09 -0.44 0.04 1.00 2123 4471
## cor(Intercept,menstruation) -0.20 0.09 -0.42 0.07 1.00 1756 3644
## cor(premenstrual_phase_fabTRUE,menstruation) 0.41 0.13 0.03 0.70 1.01 974 1963
## cor(Intercept,fertile_fab) -0.36 0.08 -0.54 -0.14 1.00 2347 4212
## cor(premenstrual_phase_fabTRUE,fertile_fab) 0.37 0.12 0.01 0.65 1.00 739 1885
## cor(menstruation,fertile_fab) 0.27 0.13 -0.13 0.56 1.00 748 1410
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept 1.71 0.03 1.64 1.79 1.00 3151
## premenstrual_phase_fabTRUE -0.02 0.03 -0.09 0.06 1.00 5265
## menstruation 0.08 0.04 -0.01 0.17 1.00 5878
## fertile_fab -0.01 0.06 -0.18 0.15 1.00 5406
## hormonal_contraceptionTRUE 0.07 0.05 -0.07 0.21 1.00 3080
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.06 0.05 -0.08 0.20 1.00 4753
## menstruation:hormonal_contraceptionTRUE 0.01 0.07 -0.16 0.18 1.00 5772
## fertile_fab:hormonal_contraceptionTRUE 0.07 0.12 -0.25 0.38 1.00 4832
## Tail_ESS
## Intercept 4803
## premenstrual_phase_fabTRUE 6135
## menstruation 6119
## fertile_fab 5329
## hormonal_contraceptionTRUE 4257
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 6272
## menstruation:hormonal_contraceptionTRUE 5629
## fertile_fab:hormonal_contraceptionTRUE 5444
##
## Family Specific Parameters:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.12 0.01 1.10 1.13 1.00 8399 4657
##
## 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).Full random slopes, no correlations
add_mod_formula = update.formula(
formula(model),
new = as.formula(. ~ . - (1 + fertile_fab | person) + (1 + premenstrual_phase_fab + menstruation + fertile_fab || person ))) # reorder so that the triptych looks nice
full_random_slopes_no_cor = update(model, recompile = TRUE, iter = 6000, warmup = 4000, cores = 4, formula = add_mod_formula, newdata = diary %>% filter(reasons_for_exclusion == ""), file= paste0("brms_fits/", model_prefix, "rob_check_full_random_slopes_no_cor"))
summary(full_random_slopes_no_cor)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: time_family ~ premenstrual_phase_fab + menstruation + fertile_fab + hormonal_contraception + (1 + premenstrual_phase_fab + menstruation + fertile_fab || person) + premenstrual_phase_fab:hormonal_contraception + menstruation:hormonal_contraception + fertile_fab:hormonal_contraception
## Data: diary %>% filter(reasons_for_exclusion == "") (Number of observations: 24450)
## Samples: 4 chains, each with iter = 6000; warmup = 4000; thin = 1;
## total post-warmup samples = 8000
##
## 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.46 0.02 0.42 0.50 1.00 3389 4898
## sd(premenstrual_phase_fabTRUE) 0.26 0.03 0.18 0.33 1.00 2140 3490
## sd(menstruation) 0.35 0.04 0.25 0.44 1.00 2645 4655
## sd(fertile_fab) 0.58 0.06 0.41 0.72 1.00 1933 3265
##
## 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 3861
## premenstrual_phase_fabTRUE -0.02 0.03 -0.09 0.05 1.00 8864
## menstruation 0.08 0.03 -0.01 0.17 1.00 8537
## fertile_fab -0.02 0.06 -0.17 0.14 1.00 8306
## hormonal_contraceptionTRUE 0.07 0.05 -0.06 0.20 1.00 3736
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.06 0.05 -0.08 0.19 1.00 7263
## menstruation:hormonal_contraceptionTRUE 0.00 0.06 -0.16 0.17 1.00 7634
## fertile_fab:hormonal_contraceptionTRUE 0.06 0.11 -0.24 0.35 1.00 7404
## Tail_ESS
## Intercept 5268
## premenstrual_phase_fabTRUE 6931
## menstruation 6777
## fertile_fab 6725
## hormonal_contraceptionTRUE 4868
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 6016
## menstruation:hormonal_contraceptionTRUE 6302
## fertile_fab:hormonal_contraceptionTRUE 6655
##
## Family Specific Parameters:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.12 0.01 1.11 1.13 1.00 10762 6072
##
## 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).Without random slopes
no_varying_slope = update(model, recompile = TRUE, formula = . ~ . - (1 + fertile_fab | person) + (1|person), cores = 4, file = paste0("brms_fits/", model_prefix, "rob_check_no_varying_slope"))
summary(no_varying_slope)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: time_family ~ premenstrual_phase_fab + menstruation + fertile_fab + hormonal_contraception + (1 | person) + premenstrual_phase_fab:hormonal_contraception + menstruation:hormonal_contraception + fertile_fab:hormonal_contraception
## 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.47 0.02 0.43 0.51 1.00 657 1096
##
## 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 870
## premenstrual_phase_fabTRUE -0.01 0.03 -0.08 0.05 1.00 2413
## menstruation 0.08 0.03 0.01 0.16 1.00 2630
## fertile_fab -0.01 0.06 -0.16 0.13 1.00 2329
## hormonal_contraceptionTRUE 0.07 0.05 -0.06 0.20 1.01 833
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.06 0.05 -0.07 0.18 1.00 1867
## menstruation:hormonal_contraceptionTRUE -0.00 0.06 -0.15 0.15 1.00 2130
## fertile_fab:hormonal_contraceptionTRUE 0.05 0.11 -0.23 0.32 1.00 1851
## Tail_ESS
## Intercept 1493
## premenstrual_phase_fabTRUE 2664
## menstruation 2884
## fertile_fab 2671
## hormonal_contraceptionTRUE 1487
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 2338
## menstruation:hormonal_contraceptionTRUE 2685
## fertile_fab:hormonal_contraceptionTRUE 2241
##
## Family Specific Parameters:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.14 0.01 1.12 1.15 1.00 7181 2760
##
## 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).Only women with outcome variability
with_outcome_var = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "") %>% group_by(session) %>% filter(sd(!!outcome, na.rm = T) > 0), file = paste0("brms_fits/", model_prefix, "rob_check_with_outcome_var"))
summary(with_outcome_var)## 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 == "") %>% (Number of observations: 24326)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 775)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.49 0.02 0.44 0.53 1.00 2216 2943
## sd(fertile_fab) 0.70 0.06 0.55 0.84 1.00 1307 2166
## cor(Intercept,fertile_fab) -0.32 0.07 -0.47 -0.13 1.00 4274 3331
##
## 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 1992
## premenstrual_phase_fabTRUE -0.01 0.03 -0.08 0.05 1.00 5374
## menstruation 0.08 0.03 0.01 0.16 1.00 6014
## fertile_fab -0.01 0.06 -0.17 0.16 1.00 4759
## hormonal_contraceptionTRUE 0.07 0.05 -0.07 0.21 1.00 1981
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.06 0.05 -0.06 0.18 1.00 4337
## menstruation:hormonal_contraceptionTRUE -0.00 0.06 -0.15 0.15 1.00 5884
## fertile_fab:hormonal_contraceptionTRUE 0.07 0.12 -0.24 0.38 1.00 4080
## Tail_ESS
## Intercept 2599
## premenstrual_phase_fabTRUE 3305
## menstruation 2977
## fertile_fab 3409
## hormonal_contraceptionTRUE 2613
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3109
## menstruation:hormonal_contraceptionTRUE 3263
## fertile_fab:hormonal_contraceptionTRUE 2981
##
## Family Specific Parameters:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.13 0.01 1.12 1.14 1.00 5312 2547
##
## 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).Ordinal model
if(identical(model$family$family, "gaussian")) {
ordinal_model = update(model, recompile = TRUE, family = cumulative(), newdata = diary %>% mutate(time_family = as.integer(1 + time_family)) %>% filter(reasons_for_exclusion == ""), file= paste0("brms_fits/", model_prefix, "rob_check_full_ordinal_model"))
summary(ordinal_model)
} ## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: time_family ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person)
## Data: diary %>% mutate(time_family = as.integer(1 + time (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.87 0.03 0.80 0.95 1.00 1631 2916
## sd(fertile_fab) 1.16 0.10 0.93 1.40 1.00 1285 2138
## cor(Intercept,fertile_fab) -0.29 0.07 -0.46 -0.11 1.00 3753 2995
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept[1] -1.50 0.05 -1.63 -1.36 1.00 2010
## Intercept[2] -0.23 0.05 -0.36 -0.10 1.00 1964
## Intercept[3] 1.22 0.05 1.09 1.34 1.00 2026
## Intercept[4] 2.60 0.05 2.47 2.74 1.00 2169
## premenstrual_phase_fabTRUE -0.03 0.04 -0.13 0.07 1.00 5368
## menstruation 0.12 0.05 0.00 0.24 1.00 6239
## fertile_fab -0.02 0.10 -0.29 0.25 1.00 3936
## hormonal_contraceptionTRUE 0.09 0.09 -0.14 0.33 1.00 2133
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.11 0.08 -0.09 0.30 1.00 5031
## menstruation:hormonal_contraceptionTRUE 0.02 0.09 -0.21 0.27 1.00 5385
## fertile_fab:hormonal_contraceptionTRUE 0.15 0.20 -0.35 0.67 1.00 3949
## Tail_ESS
## Intercept[1] 2626
## Intercept[2] 2461
## Intercept[3] 2593
## Intercept[4] 2709
## premenstrual_phase_fabTRUE 3124
## menstruation 3726
## fertile_fab 3175
## hormonal_contraceptionTRUE 2453
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3386
## menstruation:hormonal_contraceptionTRUE 3134
## fertile_fab:hormonal_contraceptionTRUE 3390
##
## Family Specific Parameters:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 1.00 4000 4000
##
## 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).Gaussian model
if(identical(model$family$family, "poisson")) {
gaussian_model = update(model, recompile = TRUE, family = gaussian(), file= paste0("brms_fits/", model_prefix, "rob_check_full_ordinal_model"))
summary(gaussian_model)
}No exclusions
no_exclusions = update(model, recompile = TRUE, cores = 4, newdata = diary, file=paste0("brms_fits/", model_prefix, "rob_check_no_exclusions"))
summary(no_exclusions)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: time_family ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person)
## Data: diary (Number of observations: 30472)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 964)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.48 0.01 0.44 0.51 1.00 1545 2553
## sd(fertile_fab) 0.68 0.05 0.54 0.81 1.00 1064 1526
## cor(Intercept,fertile_fab) -0.30 0.06 -0.45 -0.14 1.00 2830 3072
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept 1.71 0.03 1.65 1.78 1.00 1116
## premenstrual_phase_fabTRUE -0.01 0.02 -0.07 0.05 1.00 3371
## menstruation 0.10 0.03 0.03 0.17 1.00 3567
## fertile_fab 0.01 0.06 -0.14 0.15 1.00 2645
## hormonal_contraceptionTRUE 0.07 0.05 -0.04 0.18 1.00 1230
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.07 0.04 -0.04 0.17 1.00 2818
## menstruation:hormonal_contraceptionTRUE -0.05 0.05 -0.17 0.08 1.00 3716
## fertile_fab:hormonal_contraceptionTRUE 0.02 0.10 -0.23 0.29 1.00 2805
## Tail_ESS
## Intercept 2212
## premenstrual_phase_fabTRUE 2958
## menstruation 2910
## fertile_fab 2786
## hormonal_contraceptionTRUE 2107
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3429
## menstruation:hormonal_contraceptionTRUE 3289
## fertile_fab:hormonal_contraceptionTRUE 2898
##
## Family Specific Parameters:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.12 0.00 1.11 1.13 1.00 4924 2857
##
## 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).Exclude cycles shorter than 25 days
rob_check_mens_length = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "") %>% filter(cycle_length >= 25), file=paste0("brms_fits/", model_prefix, "rob_check_mens_length"))
summary(rob_check_mens_length)## 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 == "") %>% (Number of observations: 22506)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 766)
## 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 1888 2432
## sd(fertile_fab) 0.76 0.06 0.60 0.92 1.00 1090 2232
## cor(Intercept,fertile_fab) -0.31 0.06 -0.46 -0.13 1.00 3232 3064
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept 1.71 0.03 1.63 1.79 1.00 1757
## premenstrual_phase_fabTRUE -0.00 0.03 -0.07 0.06 1.00 3620
## menstruation 0.08 0.03 -0.01 0.16 1.00 3505
## fertile_fab -0.03 0.07 -0.22 0.15 1.00 2768
## hormonal_contraceptionTRUE 0.05 0.05 -0.09 0.19 1.00 1658
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.07 0.05 -0.05 0.20 1.00 3125
## menstruation:hormonal_contraceptionTRUE 0.03 0.06 -0.12 0.19 1.00 2725
## fertile_fab:hormonal_contraceptionTRUE 0.13 0.12 -0.18 0.44 1.00 2123
## Tail_ESS
## Intercept 2357
## premenstrual_phase_fabTRUE 3107
## menstruation 3057
## fertile_fab 2655
## hormonal_contraceptionTRUE 2644
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3195
## menstruation:hormonal_contraceptionTRUE 2219
## fertile_fab:hormonal_contraceptionTRUE 2889
##
## Family Specific Parameters:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.12 0.01 1.11 1.14 1.00 5271 2831
##
## 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).Exclude women who report cycles outside 25-35
rob_check_selfrepcycle_length = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "") %>% filter(menstruation_length >= 25, menstruation_length <= 35), file=paste0("brms_fits/", model_prefix, "rob_check_selfrepcycle_length"))
summary(rob_check_selfrepcycle_length)## 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 == "") %>% (Number of observations: 23048)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 745)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.49 0.02 0.44 0.53 1.00 2265 2386
## sd(fertile_fab) 0.68 0.06 0.51 0.84 1.00 1211 1922
## cor(Intercept,fertile_fab) -0.27 0.07 -0.45 -0.07 1.00 3779 3078
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept 1.71 0.03 1.63 1.79 1.00 1902
## premenstrual_phase_fabTRUE -0.00 0.03 -0.07 0.06 1.00 4988
## menstruation 0.09 0.03 0.01 0.16 1.00 4908
## fertile_fab 0.02 0.06 -0.14 0.18 1.00 3882
## hormonal_contraceptionTRUE 0.06 0.05 -0.08 0.20 1.00 2078
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.05 0.05 -0.07 0.17 1.00 4804
## menstruation:hormonal_contraceptionTRUE -0.02 0.06 -0.16 0.14 1.00 5097
## fertile_fab:hormonal_contraceptionTRUE 0.07 0.12 -0.25 0.39 1.00 4338
## Tail_ESS
## Intercept 2468
## premenstrual_phase_fabTRUE 3479
## menstruation 3269
## fertile_fab 3327
## hormonal_contraceptionTRUE 2632
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3686
## menstruation:hormonal_contraceptionTRUE 3485
## fertile_fab:hormonal_contraceptionTRUE 3224
##
## Family Specific Parameters:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.12 0.01 1.11 1.14 1.00 5560 2442
##
## 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).Exclude women who report cycles varying more than 1-2 days
rob_check_menstruation_regularity = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "") %>% filter(menstruation_regularity <= 2), file=paste0("brms_fits/", model_prefix, "rob_check_menstruation_regularity"))
summary(rob_check_menstruation_regularity)## 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 == "") %>% (Number of observations: 17089)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 557)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.48 0.02 0.43 0.53 1.00 1758 2696
## sd(fertile_fab) 0.64 0.07 0.43 0.82 1.00 1236 1814
## cor(Intercept,fertile_fab) -0.34 0.08 -0.53 -0.11 1.00 4745 3587
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept 1.69 0.04 1.60 1.78 1.00 2176
## premenstrual_phase_fabTRUE 0.01 0.03 -0.07 0.09 1.00 5131
## menstruation 0.11 0.04 0.02 0.21 1.00 6109
## fertile_fab -0.00 0.08 -0.22 0.19 1.00 4585
## hormonal_contraceptionTRUE 0.11 0.06 -0.04 0.26 1.00 2578
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.03 0.05 -0.11 0.16 1.00 4874
## menstruation:hormonal_contraceptionTRUE 0.00 0.06 -0.16 0.17 1.00 5439
## fertile_fab:hormonal_contraceptionTRUE 0.08 0.13 -0.26 0.41 1.00 4158
## Tail_ESS
## Intercept 2766
## premenstrual_phase_fabTRUE 3542
## menstruation 3485
## fertile_fab 3559
## hormonal_contraceptionTRUE 3041
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3254
## menstruation:hormonal_contraceptionTRUE 3590
## fertile_fab:hormonal_contraceptionTRUE 3359
##
## Family Specific Parameters:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.13 0.01 1.12 1.15 1.00 5616 2878
##
## 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).Exclude women who report uncertainty larger than ±2 days
Across three items (last menstrual onset, cycle length, menstruation regularity).
diary$menstruation_certainty <- rowMeans(diary %>% select(menstruation_last_certainty, menstruation_regularity_certainty, menstruation_length_certainty), na.rm = T)
rob_check_menstruation_certainty = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "") %>% filter(menstruation_certainty <= 3), file=paste0("brms_fits/", model_prefix, "rob_check_menstruation_certainty"))
summary(rob_check_menstruation_certainty)## 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 == "") %>% (Number of observations: 22385)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 727)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.50 0.02 0.45 0.55 1.00 1403 2128
## sd(fertile_fab) 0.68 0.06 0.52 0.85 1.00 934 1836
## cor(Intercept,fertile_fab) -0.31 0.07 -0.48 -0.13 1.00 3421 3223
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept 1.72 0.03 1.64 1.79 1.00 1564
## premenstrual_phase_fabTRUE -0.02 0.03 -0.09 0.05 1.00 3940
## menstruation 0.08 0.03 -0.00 0.16 1.00 4261
## fertile_fab -0.02 0.07 -0.19 0.15 1.00 3406
## hormonal_contraceptionTRUE 0.07 0.06 -0.08 0.21 1.00 1389
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.06 0.05 -0.07 0.19 1.00 3083
## menstruation:hormonal_contraceptionTRUE 0.01 0.06 -0.14 0.17 1.00 3636
## fertile_fab:hormonal_contraceptionTRUE 0.06 0.12 -0.28 0.37 1.00 2531
## Tail_ESS
## Intercept 2301
## premenstrual_phase_fabTRUE 3263
## menstruation 2987
## fertile_fab 3414
## hormonal_contraceptionTRUE 2226
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 2960
## menstruation:hormonal_contraceptionTRUE 3046
## fertile_fab:hormonal_contraceptionTRUE 2615
##
## Family Specific Parameters:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.13 0.01 1.12 1.15 1.00 6495 2619
##
## 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).Exclude women who report not feeling healthy
rob_check_health = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "") %>% filter(self_rated_health <= 2), file=paste0("brms_fits/", model_prefix, "rob_check_health"))
summary(rob_check_health)## 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 == "") %>% (Number of observations: 21452)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 687)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.49 0.02 0.44 0.54 1.00 2102 2880
## sd(fertile_fab) 0.68 0.06 0.51 0.84 1.00 1164 1963
## cor(Intercept,fertile_fab) -0.32 0.07 -0.50 -0.12 1.00 4163 3227
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept 1.73 0.03 1.65 1.81 1.00 2313
## premenstrual_phase_fabTRUE -0.02 0.03 -0.09 0.05 1.00 5064
## menstruation 0.09 0.03 0.00 0.17 1.00 5539
## fertile_fab -0.05 0.07 -0.24 0.12 1.00 4221
## hormonal_contraceptionTRUE 0.05 0.06 -0.10 0.19 1.00 2640
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.08 0.05 -0.05 0.21 1.00 5220
## menstruation:hormonal_contraceptionTRUE 0.00 0.06 -0.16 0.16 1.00 5133
## fertile_fab:hormonal_contraceptionTRUE 0.12 0.13 -0.20 0.46 1.00 4197
## Tail_ESS
## Intercept 2882
## premenstrual_phase_fabTRUE 3358
## menstruation 3276
## fertile_fab 3234
## hormonal_contraceptionTRUE 2940
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3400
## menstruation:hormonal_contraceptionTRUE 3272
## fertile_fab:hormonal_contraceptionTRUE 3284
##
## Family Specific Parameters:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.13 0.01 1.12 1.15 1.00 5857 3043
##
## 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).Exclude women trying for pregnancy
pregnant_trying = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "") %>% filter(pregnant_trying < 4), file=paste0("brms_fits/", model_prefix, "rob_check_pregnant_trying"))
summary(pregnant_trying)## 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 == "") %>% (Number of observations: 23782)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 767)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.49 0.02 0.45 0.53 1.00 1445 2508
## sd(fertile_fab) 0.69 0.06 0.54 0.85 1.00 1080 1857
## cor(Intercept,fertile_fab) -0.30 0.07 -0.46 -0.10 1.00 2968 2775
##
## 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 1272
## premenstrual_phase_fabTRUE -0.01 0.03 -0.08 0.06 1.00 3501
## menstruation 0.08 0.03 0.00 0.16 1.00 3588
## fertile_fab 0.00 0.06 -0.16 0.17 1.00 2976
## hormonal_contraceptionTRUE 0.06 0.05 -0.07 0.19 1.00 1475
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.06 0.05 -0.07 0.17 1.00 3456
## menstruation:hormonal_contraceptionTRUE -0.01 0.06 -0.15 0.14 1.00 4004
## fertile_fab:hormonal_contraceptionTRUE 0.06 0.12 -0.26 0.38 1.00 3098
## Tail_ESS
## Intercept 2203
## premenstrual_phase_fabTRUE 2905
## menstruation 2683
## fertile_fab 3421
## hormonal_contraceptionTRUE 1907
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3233
## menstruation:hormonal_contraceptionTRUE 3309
## fertile_fab:hormonal_contraceptionTRUE 3433
##
## 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 5066 2570
##
## 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).Include women living with parents
rob_check_parents = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "" | reasons_for_exclusion == "living_with_parents, "), file=paste0("brms_fits/", model_prefix, "rob_check_parents"))
summary(rob_check_parents)## 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: 26793)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 870)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.49 0.02 0.45 0.52 1.00 2079 2705
## sd(fertile_fab) 0.68 0.06 0.54 0.83 1.00 1154 2317
## cor(Intercept,fertile_fab) -0.30 0.06 -0.45 -0.13 1.00 3886 3305
##
## 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 1957
## premenstrual_phase_fabTRUE -0.01 0.02 -0.07 0.05 1.00 5248
## menstruation 0.09 0.03 0.02 0.16 1.00 5426
## fertile_fab -0.01 0.06 -0.16 0.14 1.00 4050
## hormonal_contraceptionTRUE 0.07 0.05 -0.06 0.19 1.00 2002
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.05 0.04 -0.06 0.17 1.00 4937
## menstruation:hormonal_contraceptionTRUE -0.03 0.05 -0.17 0.11 1.00 5116
## fertile_fab:hormonal_contraceptionTRUE 0.05 0.11 -0.23 0.36 1.00 4392
## Tail_ESS
## Intercept 2488
## premenstrual_phase_fabTRUE 2892
## menstruation 3268
## fertile_fab 3247
## hormonal_contraceptionTRUE 2547
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3363
## menstruation:hormonal_contraceptionTRUE 3787
## fertile_fab:hormonal_contraceptionTRUE 3099
##
## Family Specific Parameters:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.12 0.01 1.11 1.14 1.00 5170 2615
##
## 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).Exclude non-singles
if(outcome != "person_is_related_man_seen") {
rob_check_singles = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "", hetero_relationship == 0, relationship_status == 1), file=paste0("brms_fits/", model_prefix, "rob_check_singles"))
summary(rob_check_singles)
}## 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 == "", hete (Number of observations: 8512)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 262)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.48 0.03 0.41 0.56 1.00 872 2077
## sd(fertile_fab) 0.88 0.10 0.61 1.14 1.00 740 1549
## cor(Intercept,fertile_fab) -0.31 0.10 -0.55 -0.02 1.00 1189 2554
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept 1.63 0.04 1.52 1.75 1.01 861
## premenstrual_phase_fabTRUE 0.05 0.04 -0.06 0.15 1.00 2407
## menstruation 0.13 0.05 0.01 0.25 1.00 2832
## fertile_fab 0.05 0.11 -0.24 0.32 1.00 1446
## hormonal_contraceptionTRUE 0.08 0.11 -0.18 0.34 1.00 941
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.03 0.10 -0.24 0.28 1.00 1951
## menstruation:hormonal_contraceptionTRUE 0.00 0.11 -0.28 0.30 1.00 2211
## fertile_fab:hormonal_contraceptionTRUE -0.13 0.26 -0.79 0.54 1.00 1322
## Tail_ESS
## Intercept 1941
## premenstrual_phase_fabTRUE 2879
## menstruation 2449
## fertile_fab 2397
## hormonal_contraceptionTRUE 1170
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 2610
## menstruation:hormonal_contraceptionTRUE 2614
## fertile_fab:hormonal_contraceptionTRUE 1976
##
## Family Specific Parameters:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.14 0.01 1.12 1.16 1.00 4223 2649
##
## 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).Moderation by weekend
weekend = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "") %>% filter(weekend == 1), file=paste0("brms_fits/", model_prefix, "rob_check_weekend"))
weekday = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "") %>% filter(weekend == 0), file=paste0("brms_fits/", model_prefix, "rob_check_weekday"))
summary(weekend)## 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 == "") %>% (Number of observations: 10182)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 780)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.55 0.02 0.50 0.62 1.00 1697 2313
## sd(fertile_fab) 0.99 0.09 0.74 1.23 1.00 1153 2002
## cor(Intercept,fertile_fab) -0.40 0.07 -0.57 -0.19 1.00 2834 3086
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept 1.88 0.04 1.78 1.98 1.00 2259
## premenstrual_phase_fabTRUE -0.08 0.04 -0.19 0.03 1.00 4046
## menstruation 0.02 0.05 -0.10 0.15 1.00 4403
## fertile_fab -0.10 0.10 -0.37 0.17 1.00 2934
## hormonal_contraceptionTRUE 0.09 0.07 -0.10 0.29 1.00 1971
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.11 0.08 -0.08 0.30 1.00 3176
## menstruation:hormonal_contraceptionTRUE 0.07 0.09 -0.17 0.32 1.00 3675
## fertile_fab:hormonal_contraceptionTRUE -0.03 0.19 -0.51 0.49 1.00 2645
## Tail_ESS
## Intercept 2687
## premenstrual_phase_fabTRUE 3306
## menstruation 3146
## fertile_fab 2923
## hormonal_contraceptionTRUE 2302
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3135
## menstruation:hormonal_contraceptionTRUE 3052
## fertile_fab:hormonal_contraceptionTRUE 2967
##
## Family Specific Parameters:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.16 0.01 1.14 1.18 1.00 4851 3129
##
## 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(weekday)## 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 == "") %>% (Number of observations: 14268)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 792)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.50 0.02 0.45 0.55 1.00 1875 2574
## sd(fertile_fab) 0.64 0.08 0.40 0.84 1.00 824 1357
## cor(Intercept,fertile_fab) -0.21 0.09 -0.43 0.07 1.00 3142 2269
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept 1.59 0.03 1.51 1.68 1.00 2851
## premenstrual_phase_fabTRUE 0.03 0.03 -0.05 0.11 1.00 4989
## menstruation 0.11 0.04 0.01 0.21 1.00 5933
## fertile_fab 0.03 0.07 -0.17 0.22 1.00 5123
## hormonal_contraceptionTRUE 0.05 0.06 -0.10 0.20 1.00 2805
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.01 0.06 -0.13 0.17 1.00 4732
## menstruation:hormonal_contraceptionTRUE -0.04 0.07 -0.22 0.14 1.00 5095
## fertile_fab:hormonal_contraceptionTRUE 0.14 0.14 -0.21 0.52 1.00 4400
## Tail_ESS
## Intercept 2520
## premenstrual_phase_fabTRUE 3134
## menstruation 3028
## fertile_fab 3361
## hormonal_contraceptionTRUE 2676
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3487
## menstruation:hormonal_contraceptionTRUE 3578
## fertile_fab:hormonal_contraceptionTRUE 3046
##
## Family Specific Parameters:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.07 0.01 1.05 1.09 1.00 4103 2773
##
## 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).Moderation by age group
older25 = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "") %>% filter(age > 25), file=paste0("brms_fits/", model_prefix, "rob_check_older25"))
younger25 = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "") %>% filter(age <= 25), file=paste0("brms_fits/", model_prefix, "rob_check_younger25"))
summary(younger25)## 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 == "") %>% (Number of observations: 13727)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 446)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.46 0.02 0.41 0.52 1.00 1077 1563
## sd(fertile_fab) 0.62 0.09 0.37 0.85 1.00 945 1862
## cor(Intercept,fertile_fab) -0.25 0.11 -0.50 0.06 1.00 2962 2914
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept 1.72 0.04 1.61 1.83 1.00 1380
## premenstrual_phase_fabTRUE 0.01 0.04 -0.09 0.10 1.00 3323
## menstruation 0.06 0.04 -0.05 0.17 1.00 3321
## fertile_fab 0.05 0.09 -0.19 0.28 1.00 2848
## hormonal_contraceptionTRUE 0.09 0.07 -0.07 0.26 1.00 1432
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.06 0.06 -0.10 0.22 1.00 3054
## menstruation:hormonal_contraceptionTRUE 0.06 0.07 -0.13 0.24 1.00 3585
## fertile_fab:hormonal_contraceptionTRUE -0.01 0.15 -0.39 0.39 1.00 2651
## Tail_ESS
## Intercept 2066
## premenstrual_phase_fabTRUE 2931
## menstruation 2696
## fertile_fab 3104
## hormonal_contraceptionTRUE 2017
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3139
## menstruation:hormonal_contraceptionTRUE 3407
## fertile_fab:hormonal_contraceptionTRUE 3067
##
## Family Specific Parameters:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.17 0.01 1.15 1.19 1.00 4991 2874
##
## 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(older25)## 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 == "") %>% (Number of observations: 10723)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 347)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.52 0.03 0.46 0.60 1.00 1055 1763
## sd(fertile_fab) 0.76 0.08 0.55 0.98 1.00 877 1332
## cor(Intercept,fertile_fab) -0.35 0.09 -0.56 -0.11 1.00 2252 2903
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept 1.69 0.04 1.59 1.81 1.00 805
## premenstrual_phase_fabTRUE -0.03 0.03 -0.11 0.06 1.00 2781
## menstruation 0.11 0.04 0.01 0.21 1.00 3103
## fertile_fab -0.07 0.09 -0.29 0.16 1.00 1780
## hormonal_contraceptionTRUE -0.03 0.09 -0.29 0.20 1.00 711
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.04 0.08 -0.17 0.24 1.00 2519
## menstruation:hormonal_contraceptionTRUE -0.12 0.10 -0.37 0.13 1.00 2977
## fertile_fab:hormonal_contraceptionTRUE 0.21 0.21 -0.30 0.76 1.00 1525
## Tail_ESS
## Intercept 1416
## premenstrual_phase_fabTRUE 2917
## menstruation 3028
## fertile_fab 2508
## hormonal_contraceptionTRUE 916
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 2582
## menstruation:hormonal_contraceptionTRUE 2803
## fertile_fab:hormonal_contraceptionTRUE 2460
##
## Family Specific Parameters:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.07 0.01 1.06 1.09 1.00 5306 2692
##
## 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).Excluding women who used psychopharmacological, hormonal, or antibiotic medication
diary <- diary %>% filter(reasons_for_exclusion == "") %>% ungroup() %>%
mutate(medication_psychopharmacological = if_else(is.na(medication_psychopharmacological), 0, medication_psychopharmacological),
medication_hormonal = if_else(is.na(medication_hormonal), 0, medication_hormonal),
medication_antibiotics = if_else(is.na(medication_antibiotics), 0, medication_antibiotics))
exclude_meds = update(model, recompile = TRUE, cores = 4, newdata= diary %>% filter(reasons_for_exclusion == "") %>% filter(medication_psychopharmacological == 0, medication_hormonal == 0, medication_antibiotics == 0), file=paste0("brms_fits/", model_prefix, "rob_check_meds"))
summary(exclude_meds)## 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 == "") %>% (Number of observations: 21768)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 706)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.50 0.02 0.45 0.55 1.00 1564 2394
## sd(fertile_fab) 0.71 0.06 0.55 0.87 1.00 1306 2312
## cor(Intercept,fertile_fab) -0.31 0.07 -0.46 -0.11 1.00 3892 3273
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept 1.72 0.03 1.63 1.80 1.00 2331
## premenstrual_phase_fabTRUE -0.01 0.03 -0.08 0.06 1.00 5341
## menstruation 0.08 0.03 0.01 0.16 1.00 5495
## fertile_fab 0.01 0.07 -0.17 0.19 1.00 4538
## hormonal_contraceptionTRUE 0.06 0.06 -0.08 0.20 1.00 2421
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.04 0.05 -0.09 0.17 1.00 4820
## menstruation:hormonal_contraceptionTRUE -0.01 0.06 -0.17 0.13 1.00 5256
## fertile_fab:hormonal_contraceptionTRUE 0.04 0.13 -0.28 0.37 1.00 4159
## Tail_ESS
## Intercept 2722
## premenstrual_phase_fabTRUE 3590
## menstruation 3358
## fertile_fab 3581
## hormonal_contraceptionTRUE 2622
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3552
## menstruation:hormonal_contraceptionTRUE 3186
## fertile_fab:hormonal_contraceptionTRUE 3576
##
## Family Specific Parameters:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.12 0.01 1.10 1.13 1.00 5125 2741
##
## 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).Excluding hypothesis guessers
hyp_guess = update(model, recompile = TRUE, cores = 4,newdata = diary %>% filter(reasons_for_exclusion == "") %>% filter(hypothesis_guess_topic == "no_guess"), file=paste0("brms_fits/", model_prefix, "_rob_check_hyp_guess"))
summary(hyp_guess)## 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 == "") %>% (Number of observations: 20529)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 680)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.51 0.02 0.46 0.56 1.00 1291 2271
## sd(fertile_fab) 0.71 0.06 0.55 0.87 1.00 1245 2012
## cor(Intercept,fertile_fab) -0.32 0.07 -0.49 -0.15 1.00 3210 3214
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept 1.70 0.03 1.61 1.78 1.00 1520
## premenstrual_phase_fabTRUE -0.00 0.03 -0.08 0.07 1.00 3986
## menstruation 0.10 0.03 0.02 0.18 1.00 4214
## fertile_fab 0.07 0.07 -0.11 0.25 1.00 2996
## hormonal_contraceptionTRUE 0.06 0.06 -0.08 0.21 1.00 1449
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.03 0.05 -0.10 0.16 1.00 3421
## menstruation:hormonal_contraceptionTRUE -0.01 0.06 -0.16 0.15 1.00 4254
## fertile_fab:hormonal_contraceptionTRUE -0.02 0.13 -0.34 0.32 1.00 2964
## Tail_ESS
## Intercept 2483
## premenstrual_phase_fabTRUE 2688
## menstruation 3458
## fertile_fab 3086
## hormonal_contraceptionTRUE 1975
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 2985
## menstruation:hormonal_contraceptionTRUE 2953
## fertile_fab:hormonal_contraceptionTRUE 2631
##
## Family Specific Parameters:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.12 0.01 1.11 1.13 1.00 5955 2152
##
## 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).Excluding women using awareness-based methods
unaware = update(model, recompile = TRUE, cores = 4,newdata = diary %>% filter(reasons_for_exclusion == "") %>% filter(contraception_approach != "awareness", contraception_app == 0), file=paste0("brms_fits/", model_prefix, "rob_check_unaware"))
summary(unaware)## 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 == "") %>% (Number of observations: 14817)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 499)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.49 0.02 0.44 0.54 1.00 2048 2924
## sd(fertile_fab) 0.79 0.07 0.60 0.98 1.00 1384 2470
## cor(Intercept,fertile_fab) -0.31 0.08 -0.50 -0.09 1.00 3592 3343
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept 1.73 0.04 1.63 1.83 1.00 2627
## premenstrual_phase_fabTRUE -0.01 0.03 -0.10 0.08 1.00 5495
## menstruation 0.04 0.04 -0.06 0.15 1.00 5681
## fertile_fab 0.00 0.09 -0.23 0.25 1.00 4647
## hormonal_contraceptionTRUE 0.07 0.06 -0.09 0.24 1.00 2561
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.02 0.06 -0.12 0.17 1.00 5560
## menstruation:hormonal_contraceptionTRUE -0.00 0.07 -0.18 0.17 1.00 5584
## fertile_fab:hormonal_contraceptionTRUE -0.00 0.15 -0.39 0.38 1.00 4377
## Tail_ESS
## Intercept 2660
## premenstrual_phase_fabTRUE 3281
## menstruation 3384
## fertile_fab 3728
## hormonal_contraceptionTRUE 2896
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3426
## menstruation:hormonal_contraceptionTRUE 3542
## fertile_fab:hormonal_contraceptionTRUE 3472
##
## Family Specific Parameters:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.13 0.01 1.12 1.15 1.00 6175 3055
##
## 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).Don’t adjust for menstruation
dont_adjust = update.formula(formula(model), new = . ~ . - menstruation*hormonal_contraception - premenstrual_phase_fab*hormonal_contraception + hormonal_contraception)
no_mens = update(model, recompile = TRUE, cores = 4, formula = dont_adjust, file=paste0("brms_fits/", model_prefix, "rob_check_no_mens2"))
summary(no_mens)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: time_family ~ fertile_fab + (1 + fertile_fab | person) + hormonal_contraception + fertile_fab:hormonal_contraception
## 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 1737 2747
## sd(fertile_fab) 0.69 0.06 0.55 0.85 1.00 1198 2326
## cor(Intercept,fertile_fab) -0.31 0.07 -0.46 -0.12 1.00 4051 3251
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.72 0.02 1.66 1.78 1.00 1790 2557
## fertile_fab -0.03 0.06 -0.17 0.12 1.00 4775 3375
## hormonal_contraceptionTRUE 0.09 0.05 -0.02 0.21 1.00 1955 2530
## fertile_fab:hormonal_contraceptionTRUE 0.01 0.10 -0.26 0.27 1.00 4657 2926
##
## 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 6258 2421
##
## 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).Don’t adjust for menstruation interaction (bad pre-reg model)
We originally preregistered a model, in which we did not fit an interaction between hormonal contraception and (pre-)menstruation. However, if we want to adjust the effect of the probability of being in the fertile window for (pre-)menstruation, we need to include “interaction controls”. This robustness check shows that this decision makes little difference to the results.
dont_adjust = update.formula(formula(model), new = . ~ . - menstruation:hormonal_contraception - premenstrual_phase_fab:hormonal_contraception)
no_mens_int = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == ""), formula = dont_adjust, file=paste0("brms_fits/", model_prefix, "rob_check_no_mens_int"))
summary(no_mens_int)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: time_family ~ premenstrual_phase_fab + menstruation + fertile_fab + hormonal_contraception + (1 + fertile_fab | person) + fertile_fab:hormonal_contraception
## Data: diary %>% filter(reasons_for_exclusion == "") (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 2244 2645
## sd(fertile_fab) 0.69 0.06 0.53 0.84 1.00 1091 2269
## cor(Intercept,fertile_fab) -0.31 0.07 -0.47 -0.13 1.00 3812 3420
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.70 0.03 1.63 1.77 1.00 2219 2941
## premenstrual_phase_fabTRUE 0.00 0.02 -0.05 0.06 1.00 6523 3473
## menstruation 0.08 0.03 0.02 0.15 1.00 7492 3242
## fertile_fab 0.01 0.06 -0.15 0.16 1.00 5066 3319
## hormonal_contraceptionTRUE 0.09 0.04 -0.03 0.20 1.00 1912 2534
## fertile_fab:hormonal_contraceptionTRUE 0.01 0.10 -0.25 0.26 1.00 4243 3166
##
## 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 4960 2614
##
## 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).Using forward counting
forward_counted = update(model,
recompile = TRUE, cores = 4,
newdata = diary %>%
filter(reasons_for_exclusion == "fertility_never_estimable, next_menstrual_onset_unobserved, " | reasons_for_exclusion == "next_menstrual_onset_unobserved, " | reasons_for_exclusion == "") %>%
mutate(fertile_fab = prc_stirn_b_forward_counted, premenstrual_phase_fab = premenstrual_phase_forward_counted),
file=paste0("brms_fits/", model_prefix, "rob_check_forward_counted2"))
summary(forward_counted)## 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 == "fertili (Number of observations: 26008)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 841)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.50 0.02 0.46 0.54 1.00 1614 2693
## sd(fertile_fab) 0.71 0.06 0.56 0.86 1.01 1330 2257
## cor(Intercept,fertile_fab) -0.35 0.06 -0.50 -0.18 1.00 3554 3190
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept 1.69 0.03 1.62 1.77 1.00 2015
## premenstrual_phase_fabTRUE 0.05 0.03 -0.02 0.11 1.00 5552
## menstruation 0.10 0.03 0.03 0.18 1.00 6492
## fertile_fab 0.04 0.06 -0.12 0.20 1.00 4620
## hormonal_contraceptionTRUE 0.08 0.05 -0.05 0.21 1.00 2150
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.00 0.05 -0.12 0.13 1.00 5123
## menstruation:hormonal_contraceptionTRUE -0.03 0.06 -0.17 0.12 1.00 6617
## fertile_fab:hormonal_contraceptionTRUE -0.02 0.11 -0.29 0.27 1.00 4047
## Tail_ESS
## Intercept 2643
## premenstrual_phase_fabTRUE 3004
## menstruation 3503
## fertile_fab 3450
## hormonal_contraceptionTRUE 2688
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3148
## menstruation:hormonal_contraceptionTRUE 3620
## fertile_fab:hormonal_contraceptionTRUE 3238
##
## 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 5077 2716
##
## 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).Using forward counted narrow window
forward_counted_window = update(model,
recompile = TRUE, cores = 4,
newdata = diary %>%
filter(reasons_for_exclusion == "fertility_never_estimable, next_menstrual_onset_unobserved, " | reasons_for_exclusion == "next_menstrual_onset_unobserved, " | reasons_for_exclusion == "") %>%
mutate(fertile_fab = fertile_narrow_forward_counted),
formula = . ~ . - (premenstrual_phase_fab + menstruation) * hormonal_contraception,
file = paste0("brms_fits/", model_prefix, "rob_check_forward_counted_window"))
summary(forward_counted_window)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: time_family ~ fertile_fab + (1 + fertile_fab | person) + fertile_fab:hormonal_contraception
## Data: diary %>% filter(reasons_for_exclusion == "fertili (Number of observations: 12992)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 825)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.54 0.02 0.48 0.59 1.00 1787 2618
## sd(fertile_fab) 0.83 0.07 0.64 1.02 1.00 906 1347
## cor(Intercept,fertile_fab) -0.44 0.06 -0.59 -0.26 1.00 2210 2824
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.73 0.02 1.67 1.79 1.00 1913 2538
## fertile_fab -0.02 0.06 -0.19 0.15 1.00 3116 2706
## fertile_fab:hormonal_contraceptionTRUE 0.14 0.10 -0.13 0.39 1.00 3527 2748
##
## Family Specific Parameters:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.11 0.01 1.09 1.13 1.00 4224 3054
##
## 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).Using backward counted narrow window
backward_counted_window = update(model,
recompile = TRUE, cores = 4,
newdata = diary %>%
filter(reasons_for_exclusion == "") %>%
mutate(fertile_fab = fertile_narrow),
formula = . ~ . - (premenstrual_phase_fab + menstruation) * hormonal_contraception,
file = paste0("brms_fits/", model_prefix, "rob_check_backward_counted_window"))
summary(backward_counted_window)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: time_family ~ fertile_fab + (1 + fertile_fab | person) + fertile_fab:hormonal_contraception
## Data: diary %>% filter(reasons_for_exclusion == "") %>% (Number of observations: 12743)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 779)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.52 0.02 0.47 0.58 1.01 1554 2672
## sd(fertile_fab) 0.84 0.07 0.63 1.02 1.00 1037 1827
## cor(Intercept,fertile_fab) -0.40 0.07 -0.56 -0.20 1.00 1956 2413
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.74 0.02 1.68 1.80 1.00 1560 2127
## fertile_fab -0.04 0.07 -0.21 0.13 1.00 3011 2941
## fertile_fab:hormonal_contraceptionTRUE 0.24 0.11 -0.04 0.52 1.00 2329 2852
##
## Family Specific Parameters:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.12 0.01 1.10 1.13 1.00 3660 2599
##
## 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).Using backward-counting squished
bc_squished = update(model, recompile = TRUE, cores = 4,newdata = diary %>% filter(reasons_for_exclusion == "") %>% mutate(fertile_fab = prc_stirn_b_squished, premenstrual_phase_fab = premenstrual_phase_squished), file=paste0("brms_fits/", model_prefix, "rob_check_bc_squished"))
summary(bc_squished)## 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 == "") %>% (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 1374 2409
## sd(fertile_fab) 0.70 0.06 0.54 0.85 1.00 1102 1932
## cor(Intercept,fertile_fab) -0.32 0.07 -0.48 -0.14 1.00 2708 3336
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept 1.71 0.03 1.64 1.79 1.01 1527
## premenstrual_phase_fabTRUE -0.02 0.03 -0.09 0.05 1.00 3717
## menstruation 0.08 0.03 -0.00 0.16 1.00 3762
## fertile_fab -0.01 0.07 -0.19 0.15 1.00 2952
## hormonal_contraceptionTRUE 0.06 0.05 -0.08 0.19 1.00 1187
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.06 0.05 -0.06 0.18 1.00 3479
## menstruation:hormonal_contraceptionTRUE 0.01 0.06 -0.15 0.16 1.00 3686
## fertile_fab:hormonal_contraceptionTRUE 0.09 0.12 -0.24 0.40 1.00 2643
## Tail_ESS
## Intercept 1661
## premenstrual_phase_fabTRUE 3586
## menstruation 3151
## fertile_fab 2641
## hormonal_contraceptionTRUE 1175
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3066
## menstruation:hormonal_contraceptionTRUE 3003
## fertile_fab:hormonal_contraceptionTRUE 2600
##
## 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 4351 2675
##
## 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).Using inferred backward-counting
bc_inferred = update(model, recompile = TRUE, cores = 4,newdata = diary %>% filter(reasons_for_exclusion == "fertility_never_estimable, next_menstrual_onset_unobserved, " | reasons_for_exclusion == "next_menstrual_onset_unobserved, " | reasons_for_exclusion == "") %>% mutate(fertile_fab = prc_stirn_b_inferred, premenstrual_phase_fab = premenstrual_phase_inferred), file=paste0("brms_fits/", model_prefix, "rob_check_bc_inferred"))
summary(bc_inferred)## 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 == "fertili (Number of observations: 24667)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 835)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.50 0.02 0.46 0.55 1.00 1587 2545
## sd(fertile_fab) 0.71 0.06 0.56 0.86 1.00 824 1408
## cor(Intercept,fertile_fab) -0.35 0.06 -0.50 -0.18 1.00 2422 3108
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept 1.69 0.03 1.62 1.77 1.00 1218
## premenstrual_phase_fabTRUE 0.05 0.03 -0.02 0.11 1.00 2955
## menstruation 0.11 0.03 0.03 0.18 1.00 3698
## fertile_fab 0.02 0.06 -0.13 0.18 1.00 2057
## hormonal_contraceptionTRUE 0.12 0.05 -0.02 0.26 1.00 1153
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE -0.05 0.05 -0.17 0.08 1.00 2579
## menstruation:hormonal_contraceptionTRUE -0.06 0.06 -0.20 0.08 1.00 3053
## fertile_fab:hormonal_contraceptionTRUE -0.08 0.12 -0.39 0.23 1.00 2077
## Tail_ESS
## Intercept 2120
## premenstrual_phase_fabTRUE 3150
## menstruation 3189
## fertile_fab 2529
## hormonal_contraceptionTRUE 1563
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 2980
## menstruation:hormonal_contraceptionTRUE 3042
## fertile_fab:hormonal_contraceptionTRUE 2408
##
## Family Specific Parameters:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.12 0.01 1.11 1.14 1.00 5508 3063
##
## 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 plot
library(broom.mixed)
robchecks_e <- list(
`A1. Original model` = model,
`E1. No exclusions` = no_exclusions,
`E2. No hypothesis guessing` = hyp_guess,
`E3. No medication` = exclude_meds,
`E4. No cycle-aware` = unaware,
`E5a. No irregular cycles (>± 2 days)` = rob_check_menstruation_regularity,
`E5b. Avg. cycles inside 25-35 days` = rob_check_selfrepcycle_length,
`E5c. No cycles < 25 days` = rob_check_mens_length,
`E5d. Low uncertainty about cycle (<±2 days)` = rob_check_menstruation_certainty,
`E6. No planned pregnancy` = pregnant_trying,
`E7. No unhealthy women` = rob_check_health,
`E8a. Women 18-25` = younger25,
`E8b. Women older than 25` = older25,
`E9a. Friday-Saturday only` = weekend,
`E9b. Monday-Thursday only` = weekday,
`E10. Include living with parents` = rob_check_parents
)
if(outcome != "person_is_related_man_seen") {
robchecks_e[["E11. Singles only"]] = rob_check_singles
}
robchecks_e <- robchecks_e %>%
map(function(x) {
b <- broom::tidy(x, conf.level = 0.99)
b$n_days = nrow(x$data)
b$n_persons = max(n_distinct(x$data$person),
n_distinct(x$data$short))
b$category = "exclusions"
b
}) %>%
bind_rows(.id = "check")
robchecks_e %>%
filter(term == "fertile_fab") %>%
mutate(check = fct_rev(fct_inorder(paste0('atop("",atop("', check,'",',
'phantom(0)^"(',n_persons, ' women, ',n_days,' days)"))')))) %>%
ggplot(aes(check, estimate, ymin = conf.low, ymax = conf.high)) +
geom_hline(yintercept = 0, linetype = 'dashed') +
scale_x_discrete("Robustness check", labels = function(x) {parse(text=x)}) +
ylab("Fertile window effect with 99% CI") +
geom_pointrange() +
coord_flip() +
theme(axis.text.y = element_text(size = 11))
ggsave(paste0("figures/", model_prefix,"_robustness_checks_exclusions.png"), width =5, height = 7)
ggsave(paste0("figures/", model_prefix,"_robustness_checks_exclusions.pdf"), width =5, height = 7)
robchecks_p <- list(
`A1. Original model` = model,
`P1. Not adj. for (pre-)mens.` = no_mens,
`P2. No HC×(pre-)mens. interaction` = no_mens_int,
`P3. Forward counting` = forward_counted,
`P4. Squished backward-counting` = bc_squished,
`P5. Inferred backward-counting` = bc_inferred,
`P6. Forward counted window` = forward_counted_window,
`P7. Backward counted window` = backward_counted_window,
`M1. Full random slopes` = full_random_slopes,
`M2. Uncorrelated random slopes` = full_random_slopes_no_cor,
`M3. No random slopes` = no_varying_slope,
`M4. Require outcome variance` = with_outcome_var
)
if(model$family$family == "gaussian") {
robchecks_p[["M6. Ordinal regression"]] = ordinal_model
} else {
robchecks_p[["M5. Gaussian regression"]] = gaussian_model
}
if(outcome == "person_is_related_man_seen") {
robchecks_p[["O1. Inferred male relatives"]] = related_men_inferred
robchecks_p[["O2. Related men thought about"]] = related_men_thoughts
robchecks_p[["O3. Related men seen/thought about"]] = related_men
}
robchecks_p <- robchecks_p %>%
map(function(x) {
b <- broom::tidy(x, conf.level = 0.99)
b$n_days = nrow(x$data)
b$n_persons = max(n_distinct(x$data$person),
n_distinct(x$data$short))
b$category = "model"
b
}) %>%
bind_rows(.id = "check")
robchecks_p %>%
filter(term == "fertile_fab") %>%
mutate(check = fct_rev(paste0('atop("",atop("', check,'",',
'phantom(0)^"(',n_persons, ' women, ',n_days,' days)"))'))) %>%
ggplot(aes(check, estimate, ymin = conf.low, ymax = conf.high)) +
geom_hline(yintercept = 0, linetype = 'dashed') +
scale_x_discrete("Robustness check", labels = function(x) {parse(text=x)}) +
ylab("Fertile window effect with 99% CI") +
geom_pointrange() +
coord_flip() +
theme(axis.text.y = element_text(size = 11))
ggsave(paste0("figures/", model_prefix,"_robustness_checks_model.png"), width =5, height = 7)
ggsave(paste0("figures/", model_prefix,"_robustness_checks_model.pdf"), width =5, height = 7)
bind_rows(robchecks_e, robchecks_p) %>%
filter(term == "fertile_fab") %>%
mutate(check = fct_rev(fct_inorder(paste0('atop("",atop("', check,'",',
'phantom(0)^"(',n_persons, ' women, ',n_days,' days)"))')))) %>%
ggplot(aes(check, estimate, ymin = conf.low, ymax = conf.high)) +
geom_hline(yintercept = 0, linetype = 'dashed') +
scale_x_discrete("Robustness check", labels = function(x) {parse(text=x)}) +
ylab("Fertile window effect with 99% CI") +
geom_pointrange() +
coord_flip() +
facet_wrap(~ category, scales = "free_y") +
theme(axis.text.y = element_text(size = 11),
strip.text = element_blank())
ggsave(paste0("figures/", model_prefix,"_robustness_checks.png"), width = 7, height = 7)
ggsave(paste0("figures/", model_prefix,"_robustness_checks.pdf"), width = 7, height = 7)related men
model2_relatives <- brm(person_is_related_man_seen ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person), data = diary %>% filter(reasons_for_exclusion == ""), family = poisson(), cores = 4, file ="brms_fits/relatives_prob4")
do_model(model2_relatives, diary %>% filter(!str_detect(reasons_for_exclusion, "not_single, ")), model_prefix = "relatedmen_")Robustness checks
Basic model
library(effects); library(lme4); library(lmerTest); library(knitr); library(brms)
summary(model)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ (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) 3.79 0.46 2.82 5.25 1.00 875 1427
## sd(fertile_fab) 2.29 0.63 1.14 4.61 1.00 1664 2756
## cor(Intercept,fertile_fab) 0.22 0.35 -0.68 0.84 1.00 2325 2311
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept -7.23 0.60 -9.02 -5.92 1.00 1469
## premenstrual_phase_fabTRUE 0.35 0.12 0.04 0.67 1.00 5967
## menstruation 0.40 0.14 0.04 0.74 1.00 7128
## fertile_fab -0.55 1.46 -4.92 2.77 1.00 2437
## hormonal_contraceptionTRUE -0.45 0.90 -2.94 1.67 1.00 890
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.12 0.43 -0.95 1.26 1.00 4894
## menstruation:hormonal_contraceptionTRUE 0.09 0.53 -1.32 1.45 1.00 4742
## fertile_fab:hormonal_contraceptionTRUE -1.70 1.52 -5.64 2.23 1.00 4088
## Tail_ESS
## Intercept 2214
## premenstrual_phase_fabTRUE 3298
## menstruation 3669
## fertile_fab 2622
## hormonal_contraceptionTRUE 1398
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3072
## menstruation:hormonal_contraceptionTRUE 2832
## fertile_fab:hormonal_contraceptionTRUE 3349
##
## 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).outcome <- model$formula[[1]][[2]]Full random slopes
add_mod_formula = update.formula(
formula(model),
new = as.formula(. ~ . - (1 + fertile_fab | person) + (1 + premenstrual_phase_fab + menstruation + fertile_fab | person ))) # reorder so that the triptych looks nice
full_random_slopes = update(model, recompile = TRUE, iter = 6000, warmup = 4000, cores = 4, formula = add_mod_formula, newdata = diary %>% filter(reasons_for_exclusion == ""), file= paste0("brms_fits/", model_prefix, "rob_check_full_random_slopes"))
summary(full_random_slopes)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ premenstrual_phase_fab + menstruation + fertile_fab + hormonal_contraception + (1 + premenstrual_phase_fab + menstruation + fertile_fab | person) + premenstrual_phase_fab:hormonal_contraception + menstruation:hormonal_contraception + fertile_fab:hormonal_contraception
## Data: diary %>% filter(reasons_for_exclusion == "") (Number of observations: 12251)
## Samples: 4 chains, each with iter = 6000; warmup = 4000; thin = 1;
## total post-warmup samples = 8000
##
## Group-Level Effects:
## ~person (Number of levels: 263)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.87 0.47 2.81 5.23 1.00 2007 3726
## sd(premenstrual_phase_fabTRUE) 0.73 0.31 0.02 1.55 1.00 1170 1564
## sd(menstruation) 0.68 0.36 0.01 1.74 1.00 1183 1994
## sd(fertile_fab) 2.28 0.59 1.02 4.28 1.00 2031 3995
## cor(Intercept,premenstrual_phase_fabTRUE) -0.18 0.38 -0.87 0.78 1.00 5124 4706
## cor(Intercept,menstruation) 0.15 0.40 -0.83 0.92 1.00 6882 5222
## cor(premenstrual_phase_fabTRUE,menstruation) 0.35 0.39 -0.83 0.95 1.00 1957 3110
## cor(Intercept,fertile_fab) 0.14 0.33 -0.69 0.79 1.00 5242 5867
## cor(premenstrual_phase_fabTRUE,fertile_fab) 0.17 0.34 -0.81 0.82 1.00 1327 1783
## cor(menstruation,fertile_fab) 0.11 0.36 -0.86 0.84 1.00 1071 2392
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept -7.32 0.64 -9.25 -5.87 1.00 2543
## premenstrual_phase_fabTRUE 0.44 0.48 -0.92 1.92 1.00 4800
## menstruation 0.07 0.54 -1.72 1.34 1.00 5128
## fertile_fab -0.23 1.31 -4.11 2.85 1.00 4516
## hormonal_contraceptionTRUE -0.49 0.92 -2.96 1.73 1.00 1801
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.18 0.54 -1.22 1.64 1.00 5945
## menstruation:hormonal_contraceptionTRUE 0.08 0.63 -1.65 1.71 1.00 6745
## fertile_fab:hormonal_contraceptionTRUE -1.71 1.52 -5.61 2.22 1.00 7001
## Tail_ESS
## Intercept 3935
## premenstrual_phase_fabTRUE 4350
## menstruation 5057
## fertile_fab 5046
## hormonal_contraceptionTRUE 3495
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 6342
## menstruation:hormonal_contraceptionTRUE 6350
## fertile_fab:hormonal_contraceptionTRUE 6741
##
## 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).Full random slopes, no correlations
add_mod_formula = update.formula(
formula(model),
new = as.formula(. ~ . - (1 + fertile_fab | person) + (1 + premenstrual_phase_fab + menstruation + fertile_fab || person ))) # reorder so that the triptych looks nice
full_random_slopes_no_cor = update(model, recompile = TRUE, iter = 6000, warmup = 4000, cores = 4, formula = add_mod_formula, newdata = diary %>% filter(reasons_for_exclusion == ""), file= paste0("brms_fits/", model_prefix, "rob_check_full_random_slopes_no_cor"))
summary(full_random_slopes_no_cor)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ premenstrual_phase_fab + menstruation + fertile_fab + hormonal_contraception + (1 + premenstrual_phase_fab + menstruation + fertile_fab || person) + premenstrual_phase_fab:hormonal_contraception + menstruation:hormonal_contraception + fertile_fab:hormonal_contraception
## Data: diary %>% filter(reasons_for_exclusion == "") (Number of observations: 12251)
## Samples: 4 chains, each with iter = 6000; warmup = 4000; thin = 1;
## total post-warmup samples = 8000
##
## Group-Level Effects:
## ~person (Number of levels: 263)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.85 0.44 2.89 5.21 1.00 2428 3815
## sd(premenstrual_phase_fabTRUE) 0.56 0.24 0.01 1.20 1.00 1650 1876
## sd(menstruation) 0.52 0.27 0.01 1.24 1.00 1790 3535
## sd(fertile_fab) 1.98 0.44 1.00 3.28 1.00 3467 5204
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept -7.33 0.58 -9.12 -6.03 1.00 3513
## premenstrual_phase_fabTRUE 0.25 0.18 -0.27 0.69 1.00 7146
## menstruation 0.28 0.20 -0.33 0.77 1.00 6245
## fertile_fab 0.30 0.50 -1.18 1.44 1.00 5675
## hormonal_contraceptionTRUE -0.46 0.89 -2.93 1.83 1.00 2589
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.16 0.49 -1.10 1.48 1.00 10837
## menstruation:hormonal_contraceptionTRUE 0.04 0.60 -1.54 1.54 1.00 12117
## fertile_fab:hormonal_contraceptionTRUE -1.71 1.48 -5.68 1.94 1.00 10597
## Tail_ESS
## Intercept 5232
## premenstrual_phase_fabTRUE 6147
## menstruation 6138
## fertile_fab 5507
## hormonal_contraceptionTRUE 3924
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 6815
## menstruation:hormonal_contraceptionTRUE 6086
## fertile_fab:hormonal_contraceptionTRUE 6617
##
## 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).Without random slopes
no_varying_slope = update(model, recompile = TRUE, formula = . ~ . - (1 + fertile_fab | person) + (1|person), cores = 4, file = paste0("brms_fits/", model_prefix, "rob_check_no_varying_slope"))
summary(no_varying_slope)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ premenstrual_phase_fab + menstruation + fertile_fab + hormonal_contraception + (1 | person) + premenstrual_phase_fab:hormonal_contraception + menstruation:hormonal_contraception + fertile_fab:hormonal_contraception
## 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) 3.82 0.42 2.92 5.04 1.01 702 1476
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept -7.31 0.55 -9.02 -6.10 1.00 969
## premenstrual_phase_fabTRUE 0.39 0.12 0.10 0.70 1.00 4640
## menstruation 0.43 0.13 0.07 0.77 1.00 4742
## fertile_fab 0.80 0.25 0.14 1.43 1.00 4678
## hormonal_contraceptionTRUE -0.45 0.88 -2.75 1.69 1.02 684
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.13 0.42 -0.89 1.30 1.00 3228
## menstruation:hormonal_contraceptionTRUE 0.06 0.52 -1.24 1.37 1.00 3653
## fertile_fab:hormonal_contraceptionTRUE -2.01 1.16 -5.17 0.92 1.00 3642
## Tail_ESS
## Intercept 2135
## premenstrual_phase_fabTRUE 3701
## menstruation 3127
## fertile_fab 3420
## hormonal_contraceptionTRUE 1747
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3210
## menstruation:hormonal_contraceptionTRUE 3313
## fertile_fab:hormonal_contraceptionTRUE 3184
##
## 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).Only women with outcome variability
with_outcome_var = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "") %>% group_by(session) %>% filter(sd(!!outcome, na.rm = T) > 0), file = paste0("brms_fits/", model_prefix, "rob_check_with_outcome_var"))
summary(with_outcome_var)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person)
## Data: diary %>% filter(reasons_for_exclusion == "") %>% (Number of observations: 2839)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 58)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.88 0.11 0.64 1.21 1.00 1718 2337
## sd(fertile_fab) 1.89 0.40 1.01 3.07 1.00 1335 2246
## cor(Intercept,fertile_fab) -0.27 0.20 -0.69 0.27 1.00 2191 2516
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept -2.05 0.17 -2.51 -1.62 1.00 1609
## premenstrual_phase_fabTRUE 0.35 0.12 0.02 0.67 1.00 4669
## menstruation 0.40 0.14 0.03 0.77 1.00 5591
## fertile_fab 0.58 0.47 -0.77 1.72 1.00 2910
## hormonal_contraceptionTRUE -0.12 0.50 -1.46 1.20 1.00 2179
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.11 0.44 -1.05 1.23 1.00 3898
## menstruation:hormonal_contraceptionTRUE 0.01 0.53 -1.41 1.36 1.00 3863
## fertile_fab:hormonal_contraceptionTRUE -2.00 1.46 -6.17 1.60 1.00 2960
## Tail_ESS
## Intercept 2430
## premenstrual_phase_fabTRUE 3017
## menstruation 3425
## fertile_fab 2602
## hormonal_contraceptionTRUE 2192
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3068
## menstruation:hormonal_contraceptionTRUE 2930
## fertile_fab:hormonal_contraceptionTRUE 2977
##
## 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).Ordinal model
if(identical(model$family$family, "gaussian")) {
ordinal_model = update(model, recompile = TRUE, family = cumulative(), newdata = diary %>% mutate(time_family = as.integer(1 + time_family)) %>% filter(reasons_for_exclusion == ""), file= paste0("brms_fits/", model_prefix, "rob_check_full_ordinal_model"))
summary(ordinal_model)
} Gaussian model
if(identical(model$family$family, "poisson")) {
gaussian_model = update(model, recompile = TRUE, family = gaussian(), file= paste0("brms_fits/", model_prefix, "rob_check_full_ordinal_model"))
summary(gaussian_model)
}## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: person_is_related_man_seen ~ (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) 0.15 0.01 0.13 0.17 1.00 1352 2131
## sd(fertile_fab) 0.21 0.02 0.17 0.26 1.00 1945 2768
## cor(Intercept,fertile_fab) -0.35 0.07 -0.52 -0.15 1.00 3231 3099
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept 0.04 0.01 0.01 0.07 1.00 864
## premenstrual_phase_fabTRUE 0.02 0.01 0.00 0.04 1.00 5775
## menstruation 0.02 0.01 0.00 0.05 1.00 6771
## fertile_fab 0.05 0.02 -0.01 0.11 1.00 2795
## hormonal_contraceptionTRUE -0.01 0.03 -0.08 0.05 1.00 965
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE -0.01 0.02 -0.05 0.04 1.00 5559
## menstruation:hormonal_contraceptionTRUE -0.01 0.02 -0.06 0.04 1.00 5525
## fertile_fab:hormonal_contraceptionTRUE -0.07 0.05 -0.19 0.06 1.00 2369
## Tail_ESS
## Intercept 1432
## premenstrual_phase_fabTRUE 3569
## menstruation 2890
## fertile_fab 3168
## hormonal_contraceptionTRUE 1763
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3330
## menstruation:hormonal_contraceptionTRUE 3643
## fertile_fab:hormonal_contraceptionTRUE 2770
##
## Family Specific Parameters:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.25 0.00 0.25 0.25 1.00 7699 3017
##
## 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).No exclusions
no_exclusions = update(model, recompile = TRUE, cores = 4, newdata = diary, file=paste0("brms_fits/", model_prefix, "rob_check_no_exclusions"))
summary(no_exclusions)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person)
## Data: diary (Number of observations: 15568)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 329)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.81 0.43 2.91 5.11 1.01 683 1232
## sd(fertile_fab) 2.04 0.48 1.08 3.92 1.00 1617 2400
## cor(Intercept,fertile_fab) 0.19 0.34 -0.68 0.82 1.00 2222 2331
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept -7.05 0.54 -8.62 -5.86 1.00 1225
## premenstrual_phase_fabTRUE 0.29 0.10 0.03 0.55 1.00 4288
## menstruation 0.26 0.12 -0.05 0.58 1.00 4422
## fertile_fab -0.37 1.20 -4.02 2.42 1.00 2321
## hormonal_contraceptionTRUE -0.24 0.72 -2.16 1.67 1.00 645
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE -0.25 0.29 -1.01 0.48 1.00 3524
## menstruation:hormonal_contraceptionTRUE 0.07 0.33 -0.80 0.86 1.00 3962
## fertile_fab:hormonal_contraceptionTRUE -1.13 0.96 -3.68 1.33 1.00 2848
## Tail_ESS
## Intercept 2172
## premenstrual_phase_fabTRUE 3538
## menstruation 3216
## fertile_fab 2518
## hormonal_contraceptionTRUE 1317
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3267
## menstruation:hormonal_contraceptionTRUE 2912
## fertile_fab:hormonal_contraceptionTRUE 2743
##
## 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).Exclude cycles shorter than 25 days
rob_check_mens_length = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "") %>% filter(cycle_length >= 25), file=paste0("brms_fits/", model_prefix, "rob_check_mens_length"))
summary(rob_check_mens_length)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person)
## Data: diary %>% filter(reasons_for_exclusion == "") %>% (Number of observations: 10941)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 251)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.62 0.45 2.65 4.97 1.01 742 1577
## sd(fertile_fab) 2.41 0.61 1.23 4.81 1.00 1720 2669
## cor(Intercept,fertile_fab) 0.17 0.34 -0.70 0.81 1.00 2269 2648
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept -6.96 0.58 -8.61 -5.71 1.00 1459
## premenstrual_phase_fabTRUE 0.37 0.13 0.04 0.68 1.00 6078
## menstruation 0.35 0.15 -0.03 0.73 1.00 5860
## fertile_fab -0.45 1.43 -4.88 2.87 1.00 2366
## hormonal_contraceptionTRUE -0.59 0.90 -2.97 1.77 1.00 1196
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.29 0.49 -1.00 1.61 1.00 4715
## menstruation:hormonal_contraceptionTRUE 0.38 0.58 -1.10 1.85 1.00 4537
## fertile_fab:hormonal_contraceptionTRUE -1.13 1.59 -5.24 2.85 1.00 3799
## Tail_ESS
## Intercept 2106
## premenstrual_phase_fabTRUE 3206
## menstruation 3434
## fertile_fab 2658
## hormonal_contraceptionTRUE 2040
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3294
## menstruation:hormonal_contraceptionTRUE 3219
## fertile_fab:hormonal_contraceptionTRUE 3392
##
## 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).Exclude women who report cycles outside 25-35
rob_check_selfrepcycle_length = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "") %>% filter(menstruation_length >= 25, menstruation_length <= 35), file=paste0("brms_fits/", model_prefix, "rob_check_selfrepcycle_length"))
summary(rob_check_selfrepcycle_length)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person)
## Data: diary %>% filter(reasons_for_exclusion == "") %>% (Number of observations: 11413)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 245)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.84 0.49 2.83 5.32 1.00 932 1867
## sd(fertile_fab) 2.25 0.61 1.08 4.68 1.00 1568 2415
## cor(Intercept,fertile_fab) 0.15 0.37 -0.77 0.84 1.00 2460 2511
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept -7.29 0.64 -9.17 -5.87 1.00 1337
## premenstrual_phase_fabTRUE 0.35 0.12 0.03 0.64 1.00 5786
## menstruation 0.39 0.14 0.04 0.75 1.00 6238
## fertile_fab -0.20 1.47 -5.06 3.35 1.00 2446
## hormonal_contraceptionTRUE -0.23 0.90 -2.68 2.14 1.01 835
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.11 0.43 -1.00 1.27 1.00 4238
## menstruation:hormonal_contraceptionTRUE 0.09 0.52 -1.33 1.35 1.00 4684
## fertile_fab:hormonal_contraceptionTRUE -1.84 1.52 -6.10 1.99 1.00 3400
## Tail_ESS
## Intercept 1861
## premenstrual_phase_fabTRUE 3379
## menstruation 3416
## fertile_fab 2203
## hormonal_contraceptionTRUE 1464
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3134
## menstruation:hormonal_contraceptionTRUE 3534
## fertile_fab:hormonal_contraceptionTRUE 3028
##
## 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).Exclude women who report cycles varying more than 1-2 days
rob_check_menstruation_regularity = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "") %>% filter(menstruation_regularity <= 2), file=paste0("brms_fits/", model_prefix, "rob_check_menstruation_regularity"))
summary(rob_check_menstruation_regularity)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person)
## Data: diary %>% filter(reasons_for_exclusion == "") %>% (Number of observations: 7854)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 172)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.89 0.59 2.69 5.67 1.00 986 1944
## sd(fertile_fab) 2.30 0.78 0.78 5.23 1.00 1402 2415
## cor(Intercept,fertile_fab) 0.13 0.40 -0.82 0.86 1.00 2790 3033
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept -7.76 0.81 -10.17 -6.05 1.00 1490
## premenstrual_phase_fabTRUE 0.57 0.16 0.14 1.02 1.00 6257
## menstruation 0.51 0.19 0.02 0.99 1.00 6122
## fertile_fab 0.15 1.76 -5.61 4.36 1.00 2524
## hormonal_contraceptionTRUE 0.33 1.01 -2.27 2.98 1.00 902
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE -0.12 0.44 -1.22 1.06 1.00 4123
## menstruation:hormonal_contraceptionTRUE -0.01 0.54 -1.41 1.41 1.00 4604
## fertile_fab:hormonal_contraceptionTRUE -2.12 1.61 -6.39 2.08 1.00 3469
## Tail_ESS
## Intercept 2155
## premenstrual_phase_fabTRUE 3115
## menstruation 3242
## fertile_fab 2794
## hormonal_contraceptionTRUE 1356
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3148
## menstruation:hormonal_contraceptionTRUE 3283
## fertile_fab:hormonal_contraceptionTRUE 3108
##
## 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).Exclude women who report uncertainty larger than ±2 days
Across three items (last menstrual onset, cycle length, menstruation regularity).
diary$menstruation_certainty <- rowMeans(diary %>% select(menstruation_last_certainty, menstruation_regularity_certainty, menstruation_length_certainty), na.rm = T)
rob_check_menstruation_certainty = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "") %>% filter(menstruation_certainty <= 3), file=paste0("brms_fits/", model_prefix, "rob_check_menstruation_certainty"))
summary(rob_check_menstruation_certainty)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person)
## Data: diary %>% filter(reasons_for_exclusion == "") %>% (Number of observations: 10902)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 233)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.85 0.48 2.83 5.32 1.00 810 1836
## sd(fertile_fab) 2.04 0.55 0.97 4.16 1.00 1533 2381
## cor(Intercept,fertile_fab) 0.10 0.38 -0.80 0.83 1.00 2564 2623
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept -7.18 0.62 -8.85 -5.79 1.00 1470
## premenstrual_phase_fabTRUE 0.32 0.12 0.00 0.63 1.00 6734
## menstruation 0.41 0.14 0.04 0.78 1.00 6371
## fertile_fab 0.16 1.35 -4.11 3.23 1.00 2530
## hormonal_contraceptionTRUE -0.37 0.91 -2.71 1.95 1.00 969
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.16 0.44 -0.96 1.30 1.00 5053
## menstruation:hormonal_contraceptionTRUE 0.08 0.54 -1.41 1.37 1.00 4764
## fertile_fab:hormonal_contraceptionTRUE -1.85 1.47 -5.84 1.64 1.00 4058
## Tail_ESS
## Intercept 2356
## premenstrual_phase_fabTRUE 3348
## menstruation 3304
## fertile_fab 2572
## hormonal_contraceptionTRUE 1883
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3332
## menstruation:hormonal_contraceptionTRUE 3026
## fertile_fab:hormonal_contraceptionTRUE 3209
##
## 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).Exclude women who report not feeling healthy
rob_check_health = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "") %>% filter(self_rated_health <= 2), file=paste0("brms_fits/", model_prefix, "rob_check_health"))
summary(rob_check_health)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person)
## Data: diary %>% filter(reasons_for_exclusion == "") %>% (Number of observations: 11033)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 233)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.76 0.49 2.75 5.26 1.01 995 1941
## sd(fertile_fab) 2.40 0.64 1.19 4.69 1.00 1747 2516
## cor(Intercept,fertile_fab) 0.18 0.35 -0.73 0.83 1.00 2131 2655
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept -7.15 0.63 -8.96 -5.77 1.00 1466
## premenstrual_phase_fabTRUE 0.37 0.12 0.06 0.66 1.00 6209
## menstruation 0.34 0.14 -0.02 0.68 1.00 6402
## fertile_fab -0.52 1.51 -5.06 3.12 1.00 2221
## hormonal_contraceptionTRUE -0.37 0.90 -2.81 1.94 1.00 1024
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.02 0.43 -1.07 1.10 1.00 4486
## menstruation:hormonal_contraceptionTRUE 0.13 0.52 -1.16 1.45 1.00 5262
## fertile_fab:hormonal_contraceptionTRUE -1.69 1.51 -5.76 2.19 1.00 3656
## Tail_ESS
## Intercept 2530
## premenstrual_phase_fabTRUE 3168
## menstruation 3697
## fertile_fab 2793
## hormonal_contraceptionTRUE 1862
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3717
## menstruation:hormonal_contraceptionTRUE 3336
## fertile_fab:hormonal_contraceptionTRUE 3331
##
## 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).Exclude women trying for pregnancy
pregnant_trying = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "") %>% filter(pregnant_trying < 4), file=paste0("brms_fits/", model_prefix, "rob_check_pregnant_trying"))
summary(pregnant_trying)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person)
## Data: diary %>% filter(reasons_for_exclusion == "") %>% (Number of observations: 12096)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 260)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.78 0.46 2.78 5.17 1.01 957 1727
## sd(fertile_fab) 2.33 0.64 1.11 4.64 1.00 2053 2471
## cor(Intercept,fertile_fab) 0.24 0.35 -0.68 0.85 1.00 2754 2657
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept -7.23 0.61 -9.03 -5.86 1.00 1700
## premenstrual_phase_fabTRUE 0.36 0.12 0.05 0.68 1.00 6720
## menstruation 0.39 0.14 0.04 0.74 1.00 6326
## fertile_fab -0.66 1.46 -4.88 2.80 1.00 2781
## hormonal_contraceptionTRUE -0.38 0.86 -2.72 1.70 1.00 1032
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.10 0.43 -0.96 1.24 1.00 5780
## menstruation:hormonal_contraceptionTRUE 0.09 0.52 -1.33 1.44 1.00 5421
## fertile_fab:hormonal_contraceptionTRUE -1.69 1.48 -5.49 2.11 1.00 4440
## Tail_ESS
## Intercept 2177
## premenstrual_phase_fabTRUE 3067
## menstruation 3254
## fertile_fab 2518
## hormonal_contraceptionTRUE 1669
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3452
## menstruation:hormonal_contraceptionTRUE 3239
## fertile_fab:hormonal_contraceptionTRUE 3595
##
## 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).Include women living with parents
rob_check_parents = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "" | reasons_for_exclusion == "living_with_parents, "), file=paste0("brms_fits/", model_prefix, "rob_check_parents"))
summary(rob_check_parents)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person)
## Data: diary %>% filter(reasons_for_exclusion == "" | rea (Number of observations: 13727)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 296)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.84 0.44 2.89 5.22 1.00 881 1900
## sd(fertile_fab) 2.29 0.58 1.21 4.20 1.00 1559 2253
## cor(Intercept,fertile_fab) 0.26 0.34 -0.67 0.85 1.00 2288 2869
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept -7.17 0.57 -8.90 -5.93 1.00 1271
## premenstrual_phase_fabTRUE 0.36 0.11 0.09 0.64 1.00 5735
## menstruation 0.34 0.13 0.01 0.66 1.00 6030
## fertile_fab -0.62 1.37 -4.60 2.62 1.00 2267
## hormonal_contraceptionTRUE -0.10 0.77 -2.20 1.87 1.00 700
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE -0.28 0.31 -1.09 0.51 1.00 6169
## menstruation:hormonal_contraceptionTRUE -0.11 0.37 -1.05 0.86 1.00 6119
## fertile_fab:hormonal_contraceptionTRUE -1.52 1.08 -4.41 1.31 1.00 3408
## Tail_ESS
## Intercept 2235
## premenstrual_phase_fabTRUE 3384
## menstruation 3607
## fertile_fab 2853
## hormonal_contraceptionTRUE 1221
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3623
## menstruation:hormonal_contraceptionTRUE 3284
## fertile_fab:hormonal_contraceptionTRUE 2812
##
## 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).Exclude non-singles
if(outcome != "person_is_related_man_seen") {
rob_check_singles = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "", hetero_relationship == 0, relationship_status == 1), file=paste0("brms_fits/", model_prefix, "rob_check_singles"))
summary(rob_check_singles)
}Moderation by weekend
weekend = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "") %>% filter(weekend == 1), file=paste0("brms_fits/", model_prefix, "rob_check_weekend"))
weekday = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "") %>% filter(weekend == 0), file=paste0("brms_fits/", model_prefix, "rob_check_weekday"))
summary(weekend)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person)
## Data: diary %>% filter(reasons_for_exclusion == "") %>% (Number of observations: 5135)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 262)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.41 0.44 2.51 4.75 1.01 570 1132
## sd(fertile_fab) 1.43 0.65 0.03 3.36 1.00 680 798
## cor(Intercept,fertile_fab) -0.06 0.45 -0.96 0.86 1.00 1743 2238
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept -6.56 0.58 -8.26 -5.27 1.01 1086
## premenstrual_phase_fabTRUE 0.47 0.16 0.03 0.89 1.00 4507
## menstruation 0.56 0.18 0.10 1.00 1.00 4845
## fertile_fab 1.01 1.17 -2.45 4.18 1.00 1618
## hormonal_contraceptionTRUE -1.69 1.06 -4.45 0.77 1.00 1046
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 1.25 0.79 -0.75 3.41 1.00 2369
## menstruation:hormonal_contraceptionTRUE 0.89 0.91 -1.44 3.35 1.00 2693
## fertile_fab:hormonal_contraceptionTRUE 0.28 1.79 -4.27 5.01 1.00 2299
## Tail_ESS
## Intercept 1913
## premenstrual_phase_fabTRUE 2805
## menstruation 3040
## fertile_fab 1702
## hormonal_contraceptionTRUE 1371
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 2754
## menstruation:hormonal_contraceptionTRUE 2796
## fertile_fab:hormonal_contraceptionTRUE 2848
##
## 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(weekday)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person)
## Data: diary %>% filter(reasons_for_exclusion == "") %>% (Number of observations: 7116)
## 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) 3.14 0.42 2.25 4.44 1.00 740 1464
## sd(fertile_fab) 3.33 0.75 1.88 5.98 1.00 1233 1480
## cor(Intercept,fertile_fab) 0.01 0.31 -0.73 0.71 1.00 1526 2274
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept -6.50 0.57 -8.18 -5.20 1.00 1177
## premenstrual_phase_fabTRUE 0.15 0.19 -0.35 0.63 1.00 3644
## menstruation 0.03 0.23 -0.58 0.64 1.00 3901
## fertile_fab -1.53 1.85 -6.86 2.72 1.00 1295
## hormonal_contraceptionTRUE -0.15 0.87 -2.58 2.00 1.00 1079
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE -0.36 0.57 -1.78 1.10 1.00 3079
## menstruation:hormonal_contraceptionTRUE -0.15 0.73 -2.18 1.60 1.00 3240
## fertile_fab:hormonal_contraceptionTRUE -2.84 2.39 -9.24 2.96 1.00 2666
## Tail_ESS
## Intercept 2089
## premenstrual_phase_fabTRUE 3238
## menstruation 3135
## fertile_fab 1771
## hormonal_contraceptionTRUE 1767
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3245
## menstruation:hormonal_contraceptionTRUE 3413
## fertile_fab:hormonal_contraceptionTRUE 2918
##
## 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).Moderation by age group
older25 = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "") %>% filter(age > 25), file=paste0("brms_fits/", model_prefix, "rob_check_older25"))
younger25 = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == "") %>% filter(age <= 25), file=paste0("brms_fits/", model_prefix, "rob_check_younger25"))
summary(younger25)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person)
## Data: diary %>% filter(reasons_for_exclusion == "") %>% (Number of observations: 8380)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 175)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.81 0.55 2.66 5.49 1.00 984 1867
## sd(fertile_fab) 2.31 0.68 1.06 4.82 1.00 1592 2457
## cor(Intercept,fertile_fab) 0.20 0.38 -0.78 0.86 1.00 2359 2519
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept -7.01 0.72 -9.17 -5.41 1.00 1444
## premenstrual_phase_fabTRUE 0.35 0.13 0.02 0.70 1.00 7364
## menstruation 0.43 0.15 0.07 0.81 1.00 6978
## fertile_fab -0.30 1.55 -4.88 3.39 1.00 2394
## hormonal_contraceptionTRUE -1.35 1.17 -4.59 1.74 1.00 1127
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.14 0.63 -1.52 1.84 1.00 4983
## menstruation:hormonal_contraceptionTRUE -0.33 0.79 -2.60 1.55 1.00 4894
## fertile_fab:hormonal_contraceptionTRUE -1.35 2.00 -6.99 3.63 1.00 4235
## Tail_ESS
## Intercept 2321
## premenstrual_phase_fabTRUE 3297
## menstruation 3755
## fertile_fab 2487
## hormonal_contraceptionTRUE 2066
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3458
## menstruation:hormonal_contraceptionTRUE 3352
## fertile_fab:hormonal_contraceptionTRUE 2821
##
## 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(older25)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person)
## Data: diary %>% filter(reasons_for_exclusion == "") %>% (Number of observations: 3871)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 88)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.72 0.78 2.24 6.11 1.01 492 1077
## sd(fertile_fab) 2.43 1.64 0.03 8.21 1.01 267 101
## cor(Intercept,fertile_fab) 0.15 0.51 -0.95 0.98 1.00 1282 2008
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept -7.50 1.06 -10.96 -5.41 1.00 1501
## premenstrual_phase_fabTRUE 0.38 0.29 -0.32 1.10 1.00 3476
## menstruation 0.21 0.35 -0.77 1.12 1.00 3402
## fertile_fab -1.23 2.97 -12.91 4.20 1.01 315
## hormonal_contraceptionTRUE 0.96 1.40 -3.01 4.61 1.01 1213
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.03 0.65 -1.60 1.76 1.00 2794
## menstruation:hormonal_contraceptionTRUE 0.50 0.79 -1.53 2.54 1.01 2983
## fertile_fab:hormonal_contraceptionTRUE -2.21 2.61 -9.90 4.93 1.00 1511
## Tail_ESS
## Intercept 1538
## premenstrual_phase_fabTRUE 2778
## menstruation 2956
## fertile_fab 104
## hormonal_contraceptionTRUE 1958
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 2990
## menstruation:hormonal_contraceptionTRUE 3038
## fertile_fab:hormonal_contraceptionTRUE 2347
##
## 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).Excluding women who used psychopharmacological, hormonal, or antibiotic medication
diary <- diary %>% filter(reasons_for_exclusion == "") %>% ungroup() %>%
mutate(medication_psychopharmacological = if_else(is.na(medication_psychopharmacological), 0, medication_psychopharmacological),
medication_hormonal = if_else(is.na(medication_hormonal), 0, medication_hormonal),
medication_antibiotics = if_else(is.na(medication_antibiotics), 0, medication_antibiotics))
exclude_meds = update(model, recompile = TRUE, cores = 4, newdata= diary %>% filter(reasons_for_exclusion == "") %>% filter(medication_psychopharmacological == 0, medication_hormonal == 0, medication_antibiotics == 0), file=paste0("brms_fits/", model_prefix, "rob_check_meds"))
summary(exclude_meds)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person)
## Data: diary %>% filter(reasons_for_exclusion == "") %>% (Number of observations: 11155)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 238)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.75 0.46 2.77 5.09 1.00 966 1959
## sd(fertile_fab) 2.29 0.60 1.11 4.26 1.00 1760 2535
## cor(Intercept,fertile_fab) 0.18 0.36 -0.74 0.84 1.00 2388 2670
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept -7.14 0.62 -8.93 -5.80 1.00 1312
## premenstrual_phase_fabTRUE 0.44 0.12 0.15 0.75 1.00 7163
## menstruation 0.50 0.14 0.12 0.87 1.00 7094
## fertile_fab -0.16 1.43 -4.42 3.22 1.00 2490
## hormonal_contraceptionTRUE -0.49 0.92 -2.89 1.96 1.00 826
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.02 0.44 -1.03 1.20 1.00 5046
## menstruation:hormonal_contraceptionTRUE -0.01 0.53 -1.37 1.36 1.00 5291
## fertile_fab:hormonal_contraceptionTRUE -1.96 1.53 -6.16 1.72 1.00 3693
## Tail_ESS
## Intercept 2175
## premenstrual_phase_fabTRUE 3505
## menstruation 3791
## fertile_fab 2754
## hormonal_contraceptionTRUE 1556
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3332
## menstruation:hormonal_contraceptionTRUE 3540
## fertile_fab:hormonal_contraceptionTRUE 3381
##
## 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).Excluding hypothesis guessers
hyp_guess = update(model, recompile = TRUE, cores = 4,newdata = diary %>% filter(reasons_for_exclusion == "") %>% filter(hypothesis_guess_topic == "no_guess"), file=paste0("brms_fits/", model_prefix, "_rob_check_hyp_guess"))
summary(hyp_guess)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person)
## Data: diary %>% filter(reasons_for_exclusion == "") %>% (Number of observations: 10707)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 232)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.66 0.47 2.70 5.00 1.00 815 1601
## sd(fertile_fab) 2.45 0.67 1.18 4.83 1.00 1469 2504
## cor(Intercept,fertile_fab) 0.22 0.35 -0.68 0.84 1.00 2395 2585
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept -6.95 0.62 -8.81 -5.65 1.00 1501
## premenstrual_phase_fabTRUE 0.42 0.13 0.09 0.76 1.00 6505
## menstruation 0.37 0.15 -0.02 0.75 1.00 6269
## fertile_fab -0.55 1.51 -5.14 2.67 1.00 2318
## hormonal_contraceptionTRUE -0.51 0.87 -2.89 1.61 1.00 1122
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.04 0.43 -1.05 1.20 1.00 4877
## menstruation:hormonal_contraceptionTRUE 0.12 0.53 -1.27 1.52 1.00 5161
## fertile_fab:hormonal_contraceptionTRUE -1.85 1.55 -6.04 2.07 1.00 3612
## Tail_ESS
## Intercept 2420
## premenstrual_phase_fabTRUE 3537
## menstruation 3118
## fertile_fab 2603
## hormonal_contraceptionTRUE 1863
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3059
## menstruation:hormonal_contraceptionTRUE 3140
## fertile_fab:hormonal_contraceptionTRUE 3163
##
## 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).Excluding women using awareness-based methods
unaware = update(model, recompile = TRUE, cores = 4,newdata = diary %>% filter(reasons_for_exclusion == "") %>% filter(contraception_approach != "awareness", contraception_app == 0), file=paste0("brms_fits/", model_prefix, "rob_check_unaware"))
summary(unaware)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person)
## Data: diary %>% filter(reasons_for_exclusion == "") %>% (Number of observations: 7987)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 177)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 4.06 0.59 2.86 5.87 1.01 858 1673
## sd(fertile_fab) 1.81 0.66 0.53 4.09 1.00 1338 1654
## cor(Intercept,fertile_fab) -0.03 0.45 -0.94 0.85 1.00 2082 2348
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept -8.08 0.81 -10.53 -6.38 1.00 1431
## premenstrual_phase_fabTRUE 0.47 0.15 0.10 0.84 1.00 5256
## menstruation 0.64 0.17 0.18 1.07 1.00 5307
## fertile_fab 0.90 1.62 -4.01 5.05 1.00 2077
## hormonal_contraceptionTRUE 0.68 1.04 -2.17 3.49 1.01 658
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE -0.00 0.44 -1.10 1.11 1.00 3919
## menstruation:hormonal_contraceptionTRUE 0.01 0.53 -1.41 1.43 1.00 4196
## fertile_fab:hormonal_contraceptionTRUE -3.32 1.65 -7.62 0.90 1.00 3519
## Tail_ESS
## Intercept 2091
## premenstrual_phase_fabTRUE 2446
## menstruation 3105
## fertile_fab 2133
## hormonal_contraceptionTRUE 1431
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3252
## menstruation:hormonal_contraceptionTRUE 3217
## fertile_fab:hormonal_contraceptionTRUE 3413
##
## 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).Don’t adjust for menstruation
dont_adjust = update.formula(formula(model), new = . ~ . - menstruation*hormonal_contraception - premenstrual_phase_fab*hormonal_contraception + hormonal_contraception)
no_mens = update(model, recompile = TRUE, cores = 4, formula = dont_adjust, file=paste0("brms_fits/", model_prefix, "rob_check_no_mens2"))
summary(no_mens)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ fertile_fab + (1 + fertile_fab | person) + hormonal_contraception + fertile_fab:hormonal_contraception
## 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) 3.77 0.47 2.75 5.17 1.01 951 1761
## sd(fertile_fab) 2.40 0.65 1.22 4.74 1.00 1399 2046
## cor(Intercept,fertile_fab) 0.20 0.35 -0.72 0.84 1.00 2243 2359
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## Intercept -6.96 0.61 -8.84 -5.60 1.00 1477 2373
## fertile_fab -1.07 1.51 -5.92 2.64 1.00 2207 2508
## hormonal_contraceptionTRUE -0.28 0.81 -2.45 1.88 1.00 1076 1867
## fertile_fab:hormonal_contraceptionTRUE -2.13 1.33 -5.62 1.48 1.00 4537 3551
##
## 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).Don’t adjust for menstruation interaction (bad pre-reg model)
We originally preregistered a model, in which we did not fit an interaction between hormonal contraception and (pre-)menstruation. However, if we want to adjust the effect of the probability of being in the fertile window for (pre-)menstruation, we need to include “interaction controls”. This robustness check shows that this decision makes little difference to the results.
dont_adjust = update.formula(formula(model), new = . ~ . - menstruation:hormonal_contraception - premenstrual_phase_fab:hormonal_contraception)
no_mens_int = update(model, recompile = TRUE, cores = 4, newdata = diary %>% filter(reasons_for_exclusion == ""), formula = dont_adjust, file=paste0("brms_fits/", model_prefix, "rob_check_no_mens_int"))
summary(no_mens_int)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ premenstrual_phase_fab + menstruation + fertile_fab + hormonal_contraception + (1 + fertile_fab | person) + fertile_fab:hormonal_contraception
## 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) 3.77 0.45 2.77 5.11 1.00 701 1294
## sd(fertile_fab) 2.30 0.65 1.08 4.66 1.00 1478 2195
## cor(Intercept,fertile_fab) 0.22 0.36 -0.70 0.85 1.00 2216 2479
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## Intercept -7.22 0.60 -9.00 -5.92 1.01 1453 2538
## premenstrual_phase_fabTRUE 0.36 0.12 0.06 0.68 1.00 7204 3137
## menstruation 0.40 0.14 0.07 0.75 1.00 6554 2955
## fertile_fab -0.55 1.50 -5.03 2.94 1.00 2206 2119
## hormonal_contraceptionTRUE -0.35 0.82 -2.49 1.79 1.00 839 1365
## fertile_fab:hormonal_contraceptionTRUE -1.96 1.28 -5.37 1.33 1.00 4551 3132
##
## 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).Using forward counting
forward_counted = update(model,
recompile = TRUE, cores = 4,
newdata = diary %>%
filter(reasons_for_exclusion == "fertility_never_estimable, next_menstrual_onset_unobserved, " | reasons_for_exclusion == "next_menstrual_onset_unobserved, " | reasons_for_exclusion == "") %>%
mutate(fertile_fab = prc_stirn_b_forward_counted, premenstrual_phase_fab = premenstrual_phase_forward_counted),
file=paste0("brms_fits/", model_prefix, "rob_check_forward_counted2"))
summary(forward_counted)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person)
## Data: diary %>% filter(reasons_for_exclusion == "fertili (Number of observations: 12895)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 275)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.84 0.46 2.87 5.25 1.01 686 1656
## sd(fertile_fab) 2.74 0.58 1.66 4.96 1.00 1698 2086
## cor(Intercept,fertile_fab) 0.12 0.34 -0.67 0.80 1.00 1869 2293
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept -7.00 0.60 -8.86 -5.66 1.00 1458
## premenstrual_phase_fabTRUE -0.04 0.12 -0.33 0.28 1.00 5844
## menstruation 0.07 0.13 -0.28 0.41 1.00 5847
## fertile_fab -1.88 1.58 -6.86 1.69 1.00 2046
## hormonal_contraceptionTRUE -0.95 0.83 -3.25 1.10 1.01 891
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.69 0.40 -0.32 1.77 1.00 4083
## menstruation:hormonal_contraceptionTRUE 0.52 0.44 -0.61 1.71 1.00 3975
## fertile_fab:hormonal_contraceptionTRUE -0.20 1.41 -3.87 3.40 1.00 2577
## Tail_ESS
## Intercept 2049
## premenstrual_phase_fabTRUE 2660
## menstruation 3143
## fertile_fab 2117
## hormonal_contraceptionTRUE 1395
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3393
## menstruation:hormonal_contraceptionTRUE 3242
## fertile_fab:hormonal_contraceptionTRUE 3109
##
## 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).Using forward counted narrow window
forward_counted_window = update(model,
recompile = TRUE, cores = 4,
newdata = diary %>%
filter(reasons_for_exclusion == "fertility_never_estimable, next_menstrual_onset_unobserved, " | reasons_for_exclusion == "next_menstrual_onset_unobserved, " | reasons_for_exclusion == "") %>%
mutate(fertile_fab = fertile_narrow_forward_counted),
formula = . ~ . - (premenstrual_phase_fab + menstruation) * hormonal_contraception,
file = paste0("brms_fits/", model_prefix, "rob_check_forward_counted_window"))
summary(forward_counted_window)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ fertile_fab + (1 + fertile_fab | person) + fertile_fab:hormonal_contraception
## Data: diary %>% filter(reasons_for_exclusion == "fertili (Number of observations: 6293)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 271)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.20 0.46 2.22 4.63 1.00 546 1083
## sd(fertile_fab) 3.46 0.76 1.96 6.20 1.00 1152 1674
## cor(Intercept,fertile_fab) 0.07 0.30 -0.65 0.71 1.00 1268 1458
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## Intercept -6.38 0.61 -8.23 -5.04 1.00 1100 1430
## fertile_fab -2.00 1.80 -7.65 2.10 1.00 1265 1722
## fertile_fab:hormonal_contraceptionTRUE -1.20 1.56 -5.50 2.41 1.00 1919 2729
##
## 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).Using backward counted narrow window
backward_counted_window = update(model,
recompile = TRUE, cores = 4,
newdata = diary %>%
filter(reasons_for_exclusion == "") %>%
mutate(fertile_fab = fertile_narrow),
formula = . ~ . - (premenstrual_phase_fab + menstruation) * hormonal_contraception,
file = paste0("brms_fits/", model_prefix, "rob_check_backward_counted_window"))
summary(backward_counted_window)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ fertile_fab + (1 + fertile_fab | person) + fertile_fab:hormonal_contraception
## Data: diary %>% filter(reasons_for_exclusion == "") %>% (Number of observations: 6272)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 262)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 2.96 0.38 2.13 4.09 1.00 689 1230
## sd(fertile_fab) 2.26 0.60 1.07 4.24 1.00 961 1754
## cor(Intercept,fertile_fab) 0.08 0.33 -0.69 0.75 1.00 1560 2275
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## Intercept -6.20 0.52 -7.81 -5.05 1.00 1230 1759
## fertile_fab -0.28 1.30 -3.93 2.78 1.00 1443 2213
## fertile_fab:hormonal_contraceptionTRUE -1.38 1.37 -5.16 2.06 1.00 2956 2849
##
## 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).Using backward-counting squished
bc_squished = update(model, recompile = TRUE, cores = 4,newdata = diary %>% filter(reasons_for_exclusion == "") %>% mutate(fertile_fab = prc_stirn_b_squished, premenstrual_phase_fab = premenstrual_phase_squished), file=paste0("brms_fits/", model_prefix, "rob_check_bc_squished"))
summary(bc_squished)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ (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) 3.81 0.46 2.81 5.16 1.00 758 1465
## sd(fertile_fab) 2.41 0.58 1.22 4.31 1.00 1311 2587
## cor(Intercept,fertile_fab) 0.14 0.35 -0.72 0.79 1.00 1831 2178
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept -7.30 0.61 -8.92 -5.97 1.00 1198
## premenstrual_phase_fabTRUE 0.37 0.12 0.06 0.70 1.00 4610
## menstruation 0.42 0.15 0.04 0.80 1.00 4073
## fertile_fab -0.30 1.44 -4.39 3.37 1.00 2069
## hormonal_contraceptionTRUE -0.36 0.88 -2.73 1.96 1.00 790
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.09 0.44 -1.02 1.29 1.00 4241
## menstruation:hormonal_contraceptionTRUE 0.04 0.54 -1.31 1.46 1.00 4048
## fertile_fab:hormonal_contraceptionTRUE -1.76 1.54 -5.89 2.05 1.00 3210
## Tail_ESS
## Intercept 2021
## premenstrual_phase_fabTRUE 3080
## menstruation 2811
## fertile_fab 2368
## hormonal_contraceptionTRUE 1265
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3443
## menstruation:hormonal_contraceptionTRUE 3377
## fertile_fab:hormonal_contraceptionTRUE 3469
##
## 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).Using inferred backward-counting
bc_inferred = update(model, recompile = TRUE, cores = 4,newdata = diary %>% filter(reasons_for_exclusion == "fertility_never_estimable, next_menstrual_onset_unobserved, " | reasons_for_exclusion == "next_menstrual_onset_unobserved, " | reasons_for_exclusion == "") %>% mutate(fertile_fab = prc_stirn_b_inferred, premenstrual_phase_fab = premenstrual_phase_inferred), file=paste0("brms_fits/", model_prefix, "rob_check_bc_inferred"))
summary(bc_inferred)## Family: poisson
## Links: mu = log
## Formula: person_is_related_man_seen ~ (premenstrual_phase_fab + menstruation + fertile_fab) * hormonal_contraception + (1 + fertile_fab | person)
## Data: diary %>% filter(reasons_for_exclusion == "fertili (Number of observations: 12136)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~person (Number of levels: 271)
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.83 0.47 2.77 5.18 1.00 953 1879
## sd(fertile_fab) 2.65 0.56 1.54 4.63 1.00 1940 2548
## cor(Intercept,fertile_fab) 0.08 0.34 -0.73 0.79 1.00 2320 2613
##
## Population-Level Effects:
## Estimate Est.Error l-99% CI u-99% CI Rhat Bulk_ESS
## Intercept -7.12 0.62 -9.02 -5.72 1.00 1605
## premenstrual_phase_fabTRUE 0.13 0.13 -0.20 0.46 1.00 6838
## menstruation 0.23 0.14 -0.16 0.59 1.00 6713
## fertile_fab -1.24 1.51 -5.60 2.41 1.00 2419
## hormonal_contraceptionTRUE -0.77 0.87 -3.14 1.46 1.00 940
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 0.50 0.44 -0.55 1.64 1.00 4973
## menstruation:hormonal_contraceptionTRUE 0.34 0.47 -0.86 1.58 1.00 4650
## fertile_fab:hormonal_contraceptionTRUE -0.68 1.44 -4.66 3.05 1.00 3642
## Tail_ESS
## Intercept 2184
## premenstrual_phase_fabTRUE 3421
## menstruation 3027
## fertile_fab 2499
## hormonal_contraceptionTRUE 1702
## premenstrual_phase_fabTRUE:hormonal_contraceptionTRUE 3425
## menstruation:hormonal_contraceptionTRUE 3271
## fertile_fab:hormonal_contraceptionTRUE 3190
##
## 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 plot
library(broom.mixed)
robchecks_e <- list(
`A1. Original model` = model,
`E1. No exclusions` = no_exclusions,
`E2. No hypothesis guessing` = hyp_guess,
`E3. No medication` = exclude_meds,
`E4. No cycle-aware` = unaware,
`E5a. No irregular cycles (>± 2 days)` = rob_check_menstruation_regularity,
`E5b. Avg. cycles inside 25-35 days` = rob_check_selfrepcycle_length,
`E5c. No cycles < 25 days` = rob_check_mens_length,
`E5d. Low uncertainty about cycle (<±2 days)` = rob_check_menstruation_certainty,
`E6. No planned pregnancy` = pregnant_trying,
`E7. No unhealthy women` = rob_check_health,
`E8a. Women 18-25` = younger25,
`E8b. Women older than 25` = older25,
`E9a. Friday-Saturday only` = weekend,
`E9b. Monday-Thursday only` = weekday,
`E10. Include living with parents` = rob_check_parents
)
if(outcome != "person_is_related_man_seen") {
robchecks_e[["E11. Singles only"]] = rob_check_singles
}
robchecks_e <- robchecks_e %>%
map(function(x) {
b <- broom::tidy(x, conf.level = 0.99)
b$n_days = nrow(x$data)
b$n_persons = max(n_distinct(x$data$person),
n_distinct(x$data$short))
b$category = "exclusions"
b
}) %>%
bind_rows(.id = "check")
robchecks_e %>%
filter(term == "fertile_fab") %>%
mutate(check = fct_rev(fct_inorder(paste0('atop("",atop("', check,'",',
'phantom(0)^"(',n_persons, ' women, ',n_days,' days)"))')))) %>%
ggplot(aes(check, estimate, ymin = conf.low, ymax = conf.high)) +
geom_hline(yintercept = 0, linetype = 'dashed') +
scale_x_discrete("Robustness check", labels = function(x) {parse(text=x)}) +
ylab("Fertile window effect with 99% CI") +
geom_pointrange() +
coord_flip() +
theme(axis.text.y = element_text(size = 11))
ggsave(paste0("figures/", model_prefix,"_robustness_checks_exclusions.png"), width =5, height = 7)
ggsave(paste0("figures/", model_prefix,"_robustness_checks_exclusions.pdf"), width =5, height = 7)
robchecks_p <- list(
`A1. Original model` = model,
`P1. Not adj. for (pre-)mens.` = no_mens,
`P2. No HC×(pre-)mens. interaction` = no_mens_int,
`P3. Forward counting` = forward_counted,
`P4. Squished backward-counting` = bc_squished,
`P5. Inferred backward-counting` = bc_inferred,
`P6. Forward counted window` = forward_counted_window,
`P7. Backward counted window` = backward_counted_window,
`M1. Full random slopes` = full_random_slopes,
`M2. Uncorrelated random slopes` = full_random_slopes_no_cor,
`M3. No random slopes` = no_varying_slope,
`M4. Require outcome variance` = with_outcome_var
)
if(model$family$family == "gaussian") {
robchecks_p[["M6. Ordinal regression"]] = ordinal_model
} else {
robchecks_p[["M5. Gaussian regression"]] = gaussian_model
}
if(outcome == "person_is_related_man_seen") {
robchecks_p[["O1. Inferred male relatives"]] = related_men_inferred
robchecks_p[["O2. Related men thought about"]] = related_men_thoughts
robchecks_p[["O3. Related men seen/thought about"]] = related_men
}
robchecks_p <- robchecks_p %>%
map(function(x) {
b <- broom::tidy(x, conf.level = 0.99)
b$n_days = nrow(x$data)
b$n_persons = max(n_distinct(x$data$person),
n_distinct(x$data$short))
b$category = "model"
b
}) %>%
bind_rows(.id = "check")
robchecks_p %>%
filter(term == "fertile_fab") %>%
mutate(check = fct_rev(paste0('atop("",atop("', check,'",',
'phantom(0)^"(',n_persons, ' women, ',n_days,' days)"))'))) %>%
ggplot(aes(check, estimate, ymin = conf.low, ymax = conf.high)) +
geom_hline(yintercept = 0, linetype = 'dashed') +
scale_x_discrete("Robustness check", labels = function(x) {parse(text=x)}) +
ylab("Fertile window effect with 99% CI") +
geom_pointrange() +
coord_flip() +
theme(axis.text.y = element_text(size = 11))
ggsave(paste0("figures/", model_prefix,"_robustness_checks_model.png"), width =5, height = 7)
ggsave(paste0("figures/", model_prefix,"_robustness_checks_model.pdf"), width =5, height = 7)
bind_rows(robchecks_e, robchecks_p) %>%
filter(term == "fertile_fab") %>%
mutate(check = fct_rev(fct_inorder(paste0('atop("",atop("', check,'",',
'phantom(0)^"(',n_persons, ' women, ',n_days,' days)"))')))) %>%
ggplot(aes(check, estimate, ymin = conf.low, ymax = conf.high)) +
geom_hline(yintercept = 0, linetype = 'dashed') +
scale_x_discrete("Robustness check", labels = function(x) {parse(text=x)}) +
ylab("Fertile window effect with 99% CI") +
geom_pointrange() +
coord_flip() +
facet_wrap(~ category, scales = "free_y") +
theme(axis.text.y = element_text(size = 11),
strip.text = element_blank())
ggsave(paste0("figures/", model_prefix,"_robustness_checks.png"), width = 7, height = 7)
ggsave(paste0("figures/", model_prefix,"_robustness_checks.pdf"), width = 7, height = 7)