knitr::opts_chunk$set(fig.width = 9, fig.height = 9, message = F, warning = F, echo = TRUE, error = TRUE, results = "markup")
library(brms);
library(bayesplot)
library(tidybayes)
options(width=9000)
plot_varying_slopes <- function(model) {
ranefs <- posterior_samples(model, pars = "^(b_fertile|^r_woman.+fertile)", as.matrix = TRUE)
model_data <- model$data
pars <- ranefs %>%
spread_draws(b_fertile,
`b_fertile:nohc1`,
r_woman[woman, slope])
pars2 <- pars %>%
right_join(model_data %>% select(woman, nohc) %>% group_by(nohc, woman) %>% summarise(days = n()), by = "woman") %>%
mutate(woman_slope = b_fertile +
if_else(nohc == "1", `b_fertile:nohc1`, 0) +
r_woman) %>%
select(.chain, .iteration, .draw, days, nohc, woman_slope) %>%
group_by(days, nohc, woman) %>%
mean_hdi(woman_slope)
pars2 <- pars2 %>%
ungroup() %>%
arrange(woman_slope) %>%
mutate(woman = forcats::fct_inorder(woman))
ggplot(pars2 %>% filter(days > 50)) +
geom_hline(yintercept = 0, linetype = 'dashed', size = 1, color = "gray") +
geom_pointrange(aes(woman, woman_slope, ymin = .lower, ymax = .upper, color = nohc)) +
coord_flip() +
theme(axis.text.y = element_blank(), axis.ticks.y = element_blank()) +
scale_color_manual("Contraception", labels = c("0" = "HC user", "1" = "Cycling"), values = c("0" = 'black', '1' = 'red'), breaks = c("0", "1")) +
facet_wrap(~ nohc, scales = "free_y", labeller = labeller(nohc = c("0" = "HC user", "1" = "Cycling"))) +
scale_y_continuous("Fertile phase slope on general desire") +
ggtitle("Varying slopes of fertile window effect", paste("for women with more than 50 days in diary","\nModel fit in",sprintf("%.1f", model_time(model)),"hours"))
}
libido_gaussian <- readRDS("routine_and_sex/ras/libido_gaussian.rds")
summary(libido_gaussian, priors = TRUE)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: high_libido ~ (menstruation + premenstruation + fertile + weekend) * nohc + (1 + fertile | person | gr(woman, by = nohc))
## Data: desire_wide %>% filter(first_row) (Number of observations: 38254)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Priors:
## b ~ normal(0,1)
## Intercept ~ student_t(3, 2, 10)
## L ~ lkj_corr_cholesky(1)
## sd ~ normal(0,1)
## sigma ~ student_t(3, 0, 10)
##
## Group-Level Effects:
## ~woman (Number of levels: 872)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept:nohc0) 0.55 0.03 0.50 0.61 1.00 842 1692
## sd(fertile:nohc0) 0.64 0.08 0.49 0.80 1.01 807 1582
## sd(Intercept:nohc1) 0.59 0.02 0.55 0.63 1.01 660 1383
## sd(fertile:nohc1) 0.86 0.05 0.77 0.96 1.00 1063 2180
## cor(Intercept:nohc0,fertile:nohc0) 0.07 0.11 -0.14 0.30 1.00 1830 2449
## cor(Intercept:nohc1,fertile:nohc1) -0.08 0.06 -0.20 0.04 1.00 1413 2047
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 2.16 0.04 2.09 2.24 1.01 686 1462
## menstruation -0.16 0.04 -0.23 -0.09 1.00 1887 2454
## premenstruationTRUE -0.07 0.03 -0.13 -0.02 1.00 2135 2962
## fertile 0.06 0.07 -0.09 0.20 1.00 1849 2470
## weekend 0.17 0.02 0.13 0.21 1.00 3157 3049
## nohc1 -0.06 0.05 -0.15 0.04 1.01 638 1294
## menstruation:nohc1 -0.01 0.04 -0.09 0.08 1.00 1960 2421
## premenstruationTRUE:nohc1 -0.04 0.03 -0.11 0.02 1.00 2095 2613
## fertile:nohc1 0.44 0.09 0.26 0.62 1.00 1517 2200
## weekend:nohc1 -0.06 0.02 -0.10 -0.01 1.00 3084 2871
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.05 0.00 1.04 1.06 1.00 6828 2835
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot_varying_slopes(libido_gaussian)
libido_gaussian_sigma_nohc <- readRDS("routine_and_sex/ras/libido_gaussian_sigma_nohc.rds")
summary(libido_gaussian_sigma_nohc, priors = TRUE)
## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: high_libido ~ (menstruation + premenstruation + fertile + weekend) * nohc + (1 + fertile | person | gr(woman, by = nohc))
## sigma ~ 0 + nohc
## Data: desire_wide %>% filter(first_row) (Number of observations: 38254)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Priors:
## b ~ normal(0,1)
## Intercept ~ student_t(3, 2, 10)
## L ~ lkj_corr_cholesky(1)
## sd ~ normal(0,1)
##
## Group-Level Effects:
## ~woman (Number of levels: 872)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept:nohc0) 0.56 0.03 0.50 0.61 1.00 2040 3511
## sd(fertile:nohc0) 0.70 0.07 0.55 0.84 1.00 2338 3939
## sd(Intercept:nohc1) 0.59 0.02 0.55 0.63 1.00 1887 3633
## sd(fertile:nohc1) 0.86 0.05 0.77 0.97 1.00 2923 4696
## cor(Intercept:nohc0,fertile:nohc0) 0.03 0.10 -0.16 0.24 1.00 4683 5261
## cor(Intercept:nohc1,fertile:nohc1) -0.08 0.06 -0.21 0.04 1.00 4218 5148
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 2.16 0.04 2.08 2.24 1.00 1683 2972
## menstruation -0.16 0.03 -0.23 -0.09 1.00 5981 5810
## premenstruationTRUE -0.07 0.03 -0.13 -0.02 1.00 5329 5582
## fertile 0.06 0.07 -0.08 0.21 1.00 4279 5337
## weekend 0.17 0.02 0.14 0.21 1.00 9063 6483
## nohc1 -0.05 0.05 -0.15 0.04 1.00 1639 3097
## menstruation:nohc1 -0.01 0.04 -0.09 0.07 1.00 6069 6005
## premenstruationTRUE:nohc1 -0.04 0.03 -0.11 0.02 1.00 5154 6286
## fertile:nohc1 0.44 0.09 0.25 0.62 1.00 4205 5626
## weekend:nohc1 -0.06 0.02 -0.10 -0.02 1.00 8830 6575
## sigma_nohc1 0.05 0.00 0.04 0.06 1.00 14420 5859
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot_varying_slopes(libido_gaussian_sigma_nohc)
libido_gaussian_sigma_limited <- readRDS("routine_and_sex/ras/libido_gaussian_sigma_limited.rds")
summary(libido_gaussian_sigma_limited, priors = TRUE)
## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: high_libido ~ (menstruation + premenstruation + fertile + weekend) * nohc + (1 + fertile | gr(woman, by = nohc))
## sigma ~ 0 + nohc + (1 | woman)
## Data: desire_wide %>% filter(first_row) (Number of observations: 38254)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Priors:
## b ~ normal(0,1)
## b_sigma ~ normal(0,1)
## Intercept ~ student_t(3, 2, 10)
## L ~ lkj_corr_cholesky(1)
## sd ~ normal(0,1)
## sd_sigma ~ normal(0,1)
##
## Group-Level Effects:
## ~woman (Number of levels: 872)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept:nohc0) 0.42 0.25 0.00 0.61 3.30 4 13
## sd(fertile:nohc0) 0.53 0.59 0.01 1.50 4.41 4 11
## sd(Intercept:nohc1) 0.55 0.05 0.48 0.62 3.48 4 18
## sd(fertile:nohc1) 0.13 0.13 0.00 0.29 3.39 4 15
## sd(sigma_Intercept) 0.69 0.09 0.61 0.84 3.90 4 12
## cor(Intercept:nohc0,fertile:nohc0) -0.01 0.61 -0.79 0.92 3.73 4 11
## cor(Intercept:nohc1,fertile:nohc1) 0.17 0.85 -1.00 1.00 3.26 4 11
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 2.13 0.13 1.91 2.26 3.00 5 20
## menstruation -0.05 0.08 -0.19 0.00 2.19 6 11
## premenstruationTRUE -0.03 0.06 -0.13 0.00 2.07 6 28
## fertile 0.27 0.30 -0.00 0.78 3.06 5 24
## weekend 0.03 0.06 -0.00 0.14 1.89 7 11
## nohc1 -0.01 0.15 -0.13 0.25 4.23 4 11
## menstruation:nohc1 0.05 0.08 -0.00 0.19 2.19 6 11
## premenstruationTRUE:nohc1 0.03 0.06 -0.00 0.13 2.07 6 27
## fertile:nohc1 0.03 0.53 -0.78 0.61 2.96 5 24
## weekend:nohc1 -0.03 0.06 -0.14 0.00 1.89 7 11
## sigma_nohc1 -0.04 0.03 -0.08 -0.01 3.50 4 14
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot_varying_slopes(libido_gaussian_sigma_limited)
libido_gaussian_sigma_limitedinits0 <- readRDS("routine_and_sex/ras/libido_gaussian_sigma_limitedinits0.rds")
summary(libido_gaussian_sigma_limitedinits0, priors = TRUE)
## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: high_libido ~ (menstruation + premenstruation + fertile + weekend) * nohc + (1 + fertile | gr(woman, by = nohc))
## sigma ~ 0 + nohc + (1 | woman)
## Data: desire_wide %>% filter(first_row) (Number of observations: 38254)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Priors:
## b ~ normal(0,1)
## b_sigma ~ normal(0,1)
## Intercept ~ student_t(3, 2, 10)
## L ~ lkj_corr_cholesky(1)
## sd ~ normal(0,1)
## sd_sigma ~ normal(0,1)
##
## Group-Level Effects:
## ~woman (Number of levels: 872)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept:nohc0) 0.74 0.06 0.64 0.79 3.29 4 15
## sd(fertile:nohc0) 0.29 0.22 0.01 0.61 4.07 4 11
## sd(Intercept:nohc1) 0.77 0.09 0.68 0.92 4.19 4 12
## sd(fertile:nohc1) 0.65 0.19 0.36 0.86 3.88 4 12
## sd(sigma_Intercept) 1.03 0.13 0.81 1.16 3.69 4 11
## cor(Intercept:nohc0,fertile:nohc0) 0.49 0.26 0.09 0.78 3.50 4 11
## cor(Intercept:nohc1,fertile:nohc1) 0.35 0.26 0.06 0.74 3.38 4 11
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 2.15 0.08 2.05 2.26 3.38 4 16
## menstruation -0.00 0.00 -0.00 0.00 2.52 5 12
## premenstruationTRUE 0.00 0.00 -0.00 0.00 1.95 6 14
## fertile 0.17 0.09 0.01 0.25 3.94 4 12
## weekend -0.00 0.00 -0.00 0.00 2.83 5 18
## nohc1 -0.10 0.10 -0.25 0.03 3.41 4 13
## menstruation:nohc1 0.00 0.00 -0.00 0.00 2.53 5 12
## premenstruationTRUE:nohc1 -0.00 0.00 -0.00 0.00 1.95 6 14
## fertile:nohc1 0.48 0.12 0.38 0.68 3.27 4 16
## weekend:nohc1 0.00 0.00 -0.00 0.00 2.83 5 18
## sigma_nohc1 -0.05 0.02 -0.09 -0.04 3.70 4 11
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot_varying_slopes(libido_gaussian_sigma_limitedinits0)
plot(libido_gaussian_sigma_limitedinits0)
libido_gaussian_sigma_full <- readRDS("routine_and_sex/ras/libido_gaussian_sigma_full.rds")
summary(libido_gaussian_sigma_full, priors = TRUE)
## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: high_libido ~ (menstruation + premenstruation + fertile + weekend) * nohc + (1 + fertile | person | gr(woman, by = nohc))
## sigma ~ (menstruation + premenstruation + fertile + weekend) * nohc + (1 + fertile | person | gr(woman, by = nohc))
## Data: desire_wide %>% filter(first_row) (Number of observations: 38254)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Priors:
## b ~ normal(0,1)
## Intercept ~ student_t(3, 2, 10)
## Intercept_sigma ~ student_t(3, 0, 10)
## L ~ lkj_corr_cholesky(1)
## sd ~ normal(0,1)
## sd_sigma ~ student_t(3, 0, 10)
##
## Group-Level Effects:
## ~woman (Number of levels: 872)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept:nohc0) 0.45 0.08 0.34 0.57 3.38 4 12
## sd(fertile:nohc0) 0.15 0.10 0.06 0.32 3.56 4 15
## sd(sigma_Intercept:nohc0) 0.39 0.25 0.21 0.82 3.08 5 26
## sd(sigma_fertile:nohc0) 0.06 0.03 0.03 0.10 2.97 5 19
## sd(Intercept:nohc1) 0.34 0.21 0.00 0.56 3.49 4 11
## sd(fertile:nohc1) 0.55 0.59 0.17 1.59 3.49 4 11
## sd(sigma_Intercept:nohc1) 0.63 0.27 0.19 0.92 4.11 4 11
## sd(sigma_fertile:nohc1) 0.65 0.32 0.14 1.03 3.83 4 11
## cor(Intercept:nohc0,fertile:nohc0) 0.44 0.49 -0.18 1.00 3.34 4 11
## cor(Intercept:nohc0,sigma_Intercept:nohc0) 0.65 0.09 0.50 0.74 3.18 5 19
## cor(fertile:nohc0,sigma_Intercept:nohc0) 0.51 0.08 0.39 0.63 3.41 4 15
## cor(Intercept:nohc0,sigma_fertile:nohc0) 0.13 0.26 -0.11 0.56 3.22 4 11
## cor(fertile:nohc0,sigma_fertile:nohc0) -0.07 0.51 -0.59 0.76 3.11 4 15
## cor(sigma_Intercept:nohc0,sigma_fertile:nohc0) 0.03 0.34 -0.41 0.52 3.24 4 16
## cor(Intercept:nohc1,fertile:nohc1) 0.31 0.79 -0.93 1.00 4.26 4 11
## cor(Intercept:nohc1,sigma_Intercept:nohc1) 0.15 0.59 -0.85 0.63 3.22 4 11
## cor(fertile:nohc1,sigma_Intercept:nohc1) 0.51 0.17 0.33 0.79 3.67 4 11
## cor(Intercept:nohc1,sigma_fertile:nohc1) -0.49 0.19 -0.75 -0.29 3.01 5 16
## cor(fertile:nohc1,sigma_fertile:nohc1) -0.19 0.34 -0.58 0.34 3.09 5 16
## cor(sigma_Intercept:nohc1,sigma_fertile:nohc1) -0.59 0.63 -0.99 0.51 2.48 5 11
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 2.13 0.02 2.11 2.18 2.96 5 14
## sigma_Intercept -0.00 0.03 -0.04 0.05 3.56 4 12
## menstruation -0.07 0.04 -0.11 0.00 3.33 4 12
## premenstruationTRUE -0.04 0.02 -0.06 0.00 3.19 4 12
## fertile 0.07 0.06 -0.03 0.11 3.39 4 14
## weekend 0.12 0.07 -0.00 0.19 2.90 5 23
## nohc1 -0.08 0.08 -0.22 -0.02 3.90 4 13
## menstruation:nohc1 0.03 0.11 -0.17 0.11 3.69 4 12
## premenstruationTRUE:nohc1 0.01 0.07 -0.11 0.06 3.63 4 11
## fertile:nohc1 0.63 0.14 0.44 0.85 3.79 4 11
## weekend:nohc1 -0.10 0.11 -0.19 0.09 3.14 5 16
## sigma_menstruation -0.05 0.03 -0.10 -0.02 2.63 5 13
## sigma_premenstruationTRUE -0.01 0.02 -0.03 0.03 3.60 4 11
## sigma_fertile -0.03 0.06 -0.12 0.06 3.14 4 12
## sigma_weekend 0.10 0.01 0.08 0.12 3.52 4 11
## sigma_nohc1 0.05 0.07 -0.02 0.18 2.89 5 20
## sigma_menstruation:nohc1 0.01 0.04 -0.03 0.08 3.87 4 11
## sigma_premenstruationTRUE:nohc1 -0.01 0.04 -0.07 0.03 3.05 5 14
## sigma_fertile:nohc1 0.06 0.18 -0.22 0.27 4.00 4 11
## sigma_weekend:nohc1 -0.07 0.02 -0.10 -0.04 3.78 4 11
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot_varying_slopes(libido_gaussian_sigma_full)
libido_cumulative <- readRDS("routine_and_sex/ras/libido_cumulative.rds")
summary(libido_cumulative, priors = TRUE)
## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: high_libido ~ (menstruation + premenstruation + fertile + weekend) * nohc + (1 + fertile | person | gr(woman, by = nohc))
## Data: desire_wide %>% filter(first_row) (Number of observations: 38254)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Priors:
## b ~ normal(0,1)
## Intercept ~ student_t(3, 0, 10)
## L ~ lkj_corr_cholesky(1)
## sd ~ normal(0,1)
##
## Group-Level Effects:
## ~woman (Number of levels: 872)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept:nohc0) 1.05 0.06 0.94 1.16 1.00 843 1725
## sd(fertile:nohc0) 1.06 0.14 0.77 1.34 1.00 1046 1794
## sd(Intercept:nohc1) 1.22 0.04 1.14 1.31 1.00 773 1575
## sd(fertile:nohc1) 1.53 0.10 1.34 1.73 1.00 1061 1609
## cor(Intercept:nohc0,fertile:nohc0) 0.08 0.13 -0.17 0.34 1.00 2030 2146
## cor(Intercept:nohc1,fertile:nohc1) -0.26 0.07 -0.38 -0.12 1.00 1742 2379
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -0.39 0.07 -0.53 -0.24 1.00 593 1438
## Intercept[2] 0.66 0.07 0.52 0.81 1.00 592 1404
## Intercept[3] 1.97 0.08 1.83 2.12 1.00 602 1407
## Intercept[4] 3.62 0.08 3.47 3.78 1.00 637 1459
## menstruation -0.27 0.06 -0.40 -0.15 1.00 2258 2789
## premenstruationTRUE -0.13 0.05 -0.24 -0.03 1.00 2126 2733
## fertile 0.06 0.13 -0.19 0.32 1.00 1942 2517
## weekend 0.31 0.04 0.23 0.37 1.00 3002 2712
## nohc1 -0.15 0.10 -0.35 0.03 1.01 510 1009
## menstruation:nohc1 -0.04 0.08 -0.19 0.11 1.00 2138 2572
## premenstruationTRUE:nohc1 -0.09 0.07 -0.22 0.03 1.00 2081 2963
## fertile:nohc1 0.85 0.16 0.52 1.17 1.00 1900 2355
## weekend:nohc1 -0.10 0.04 -0.18 -0.01 1.00 2951 2717
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot_varying_slopes(libido_cumulative)
libido_cumulative_disc <- readRDS("routine_and_sex/ras/libido_cumulative_disc.rds")
summary(libido_cumulative_disc, priors = TRUE)
## Family: cumulative
## Links: mu = logit; disc = log
## Formula: high_libido ~ (menstruation + premenstruation + fertile + weekend) * nohc + (1 + fertile | person | gr(woman, by = nohc))
## disc ~ (menstruation + premenstruation + fertile + weekend) * nohc + (1 + fertile | person | gr(woman, by = nohc))
## Data: desire_wide %>% filter(first_row) (Number of observations: 38254)
## Samples: 3 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 3000
##
## Priors:
## b ~ normal(0,1)
## Intercept ~ student_t(3, 0, 10)
## Intercept_disc ~ normal(0, 1)
## L ~ lkj_corr_cholesky(1)
## sd ~ normal(0,1)
## sd_disc ~ student_t(3, 0, 10)
##
## Group-Level Effects:
## ~woman (Number of levels: 872)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept:nohc0) 2.38 0.44 1.57 3.29 1.03 269 644
## sd(fertile:nohc0) 2.04 0.43 1.26 2.95 1.01 378 971
## sd(disc_Intercept:nohc0) 0.40 0.02 0.36 0.45 1.00 737 1467
## sd(disc_fertile:nohc0) 0.11 0.09 0.00 0.32 1.01 506 735
## sd(Intercept:nohc1) 2.94 0.55 1.95 4.09 1.03 242 545
## sd(fertile:nohc1) 3.23 0.63 2.08 4.54 1.03 288 572
## sd(disc_Intercept:nohc1) 0.41 0.02 0.38 0.45 1.00 1021 1473
## sd(disc_fertile:nohc1) 0.60 0.07 0.47 0.73 1.00 462 823
## cor(Intercept:nohc0,fertile:nohc0) 0.13 0.14 -0.15 0.42 1.00 1498 1431
## cor(Intercept:nohc0,disc_Intercept:nohc0) 0.04 0.08 -0.12 0.19 1.00 459 1003
## cor(fertile:nohc0,disc_Intercept:nohc0) -0.12 0.13 -0.37 0.13 1.01 204 372
## cor(Intercept:nohc0,disc_fertile:nohc0) 0.05 0.39 -0.70 0.79 1.00 2621 1895
## cor(fertile:nohc0,disc_fertile:nohc0) -0.03 0.44 -0.82 0.80 1.00 1678 2073
## cor(disc_Intercept:nohc0,disc_fertile:nohc0) -0.09 0.39 -0.77 0.70 1.00 2481 2204
## cor(Intercept:nohc1,fertile:nohc1) -0.27 0.08 -0.42 -0.12 1.00 1189 1852
## cor(Intercept:nohc1,disc_Intercept:nohc1) 0.21 0.06 0.10 0.33 1.01 629 885
## cor(fertile:nohc1,disc_Intercept:nohc1) -0.12 0.08 -0.28 0.04 1.02 211 439
## cor(Intercept:nohc1,disc_fertile:nohc1) 0.03 0.10 -0.16 0.21 1.00 1581 2141
## cor(fertile:nohc1,disc_fertile:nohc1) -0.27 0.12 -0.52 -0.03 1.00 594 1449
## cor(disc_Intercept:nohc1,disc_fertile:nohc1) -0.25 0.09 -0.42 -0.06 1.00 1316 2274
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -0.85 0.22 -1.32 -0.45 1.01 335 633
## Intercept[2] 1.56 0.34 0.94 2.27 1.03 249 542
## Intercept[3] 4.63 0.87 3.03 6.42 1.03 237 518
## Intercept[4] 9.13 1.69 6.07 12.58 1.03 248 523
## disc_Intercept -0.81 0.19 -1.15 -0.41 1.03 247 570
## menstruation -0.53 0.16 -0.88 -0.27 1.01 646 1025
## premenstruationTRUE -0.24 0.12 -0.50 -0.02 1.00 955 1209
## fertile 0.45 0.31 -0.11 1.10 1.01 690 1421
## weekend 0.55 0.12 0.34 0.81 1.02 355 857
## nohc1 -0.46 0.22 -0.93 -0.06 1.01 369 640
## menstruation:nohc1 -0.09 0.16 -0.45 0.22 1.01 1268 1345
## premenstruationTRUE:nohc1 -0.27 0.15 -0.60 -0.00 1.01 812 972
## fertile:nohc1 1.92 0.46 1.08 2.89 1.01 492 1105
## weekend:nohc1 -0.09 0.09 -0.28 0.09 1.00 1952 2149
## disc_menstruation 0.06 0.04 -0.03 0.13 1.00 1265 2241
## disc_premenstruationTRUE 0.02 0.03 -0.05 0.09 1.00 1114 1585
## disc_fertile 0.10 0.08 -0.06 0.25 1.00 1334 1920
## disc_weekend -0.09 0.02 -0.14 -0.04 1.00 1850 2003
## disc_nohc1 -0.07 0.05 -0.16 0.02 1.01 679 1031
## disc_menstruation:nohc1 -0.02 0.05 -0.11 0.08 1.00 1237 1803
## disc_premenstruationTRUE:nohc1 0.01 0.04 -0.07 0.08 1.00 1166 1620
## disc_fertile:nohc1 0.03 0.09 -0.16 0.21 1.00 1250 1691
## disc_weekend:nohc1 0.07 0.03 0.02 0.13 1.00 1858 1930
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot_varying_slopes(libido_cumulative_disc)
libido_cumulative_disc_limited <- readRDS("routine_and_sex/ras/libido_cumulative_disc_limited.rds")
summary(libido_cumulative_disc_limited, priors = TRUE)
## Family: cumulative
## Links: mu = logit; disc = log
## Formula: high_libido ~ (menstruation + premenstruation + fertile + weekend) * nohc + (1 + fertile | person | gr(woman, by = nohc))
## disc ~ nohc + (1 | person | gr(woman, by = nohc))
## Data: desire_wide %>% filter(first_row) (Number of observations: 38254)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Priors:
## b ~ normal(0,1)
## Intercept ~ student_t(3, 0, 10)
## Intercept_disc ~ normal(0, 1)
## L ~ lkj_corr_cholesky(1)
## sd ~ normal(0,1)
## sd_disc ~ student_t(3, 0, 10)
##
## Group-Level Effects:
## ~woman (Number of levels: 872)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept:nohc0) 2.35 0.44 1.56 3.27 1.01 546 971
## sd(fertile:nohc0) 2.05 0.43 1.30 2.94 1.01 652 1215
## sd(disc_Intercept:nohc0) 0.40 0.02 0.36 0.44 1.00 1234 2358
## sd(Intercept:nohc1) 2.92 0.54 1.95 4.01 1.01 472 1024
## sd(fertile:nohc1) 3.31 0.64 2.16 4.66 1.01 542 966
## sd(disc_Intercept:nohc1) 0.39 0.02 0.36 0.42 1.00 1053 2052
## cor(Intercept:nohc0,fertile:nohc0) 0.12 0.14 -0.14 0.39 1.00 2092 2805
## cor(Intercept:nohc0,disc_Intercept:nohc0) 0.04 0.08 -0.10 0.20 1.00 834 1740
## cor(fertile:nohc0,disc_Intercept:nohc0) -0.11 0.13 -0.35 0.15 1.04 175 410
## cor(Intercept:nohc1,fertile:nohc1) -0.26 0.07 -0.39 -0.12 1.00 1944 2526
## cor(Intercept:nohc1,disc_Intercept:nohc1) 0.23 0.05 0.13 0.34 1.01 1027 1953
## cor(fertile:nohc1,disc_Intercept:nohc1) -0.17 0.07 -0.32 -0.03 1.02 304 636
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -0.84 0.23 -1.35 -0.45 1.01 479 1062
## Intercept[2] 1.55 0.34 0.97 2.27 1.01 555 1117
## Intercept[3] 4.58 0.86 3.05 6.32 1.01 525 1024
## Intercept[4] 8.95 1.64 6.00 12.31 1.01 513 989
## disc_Intercept -0.80 0.19 -1.15 -0.41 1.01 523 1010
## menstruation -0.56 0.16 -0.90 -0.28 1.00 1102 1905
## premenstruationTRUE -0.25 0.11 -0.48 -0.04 1.00 1728 1856
## fertile 0.38 0.30 -0.16 1.02 1.01 1035 1466
## weekend 0.60 0.13 0.37 0.87 1.01 641 1659
## nohc1 -0.43 0.22 -0.89 -0.02 1.00 507 1094
## menstruation:nohc1 -0.07 0.16 -0.39 0.22 1.00 2495 2375
## premenstruationTRUE:nohc1 -0.27 0.14 -0.59 -0.01 1.01 1446 2007
## fertile:nohc1 1.79 0.44 0.99 2.70 1.00 852 1982
## weekend:nohc1 -0.14 0.09 -0.33 0.03 1.00 2700 2433
## disc_nohc1 -0.04 0.03 -0.10 0.03 1.00 1038 1426
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot_varying_slopes(libido_cumulative_disc_limited)
libido_cumulative_disc_limited <- readRDS("routine_and_sex/ras/libido_cumulative_disccmc.rds")
## Error in gzfile(file, "rb"): cannot open the connection
summary(libido_cumulative_disc_limited, priors = TRUE)
## Family: cumulative
## Links: mu = logit; disc = log
## Formula: high_libido ~ (menstruation + premenstruation + fertile + weekend) * nohc + (1 + fertile | person | gr(woman, by = nohc))
## disc ~ nohc + (1 | person | gr(woman, by = nohc))
## Data: desire_wide %>% filter(first_row) (Number of observations: 38254)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Priors:
## b ~ normal(0,1)
## Intercept ~ student_t(3, 0, 10)
## Intercept_disc ~ normal(0, 1)
## L ~ lkj_corr_cholesky(1)
## sd ~ normal(0,1)
## sd_disc ~ student_t(3, 0, 10)
##
## Group-Level Effects:
## ~woman (Number of levels: 872)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept:nohc0) 2.35 0.44 1.56 3.27 1.01 546 971
## sd(fertile:nohc0) 2.05 0.43 1.30 2.94 1.01 652 1215
## sd(disc_Intercept:nohc0) 0.40 0.02 0.36 0.44 1.00 1234 2358
## sd(Intercept:nohc1) 2.92 0.54 1.95 4.01 1.01 472 1024
## sd(fertile:nohc1) 3.31 0.64 2.16 4.66 1.01 542 966
## sd(disc_Intercept:nohc1) 0.39 0.02 0.36 0.42 1.00 1053 2052
## cor(Intercept:nohc0,fertile:nohc0) 0.12 0.14 -0.14 0.39 1.00 2092 2805
## cor(Intercept:nohc0,disc_Intercept:nohc0) 0.04 0.08 -0.10 0.20 1.00 834 1740
## cor(fertile:nohc0,disc_Intercept:nohc0) -0.11 0.13 -0.35 0.15 1.04 175 410
## cor(Intercept:nohc1,fertile:nohc1) -0.26 0.07 -0.39 -0.12 1.00 1944 2526
## cor(Intercept:nohc1,disc_Intercept:nohc1) 0.23 0.05 0.13 0.34 1.01 1027 1953
## cor(fertile:nohc1,disc_Intercept:nohc1) -0.17 0.07 -0.32 -0.03 1.02 304 636
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -0.84 0.23 -1.35 -0.45 1.01 479 1062
## Intercept[2] 1.55 0.34 0.97 2.27 1.01 555 1117
## Intercept[3] 4.58 0.86 3.05 6.32 1.01 525 1024
## Intercept[4] 8.95 1.64 6.00 12.31 1.01 513 989
## disc_Intercept -0.80 0.19 -1.15 -0.41 1.01 523 1010
## menstruation -0.56 0.16 -0.90 -0.28 1.00 1102 1905
## premenstruationTRUE -0.25 0.11 -0.48 -0.04 1.00 1728 1856
## fertile 0.38 0.30 -0.16 1.02 1.01 1035 1466
## weekend 0.60 0.13 0.37 0.87 1.01 641 1659
## nohc1 -0.43 0.22 -0.89 -0.02 1.00 507 1094
## menstruation:nohc1 -0.07 0.16 -0.39 0.22 1.00 2495 2375
## premenstruationTRUE:nohc1 -0.27 0.14 -0.59 -0.01 1.01 1446 2007
## fertile:nohc1 1.79 0.44 0.99 2.70 1.00 852 1982
## weekend:nohc1 -0.14 0.09 -0.33 0.03 1.00 2700 2433
## disc_nohc1 -0.04 0.03 -0.10 0.03 1.00 1038 1426
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot_varying_slopes(libido_cumulative_disc_limited)
libido_cumulative_disc_limited <- readRDS("routine_and_sex/ras/libido_cumulative_disccmc_limited.rds")
summary(libido_cumulative_disc_limited, priors = TRUE)
## Family: cumulative
## Links: mu = logit; disc = log
## Formula: high_libido ~ (menstruation + premenstruation + fertile + weekend) * nohc + (1 + fertile | person | gr(woman, by = nohc))
## disc ~ 0 + nohc + (1 | person | gr(woman, by = nohc))
## Data: desire_wide %>% filter(first_row) (Number of observations: 38254)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Priors:
## b ~ normal(0,1)
## Intercept ~ student_t(3, 0, 10)
## L ~ lkj_corr_cholesky(1)
## sd ~ normal(0,1)
## sd_disc ~ student_t(3, 0, 10)
##
## Group-Level Effects:
## ~woman (Number of levels: 872)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept:nohc0) 1.06 0.07 0.94 1.20 1.00 1483 2250
## sd(fertile:nohc0) 1.00 0.14 0.71 1.28 1.00 1253 2103
## sd(disc_Intercept:nohc0) 0.40 0.02 0.35 0.44 1.00 2027 2823
## sd(Intercept:nohc1) 1.30 0.06 1.18 1.43 1.00 1369 2232
## sd(fertile:nohc1) 1.48 0.11 1.28 1.70 1.00 1632 2572
## sd(disc_Intercept:nohc1) 0.39 0.02 0.36 0.43 1.00 1535 2412
## cor(Intercept:nohc0,fertile:nohc0) 0.11 0.14 -0.16 0.38 1.00 3703 3354
## cor(Intercept:nohc0,disc_Intercept:nohc0) 0.04 0.08 -0.12 0.19 1.00 1295 2186
## cor(fertile:nohc0,disc_Intercept:nohc0) -0.07 0.12 -0.31 0.17 1.04 99 716
## cor(Intercept:nohc1,fertile:nohc1) -0.26 0.07 -0.39 -0.12 1.00 3490 3244
## cor(Intercept:nohc1,disc_Intercept:nohc1) 0.24 0.05 0.14 0.34 1.00 1720 2761
## cor(fertile:nohc1,disc_Intercept:nohc1) -0.18 0.07 -0.32 -0.04 1.01 544 1201
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -0.39 0.08 -0.55 -0.24 1.00 1011 1723
## Intercept[2] 0.67 0.08 0.52 0.83 1.00 1071 1855
## Intercept[3] 2.02 0.09 1.85 2.22 1.00 1209 2109
## Intercept[4] 3.97 0.13 3.72 4.24 1.00 1396 1911
## menstruation -0.26 0.06 -0.38 -0.14 1.00 4940 3064
## premenstruationTRUE -0.12 0.05 -0.21 -0.02 1.00 5130 2992
## fertile 0.11 0.13 -0.16 0.37 1.00 2004 2922
## weekend 0.27 0.03 0.20 0.33 1.00 5722 3184
## nohc1 -0.22 0.10 -0.42 -0.01 1.01 889 1517
## menstruation:nohc1 -0.02 0.07 -0.16 0.13 1.00 6215 3219
## premenstruationTRUE:nohc1 -0.11 0.06 -0.23 0.01 1.00 5417 3802
## fertile:nohc1 0.90 0.17 0.57 1.24 1.00 2199 2397
## weekend:nohc1 -0.06 0.04 -0.14 0.01 1.00 6768 3247
## disc_nohc1 -0.05 0.03 -0.12 0.02 1.00 1494 2155
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot_varying_slopes(libido_cumulative_disc_limited)
Model comparison
loo_compare(libido_cumulative, libido_cumulative_disc, libido_cumulative_disc_limited)
## elpd_diff se_diff
## libido_cumulative_disc 0.0 0.0
## libido_cumulative_disc_limited -39.2 203.8
## libido_cumulative -1493.7 54.9
Varying slope comparison
model <- libido_cumulative
ranefs <- posterior_samples(model, pars = "^(b_fertile|^r_woman.+fertile)", as.matrix = TRUE)
model_data <- model$data
pars <- ranefs %>%
spread_draws(b_fertile,
`b_fertile:nohc1`,
r_woman[woman, slope])
pars2 <- pars %>%
right_join(model_data %>% select(woman, nohc) %>% group_by(nohc, woman) %>% summarise(days = n()), by = "woman") %>%
mutate(woman_slope = b_fertile +
if_else(nohc == "1", `b_fertile:nohc1`, 0) +
r_woman) %>%
select(.chain, .iteration, .draw, days, nohc, woman_slope) %>%
group_by(days, nohc, woman) %>%
mean_hdi(woman_slope)
libido_cumulative_slopes <- pars2 %>%
ungroup() %>%
arrange(woman_slope) %>%
mutate(woman = forcats::fct_inorder(woman))
model <- libido_cumulative_disc_limited
ranefs <- posterior_samples(model, pars = "^(b_fertile|^r_woman.+fertile)", as.matrix = TRUE)
model_data <- model$data
pars <- ranefs %>%
spread_draws(b_fertile,
`b_fertile:nohc1`,
r_woman[woman, slope])
pars2 <- pars %>%
right_join(model_data %>% select(woman, nohc) %>% group_by(nohc, woman) %>% summarise(days = n()), by = "woman") %>%
mutate(woman_slope = b_fertile +
if_else(nohc == "1", `b_fertile:nohc1`, 0) +
r_woman) %>%
select(.chain, .iteration, .draw, days, nohc, woman_slope) %>%
group_by(days, nohc, woman) %>%
mean_hdi(woman_slope)
libido_cumulative_disc_slopes <- pars2 %>%
ungroup() %>%
arrange(woman_slope) %>%
mutate(woman = forcats::fct_inorder(woman))
libido_cumulative_slopes %>% full_join(libido_cumulative_disc_slopes, by = c("woman", "nohc", "days"), suffix = c("_mlm", "_mlsem")) %>%
arrange(woman_slope_mlm) %>%
mutate(woman = forcats::fct_inorder(woman)) %>%
filter(days > 50) %>%
ggplot(.) +
geom_hline(yintercept = 0, linetype = 'dashed', size = 1, color = "gray") +
geom_pointrange(aes(woman, woman_slope_mlm, ymin = .lower_mlm, ymax = .upper_mlm), color = "#689A65", linetype = 'dashed', shape = 3, alpha = 0.7, position = position_dodge(width = 0.2)) +
geom_pointrange(aes(woman, woman_slope_mlsem, ymin = .lower_mlsem, ymax = .upper_mlsem), color = "#626B9D", position = position_dodge(width = 0.2), alpha = 0.7) +
coord_flip() +
theme(axis.text.y = element_blank(), axis.ticks.y = element_blank()) +
scale_color_manual("Contraception", labels = c("0" = "HC user", "1" = "Cycling"), values = c("0" = 'black', '1' = 'red'), breaks = c("0", "1")) +
facet_wrap(~ nohc, scales = "free_y", labeller = labeller(nohc = c("0" = "HC user", "1" = "Cycling"))) +
scale_y_continuous("Fertile phase slope on general desire") +
ggtitle("Varying slopes of fertile window effect", paste("for women with more than 50 days in diary","\nModel fit in",sprintf("%.1f", model_time(libido_gaussian)),"hours"))
Outcome distribution
ggplot(mapping = aes(x = libido_gaussian$data[,1])) + geom_bar()
