Marginal Effects in Controlled Models

Data

source("0_helpers.R")
library(modelr)
library(tidybayes)

load("data/cleaned_selected_wrangled.rdata")

Sexual Frequency (Penetrative Intercourse)

Uncontrolled Model

m_hc_sexfreqpen = brm(diary_sex_active_sex_sum ~
                        offset(log(number_of_days)) +
                        contraception_hormonal,
                data = data, family = poisson(),
                file = "m_hc_sexfreqpen")

Data Grid

grid = data %>%
  data_grid(number_of_days = 1,
            contraception_hormonal)

Add Predicted Values

data_pred = add_epred_draws(grid, m_hc_sexfreqpen,
                   value = "E[y|A,B]",
                   ndraws = 1000)

Average Marginal Effect

ame = data_pred %>%
  group_by(contraception_hormonal, .draw) %>%  # group by predictors to keep
  summarise(`E[y|A]` = mean(`E[y|A,B]`, na.rm = T)) %>%
  compare_levels(`E[y|A]`, by = contraception_hormonal) %>%  
  rename(`mean difference` = `E[y|A]`) %>%
  mean_qi(`mean difference`, .width = .90)

## `summarise()` has grouped output by 'contraception_hormonal'. You can override using the `.groups` argument.

ame

## # A tibble: 1 x 7
##   contraception_hormonal `mean difference` .lower .upper .width .point .interval
##   <chr>                              <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1 yes - no                          0.0350 0.0288 0.0410    0.9 mean   qi

Compare Means

means = data_pred %>%
  group_by(contraception_hormonal) %>%
  summarise(mean_fit = mean(`E[y|A,B]`, na.rm = T))
print(data.frame(means), digits = 4)

##   contraception_hormonal mean_fit
## 1                     no    0.124
## 2                    yes    0.159

Controlled Model

m_hc_sexfreqpen_controlled = brm(diary_sex_active_sex_sum ~
                        offset(log(number_of_days)) +
                        contraception_hormonal +
                        age + net_income + relationship_duration_factor +
                              education_years +
                              bfi_extra + bfi_neuro + bfi_agree + bfi_consc + bfi_open +
                              religiosity,
                data = data, family = poisson(),
                file = "m_hc_sexfreqpen_controlled")

Data Grid

grid = data %>%
  data_grid(number_of_days = 1,
            contraception_hormonal,
            age = seq_range(age, n = 5), # measured in full years
            net_income,
            relationship_duration_factor,
            education_years = seq_range(education_years, n = 5), # measured in full years
            bfi_extra = mean(data$bfi_extra),
            bfi_neuro = mean(data$bfi_neuro), 
            bfi_agree = mean(data$bfi_agree), 
            bfi_consc = mean(data$bfi_consc), 
            bfi_open = mean(data$bfi_open),
            religiosity = mean(data$religiosity)
  )

Add Predicted Values

data_pred = add_epred_draws(grid, m_hc_sexfreqpen_controlled,
                   value = "E[y|A,B]",
                   ndraws = 1000)

Average Marginal Effect

ame = data_pred %>%
  group_by(contraception_hormonal, .draw) %>%  # group by predictors to keep
  summarise(`E[y|A]` = mean(`E[y|A,B]`, na.rm = T)) %>%
  compare_levels(`E[y|A]`, by = contraception_hormonal) %>%  
  rename(`mean difference` = `E[y|A]`) %>%
  mean_qi(`mean difference`, .width = .90)

## `summarise()` has grouped output by 'contraception_hormonal'. You can override using the `.groups` argument.

ame

## # A tibble: 1 x 7
##   contraception_hormonal `mean difference` .lower .upper .width .point .interval
##   <chr>                              <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1 yes - no                          0.0272 0.0194 0.0348    0.9 mean   qi

Compare Means

means = data_pred %>%
  group_by(contraception_hormonal) %>%
  summarise(mean_fit = mean(`E[y|A,B]`, na.rm = T))
print(as.data.frame(means), digits = 4)

##   contraception_hormonal mean_fit
## 1                     no   0.1525
## 2                    yes   0.1797

