When a cell is empty, this reflects no difference to M_1.
To see predictor definitions, see this table (Table 3 in the manuscript).
Stricter exclusion criteria (further down), subsume the ones above.
See also 0_helpers.R
robustness_check_ovu_shift = function(obj, diary) {
library(lme4); library(lmerTest); library(sjPlot)
get_coefs = function(fit, model_name) {
library(broom)
if (class(fit) == "lme") {
obj_coef = tryCatch({my_tidy.lme(fit, effects = 'fixed') }, error = function(e) { cat_message(e, "danger") })
} else {
obj_coef = tryCatch({tidy(fit, conf.int = TRUE, effects = "fixed")}, error = function(e) { cat_message(e, "danger") })
}
obj_coef$model = model_name
obj_coef
}
outcome = as.character(formula(obj)[2])
outcome_label = recode(str_replace_all(str_replace_all(str_replace_all(outcome, "_", " "), " pair", "-pair"), " 1", ""),
"desirability" = "self-perceived desirability",
"NARQ admiration" = "narcissistic admiration",
"NARQ rivalry" = "narcissistic rivalry",
"extra-pair" = "extra-pair desire & behaviour",
"had sexual intercourse" = "sexual intercourse")
less_ctrl_formula = update.formula(formula(obj), new = . ~ . - included * menstruation - fertile_mean)
tryCatch({
update(obj, formula = less_ctrl_formula) -> fewer_controls
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
one_ctrl_formula = update.formula(formula(obj), new = . ~ . - fertile_mean)
update(obj, formula = one_ctrl_formula) -> dontcontrolavg
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
diary_nododgy = diary %>% filter(dodgy_data == F)
update(obj, data = diary_nododgy) -> nododgy
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
diary_nopreg = diary %>% filter(trying_to_get_pregnant == 'no')
update(obj, data = diary_nopreg) -> notrypreg
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
diary_broad = diary %>% mutate(fertile = fertile_broad)
update(obj, formula = less_ctrl_formula, data = diary_broad) -> broad_window
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
diary_broad = diary %>% mutate(fertile = fertile_broad)
update(obj, data = diary_broad) -> broad_window_ctrl
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
diary_bci = diary %>% mutate(fertile = prc_stirn_b_backward_inferred)
update(obj, formula = less_ctrl_formula, data = diary_bci) -> backward_inferred
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
diary_narrow = diary %>% mutate(fertile = fertile_narrow)
update(obj, formula = less_ctrl_formula, data = diary_narrow) -> narrow_window
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
diary_cont_fc = diary %>% mutate(fertile = prc_stirn_b_forward_counted)
update(obj, data = diary_cont_fc) -> forward_counting
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
diary_narrow_fc = diary %>% mutate(fertile = fertile_narrow_forward_counted)
update(obj, formula = less_ctrl_formula, data = diary_narrow_fc) -> forward_counting_narrow
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
diary_broad_fc = diary %>% mutate(fertile = fertile_broad_forward_counted)
update(obj, formula = less_ctrl_formula, data = diary_broad_fc) -> forward_counting_broad
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
diary_cont = diary %>% mutate(fertile = prc_stirn_b)
update(obj, data = diary_cont) -> backward
}, error = function(e) { cat_message(e, "danger") })
# tryCatch({
# diary %>% mutate(fertile = prc_wcx_b) %>%
# update(obj, data = .) -> wilcox
# }, error = function(e) { cat_message(e, "danger") })
tryCatch({
ilax = diary %>% mutate(included = included_lax)
update(obj, data = ilax) -> included_lax
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
if(outcome != 'self_esteem_1') {
update(obj, formula = update.formula(formula(obj), new = . ~ . + self_esteem_1)) -> control_self_esteem
}
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
iconser = diary %>% mutate(included = included_conservative)
update(obj, data = iconser) -> included_conservative
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
istrict = diary %>% mutate(included = included_strict)
update(obj, data = istrict) -> included_strict
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
no_short = diary %>% filter(is.na(cycle_length_diary) | cycle_length_diary >= 20)
update(obj, data = no_short) -> no_cycles_shorter_than_20
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
no_long = diary %>% filter(is.na(cycle_length_diary) | (cycle_length_diary >= 20 & cycle_length_diary <= 40))
update(obj, data = no_long) -> cycles_between_20_40
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
update(obj, formula = update.formula(formula(obj), new = . ~ . + weekday + week_number)) -> control_week
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
update(obj, formula = update.formula(formula(obj), new = . ~ . + time_of_response + log10(time_for_response+1))) -> control_time_of_response
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
diary %>%
mutate(fertile = prc_stirn_b_forward_counted) %>%
filter(!is.na(fertile) & !is.na(included) & menstruation == "no") %>%
filter_(lazyeval::interp( ~ !is.na(outcome), outcome = as.name(obj@frame %>% names() %>% .[1]))) %>%
filter(!duplicated(person)) -> diary_tmp2
form = update.formula(formula(obj), new = . ~ . - fertile_mean - menstruation - included:menstruation - (1 | person))
environment(form) = environment()
lm(formula = form, data = diary_tmp2) -> between_person
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
four_occ = diary %>%
mutate(fertile = fertile_narrow_forward_counted,
person_occasion = paste(person, fertile)) %>%
filter(!is.na(fertile) & !is.na(included) & menstruation == "no") %>%
filter_(lazyeval::interp( ~ !is.na(outcome), outcome = as.name(obj@frame %>% names() %>% .[1]))) %>%
mutate(row_nr_even = row_number(person) %% 2,
person_occasion = paste(person_occasion, row_nr_even)) %>%
filter(!duplicated(person_occasion)) %>%
group_by(person) %>%
mutate(nr_days = n()) %>%
filter(nr_days == 4)
update(obj, formula = update.formula(formula(obj), new = . ~ . - fertile_mean - included:menstruation - menstruation) , data = four_occ) -> within_person_four_occasions
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
avg_twotime = diary %>%
mutate(fertile = fertile_broad) %>%
filter(!is.na(fertile) & !is.na(included) & menstruation == "no") %>%
select(person, fertile, included, one_of(outcome)) %>%
group_by(included, person, fertile) %>%
summarise_all(funs(mean(.,na.rm = T)))
update(obj, formula = less_ctrl_formula , data = avg_twotime) -> within_person_avg_hilow
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
two_occ = diary %>%
mutate(fertile = fertile_narrow_forward_counted,
person_occasion = paste(person, fertile)) %>%
filter(!is.na(fertile) & !is.na(included) & menstruation == "no") %>%
filter_(lazyeval::interp( ~ !is.na(outcome), outcome = as.name(obj@frame %>% names() %>% .[1]))) %>%
filter(!duplicated(person_occasion)) %>%
group_by(person) %>%
mutate(nr_days = n()) %>%
filter(nr_days == 2)
update(obj, formula = update.formula(formula(obj), new = . ~ . - fertile_mean - included:menstruation - menstruation) , data = two_occ) -> within_person_two_occasions
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
# findbars(formula(obj))[[1]]
reactive = diary %>% group_by(person) %>% mutate(days_filled_out = row_number(day_number)) %>% filter(day_number < 40, days_filled_out < 37)
new_form = update.formula(nobars(formula(obj)), . ~ . + s(day_number, by = included) + s(days_filled_out, by = included))
gamm4::gamm4(new_form, random = ~ (1 | person), data = reactive, family = family(obj)) -> control_reactivity
control_reactivity_mer = control_reactivity$mer
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
# findbars(formula(obj))[[1]]
if (class(obj) != "glmerMod") {
nlme::lme(fixed = nobars(formula(obj)), random = ~ 1 | person, data = diary, na.action = na.exclude, correlation = nlme::corAR1(form = ~ day_number | person), method = "REML") -> control_autocorrelation
control_autocorrelation$call$fixed = control_autocorrelation$terms
}
}, error = function(e) { cat_message(e, "danger") })
tryCatch({
if (class(obj) != "glmerMod") {
nlme::lme(fixed = nobars(formula(obj)), random = ~ 1 | person, data = diary, na.action = na.exclude, correlation = nlme::corARMA(form = ~ day_number | person, p = 1, q = 1), method = "REML") -> autocorrelation_moving_avg
autocorrelation_moving_avg$call$fixed = autocorrelation_moving_avg$terms
}
}, error = function(e) { cat_message(e, "danger") })
models_in_function = ls_type(c('lmerMod','glmerMod','bglmerMod','blmerMod','merModLmerTest','lme','lm'))
coefs = rbindlist(lapply(models_in_function, FUN = function(x) { get_coefs(get(x), x) }), fill = TRUE)
# eff_coefs = rbindlist(lapply(models_in_function, FUN = function(x) { get_eff_coefs(get(x), x) }), fill = TRUE)
coefs$model = dplyr::recode_factor(factor(coefs$model),
'obj' = 'M_1. Main model (all), BC+BCi',
'included_lax' = 'M_e2. Lax exclusion criteria',
'included_conservative' = 'M_e3. Conservative exclusion criteria',
'included_strict' = 'M_e4. Strict exclusion criteria',
'nododgy' = 'M_e5. Exclude potentially dodgy data',
'notrypreg' = 'M_e6. Exclude those trying to get pregnant',
'backward' = 'M_p1. Continuous, BC',
'forward_counting' = 'M_p2. Continuous, FC',
'broad_window' = 'M_p3. Broad window, BC',
'narrow_window' = 'M_p4. Narrow window, BC',
'forward_counting_narrow' = 'M_p5. Narrow window, FC',
'forward_counting_broad' = 'M_p6. Broad window, FC',
'backward_inferred' = 'M_p7. BC from rep. cycle length, when onset unknown',
'cycles_between_20_40' = 'M_p8. Only cycles 20-40 days long',
'between_person' = 'M_d1. Between, FC (first obs. only)',
'within_person_two_occasions' = 'M_d2. Within, FC, high/low 2 obs.',
'within_person_four_occasions' = 'M_d3. Within, FC, two high/low 4 obs.',
'no_cycles_shorter_than_20' = 'M_d4. No cycles shorter than 20 days',
'within_person_avg_hilow' = 'M_d5. Within, BC, average high/low all',
'control_self_esteem' = 'M_c1. Adjust for self esteem',
'dontcontrolavg' = 'M_c2. No adjustment for avg. fertility',
'fewer_controls' = 'M_c3. No adjustment for menstruation & avg. fertility',
'control_week' = 'M_c4. Control week day and number',
'control_time_of_response' = 'M_c5. Adjust for time of/for response',
'control_autocorrelation' = 'M_c6. Model autocorrelation, lag 1',
'autocorrelation_moving_avg' = 'M_c7. Model autocorrelation, moving avg., lag 1',
'broad_window_ctrl' = 'M_c8. Broad BC, adj. menstruation',
'control_reactivity_mer' = 'M_c9. Adjust for reactivity', .ordered = T)
coefs$model = factor(coefs$model, levels = rev(levels(coefs$model)))
coefs = coefs %>% mutate(term = recode(term, "fertile:includedhorm_contra" = "includedhorm_contra:fertile", "Xfertile" = "fertile", "Xincludedhorm_contra:fertile" = "includedhorm_contra:fertile"))
# eff_coefs$model = factor( car::Recode(eff_coefs$model, modeltranslate ))
coefs_all = coefs
coefs = coefs %>% filter(term == "fertile" | term == "includedhorm_contra:fertile")
cat("\n\n\n#### M_e: Exclusion criteria \n\n")
plote = ggplot(coefs %>% filter(model %begins_with% "M_e" | model %begins_with% "M_1. "), aes(x = model, y = estimate, ymax = conf.high, ymin = conf.low, colour = term), group = 1) +
geom_hline(yintercept = 0, linetype = "dotted", color = "gray70") +
geom_text(aes(label = round(estimate,2), y = estimate), position = position_dodge(width = 0.6), vjust = -0.7) +
geom_pointrange( position = position_dodge(width = 0.6), size = 1) +
scale_color_manual("Contraception status", values = c("includedhorm_contra:fertile"="black","fertile" = "red"), labels = c("includedhorm_contra:fertile"="hormonally\ncontracepting","fertile" = "fertile"), guide = F) +
xlab("Model") +
ylab(paste("Regression slope + 95 CI%")) +
ggtitle(outcome_label) +
coord_flip()
print(plote)
# cat("\n\n\n")
# coefs_all %>%
# filter(model %begins_with% "M_e" | model %begins_with% "M_1. ") %>%
# select(model, term, estimate, conf.low, conf.high) %>%
# pander() %>%
# cat()
cat("\n\n\n#### M_p: Predictors \n\n")
plotp = ggplot( coefs %>% filter(model %begins_with% "M_p" | model %begins_with% "M_1. "), aes(x = model, y = estimate, ymax = conf.high, ymin = conf.low, colour = term), group = 1) +
geom_hline(yintercept = 0, linetype = "dotted", color = "gray70") +
geom_text(aes(label = round(estimate,2), y = estimate), position = position_dodge(width = 0.6), vjust = -0.7) +
geom_pointrange( position = position_dodge(width = 0.6), size = 1) +
scale_color_manual("Contraception status", values = c("includedhorm_contra:fertile"="black","fertile"= "red"), labels = c("includedhorm_contra:fertile"="hormonally\ncontracepting","fertile"="fertile"), guide = F) +
xlab("Model") +
ylab(paste("Regression slope + 95 CI%")) +
ggtitle(outcome_label) +
coord_flip()
print(plotp)
#
cat("\n\n\n#### M_c: Covariates, controls, autocorrelation \n\n")
plotc = ggplot(coefs %>% filter(model %begins_with% "M_c" | model %begins_with% "M_1. "), aes(x = model, y = estimate, ymax = conf.high, ymin = conf.low, colour = term), group = 1) +
geom_hline(yintercept = 0, linetype = "dotted", color = "gray70") +
geom_text(aes(label = round(estimate,2), y = estimate), position = position_dodge(width = 0.6), vjust = -0.7) +
geom_pointrange( position = position_dodge(width = 0.6), size = 1) +
scale_color_manual("Contraception status", values = c("includedhorm_contra:fertile"="black","fertile"= "red"), labels = c("includedhorm_contra:fertile"="hormonally\ncontracepting","fertile"="fertile"), guide = F) +
xlab("Model") +
ylab(paste("Regression slope + 95 CI%")) +
ggtitle(outcome_label) +
coord_flip()
print(plotc)
tryCatch({
print_summary(control_reactivity$mer)
print_summary(control_reactivity$gam)
plot(control_reactivity$gam, pages = 1)
if (class(obj) != "glmerMod") {
print_summary(control_autocorrelation)
print_summary(autocorrelation_moving_avg)
} else {
cat_message("No AR1/ARMA autocorrelation models were fitted for binomial outcomes.", "info")
}
}, error = function(e) { cat_message(e, "danger") })
cat("\n\n\n#### M_d: Other designs \n\n")
plotd = ggplot( coefs %>% filter(model %begins_with% "M_d" | model %begins_with% "M_1. "), aes(x = model, y = estimate, ymax = conf.high, ymin = conf.low, colour = term), group = 1) +
geom_hline(yintercept = 0, linetype = "dotted", color = "gray70") +
geom_text(aes(label = round(estimate,2), y = estimate), position = position_dodge(width = 0.6), vjust = -0.7) +
geom_pointrange( position = position_dodge(width = 0.6), size = 1) +
scale_color_manual("Contraception status", values = c("includedhorm_contra:fertile"="black","fertile"= "red"), labels = c("includedhorm_contra:fertile"="hormonally\ncontracepting","fertile"="fertile"), guide = F) +
xlab("Model") +
ylab(paste("Estimate of effect on", outcome)) +
coord_flip()
print(plotd)
cat("\n\n\n#### _M_m1_: Moderation by contraceptive method \n\n
Based on the sample with lax exclusion criteria. Users who used any hormonal contraception are classified as hormonal, users who use any awareness-based methods (counting, temperature-based) are classified as 'fertility-awareness', women who don't fall into the before groups and use condoms, pessars, coitus interruptus etc. are classified as 'barrie or abstinence'. Women who don't use contraception or use other methods such as sterilisation are excluded from this analysis.
\n\n")
tryCatch({
update(obj, formula = . ~ . - included * (menstruation + fertile) + contraceptive_methods + (fertile + menstruation) + included:(fertile + menstruation),
data = diary, subset = !is.na(included_lax) & contraceptive_method != "other") -> add_main
update(obj, formula = . ~ . - included * (menstruation + fertile) + contraceptive_methods * (fertile + menstruation),
data = diary, subset = !is.na(included_lax) & contraceptive_method != "other") -> by_method
method_eff_coefs = sjp.int(by_method, type = "eff", showCI = F, printPlot = F)$data.list[[1]]
rec = c("hormonal" = "hormonally\ncontracepting","barrier_or_abstinence" = "only barrier\n(condoms, ...)\nor abstinence", "fertility_awareness" = "potentially\nfertility aware", "none" = "not using contraception")
method_eff_coefs$method = rec[as.character(method_eff_coefs$grp)]
eff_plot = ggplot(method_eff_coefs,
aes(x = x, y = y, ymax = conf.high, ymin = conf.low)) +
geom_smooth( size = 1, stat = "identity", color = 'black') +
xlab("Conception probability") +
facet_wrap(~ method) +
ylab(outcome)
print(eff_plot)
coefs = get_coefs(by_method, "by method") %>% filter(term != "(Intercept)")
plot = ggplot(coefs, aes(x = term, y = estimate, ymax = conf.high, ymin = conf.low, group = term)) +
geom_hline(yintercept = 0, linetype = "dotted", color = "gray70") +
geom_text(aes(label = round(estimate,2), y = estimate), position = position_dodge(width = 0.6), vjust = -0.7) +
geom_pointrange( position = position_dodge(width = 0.6), size = 1) +
xlab("Model") +
ylab(paste("Estimate of effect on", outcome)) +
coord_flip()
print(plot)
print_summary(by_method)
cat(pander(anova(add_main, by_method)))
}, error = function(e) { cat_message(e, "danger") })
invisible(obj)
}
#' ## Test a moderator
test_moderator = function (obj, moderator, diary, multiline = F, xlevels = 3) {
tryCatch({
add_main = update.formula(formula(obj), new = as.formula(paste0(". ~ . + ", moderator, " * included "))) # reorder so that the triptych looks nice
add_mod_formula = update.formula(update.formula(formula(obj), new = . ~ . - included * fertile), new = as.formula(paste0(". ~ . + ", moderator, " * included * fertile"))) # reorder so that the triptych looks nice
update(obj, formula = add_main) -> with_main
update(obj, formula = add_mod_formula) -> with_mod
cat(pander(anova(with_main, with_mod)))
plot_triptych(with_mod, x.var = "fertile", multiline = multiline, xlevels = xlevels)
print_summary(with_mod)
}, error = function(e) { cat_message(e, "danger") })
invisible(obj)
}