Hormonal Birth Control and Estrogen Levels Analysis

Introduction

This report investigates whether women on hormonal birth control have higher basal levels of estrogen than those who are not.

Setup & Data Import

Load required libraries

rm(list=ls())

library(tidyverse)
library(readxl)
library(car)
library(FSA)

Load BC & Aurora data

bc_data <- read_excel("AURORA_Contraceptives Progression.xlsx", sheet = "Pre_Contcp") %>%
  select(PID, `Progestin only / Estrogen containing`) %>%
  rename(BC = `Progestin only / Estrogen containing`)

aurora_data <- read_csv("full_aurora.csv")

Load Estrogen data

# Load estrogen concentration data
e2_conc <- read.csv("E2_concentrated.csv") %>%
  mutate(
    mean_conc = as.numeric(Mean.Concentration..pg.mL.),
    CV = as.numeric(`X..CV`)
  ) %>%
  filter(
    mean_conc > 0,
    mean_conc < 420,
    CV > 0
  ) %>%
  select(Inventory.Code, mean_conc, CV) %>%
  mutate(sqrt_conc = sqrt(mean_conc))

# Load PID key for estrogen data
pid_key <- read_excel("e2_pid_key.xlsx")

Data Cleaning & Prep

Key E2 Data with PID

e2_keyed <- e2_conc %>%
  inner_join(pid_key, by = c("Inventory.Code" = "InventoryCode"))

# Filter for specific time points
e2_time_filtered <- e2_keyed %>%
  filter(AlternateID == "T0")  # Basal level at time of ED visit

Merge E2 + Aurora

e2_full <- e2_time_filtered %>%
  left_join(aurora_data, by = "PID") %>%
  select(
    PID,
    ED_Age,
    ED_HighestGrade,
    BMI,
    ED_RaceEthCode,
    mean_conc,
    ED_NowPain_C,
    WK8_Pain_C,
    M3_Pain_C,
    M6_Pain_C,
    AlternateID
  ) %>%
  rename(
    ED_Pain = ED_NowPain_C,
    WK8_Pain = WK8_Pain_C,
    M3_Pain = M3_Pain_C,
    M6_Pain = M6_Pain_C,
    Race = ED_RaceEthCode
  ) %>%
  filter(!is.na(mean_conc))

Merge BC + E2

# Merge birth control data with estrogen and Aurora data
combined_data <- bc_data %>%
  full_join(e2_full, by = "PID") %>%
  filter(AlternateID == "T0")  # Ensure basal estrogen levels

# Recode NA values in the 'BC' column to 'None'
combined_data <- combined_data %>%
  mutate(BC = ifelse(is.na(BC), "None", BC))

# Convert 'BC' to a factor with explicit levels
combined_data <- combined_data %>%
  mutate(BC = factor(BC, levels = c("Progestin only", "Estrogen containing", "None")))

Analysis

Summary Statistics

# Calculate mean mean_conc for each BC category
mean_estrogen <- combined_data %>%
  group_by(BC) %>%
  summarise(mean_sqrt_conc = mean(mean_conc, na.rm = TRUE),
            count = n())

# Display the summary
print(mean_estrogen)

Visualizations

gghisto <- list(
  theme(
    axis.text.x = element_text(face = "bold", color = "cornflowerblue", size = 12),
    axis.text.y = element_text(face = "bold", color = "royalblue4", size = 12),
    axis.title = element_text(size = 14, face = "bold"),
    plot.title = element_text(size = 14, face = "bold.italic"),
    panel.border = element_rect(color = "black", fill = NA, size = 1)))

ggplot(combined_data, aes(x = BC, y = mean_conc, color = BC)) +
  geom_boxplot(color = "black", alpha = 0.8, width = 0.6) +
  geom_jitter(aes(fill = BC),shape = 21, color = "black",  
    width = 0.2, alpha = 0.7, size = 4      
  ) +
  labs(
    title = "Basal Estrogen Levels by Birth Control Category",
    x = "Birth Control Category",
    y = "E2 Concentration (pg/mL)"
  ) +
  gghisto +
  theme(legend.position = "none")

Normality/Homogeneity

analysis_data <- combined_data

# Shapiro-Wilk test for each group
shapiro_progestin <- shapiro.test(analysis_data$mean_conc[analysis_data$BC == "Progestin only"])
shapiro_estrogen <- shapiro.test(analysis_data$mean_conc[analysis_data$BC == "Estrogen containing"])

# Print results
print(shapiro_progestin)

    Shapiro-Wilk normality test

data:  analysis_data$mean_conc[analysis_data$BC == "Progestin only"]
W = 0.66494, p-value = 2.546e-06

print(shapiro_estrogen)

    Shapiro-Wilk normality test

data:  analysis_data$mean_conc[analysis_data$BC == "Estrogen containing"]
W = 0.41389, p-value = 4.381e-06

# Q-Q Plots
ggplot(analysis_data, aes(sample = mean_conc, color = BC)) +
  stat_qq() +
  stat_qq_line() +
  facet_wrap(~ BC) +
  theme_minimal() +
  labs(title = "Q-Q Plots for Estrogen Concentration by BC Category") + gghisto+
  theme(legend.position = "none")

# Levene's Test
levene_test <- leveneTest(mean_conc ~ BC, data = analysis_data)
print(levene_test)

Levene's Test for Homogeneity of Variance (center = median)
       Df F value Pr(>F)
group   2  1.0734 0.3425
      555               

KW Test & Dunn’s

kruskal_result <- kruskal.test(mean_conc ~ BC, data = analysis_data)
dunn_test <- dunnTest(mean_conc ~ BC, data = analysis_data)

print(kruskal_result);print(dunn_test)

    Kruskal-Wallis rank sum test

data:  mean_conc by BC
Kruskal-Wallis chi-squared = 0.63846, df = 2, p-value = 0.7267

Dunn (1964) Kruskal-Wallis multiple comparison
  p-values adjusted with the Holm method.