Masturbation Frequency

Uncontrolled Model

m_hc_masfreq = brm(diary_masturbation_sum ~ offset(log(number_of_days)) +
                     contraception_hormonal,
                data = data, family = poisson(),
                file = "m_hc_masfreq")

Data Grid

grid = data %>%
  data_grid(number_of_days = 1,
            contraception_hormonal)

Add Predicted Values

data_pred = add_epred_draws(grid, m_hc_masfreq,
                   value = "E[y|A,B]",
                   ndraws = 1000)

Average Marginal Effect

ame = data_pred %>%
  group_by(contraception_hormonal, .draw) %>%  # group by predictors to keep
  summarise(`E[y|A]` = mean(`E[y|A,B]`, na.rm = T)) %>%
  compare_levels(`E[y|A]`, by = contraception_hormonal) %>%  
  rename(`mean difference` = `E[y|A]`) %>%
  mean_qi(`mean difference`, .width = .90)

## `summarise()` has grouped output by 'contraception_hormonal'. You can override using the `.groups` argument.

ame

## # A tibble: 1 x 7
##   contraception_hormonal `mean difference`  .lower  .upper .width .point .interval
##   <chr>                              <dbl>   <dbl>   <dbl>  <dbl> <chr>  <chr>    
## 1 yes - no                         -0.0502 -0.0556 -0.0450    0.9 mean   qi

Compare Means

means = data_pred %>%
  group_by(contraception_hormonal) %>%
  summarise(mean_fit = mean(`E[y|A,B]`, na.rm = T))
print(as.data.frame(means), digits = 4)

##   contraception_hormonal mean_fit
## 1                     no   0.1531
## 2                    yes   0.1029

Controlled Model

m_hc_masfreq_controlled = brm(diary_masturbation_sum ~
                        offset(log(number_of_days)) +
                        contraception_hormonal +
                        age + net_income + relationship_duration_factor +
                              education_years +
                              bfi_extra + bfi_neuro + bfi_agree + bfi_consc + bfi_open +
                              religiosity,
                data = data, family = poisson(),
                file = "m_hc_masfreq_controlled")

Data Grid

grid = data %>%
  data_grid(number_of_days = 1,
            contraception_hormonal,
            age = seq_range(age, n = 5), # measured in full years
            net_income,
            relationship_duration_factor,
            education_years = seq_range(education_years, n = 5), # measured in full years
            bfi_extra = mean(data$bfi_extra),
            bfi_neuro = mean(data$bfi_neuro), 
            bfi_agree = mean(data$bfi_agree), 
            bfi_consc = mean(data$bfi_consc), 
            bfi_open = mean(data$bfi_open),
            religiosity = mean(data$religiosity)
  )

Add Predicted Values

data_pred = add_epred_draws(grid, m_hc_masfreq_controlled,
                   value = "E[y|A,B]",
                   ndraws = 1000)

Average Marginal Effect

ame = data_pred %>%
  group_by(contraception_hormonal, .draw) %>%  # group by predictors to keep
  summarise(`E[y|A]` = mean(`E[y|A,B]`, na.rm = T)) %>%
  compare_levels(`E[y|A]`, by = contraception_hormonal) %>%  
  rename(`mean difference` = `E[y|A]`) %>%
  mean_qi(`mean difference`, .width = .90)

## `summarise()` has grouped output by 'contraception_hormonal'. You can override using the `.groups` argument.

ame

## # A tibble: 1 x 7
##   contraception_hormonal `mean difference`  .lower  .upper .width .point .interval
##   <chr>                              <dbl>   <dbl>   <dbl>  <dbl> <chr>  <chr>    
## 1 yes - no                         -0.0311 -0.0358 -0.0266    0.9 mean   qi

Compare Means

means = data_pred %>%
  group_by(contraception_hormonal) %>%
  summarise(mean_fit = mean(`E[y|A,B]`, na.rm = T))
print(as.data.frame(means), digits = 4)

##   contraception_hormonal mean_fit
## 1                     no  0.12145
## 2                    yes  0.09033