# ==============================================================================
# BRFSS 2024 DATA CLEANING & DESCRIPTIVE REPORT
# Author: Herve Fossou / Samar Ikram
# Focus: Menopausal Status and Osteoarthritis
# ==============================================================================

# 1. Clear environment memory
rm(list = ls())

# 2. Load required packages
library(tidyverse)
library(haven)

# 3. Optimized Data Ingestion
brfss_subset <- read_xpt(
  "LLCP2024.XPT", 
  col_select = c("SEXVAR", "_LLCPWT", "_STATE", "_PSU", "HAVARTH4", "_AGEG5YR", "_BMI5")
) %>% 
  # Filter immediately for female respondents (SEXVAR == 2)
  filter(SEXVAR == 2)

# 4. Handling Missing Values and Exclusions
brfss_missing_handled <- brfss_subset %>%
  filter(
    HAVARTH4 %in% c(1, 2),                  # Keep only explicit 'Yes' or 'No'
    `_AGEG5YR` %in% 1:14,                   # Exclude 'Refused/Missing' age categories
    `_BMI5` > 0 & `_BMI5` != 99999          # Exclude missing BMI values
  )

# 5. Variable Recoding and Collapsing Categories
brfss_clean <- brfss_missing_handled %>%
  mutate(
    # Recode Osteoarthritis binary outcome
    osteoarthritis = case_when(
      HAVARTH4 == 1 ~ "Yes",
      HAVARTH4 == 2 ~ "No"
    ),
    osteoarthritis = factor(osteoarthritis, levels = c("No", "Yes")),
    
    # Collapse 14 raw age brackets into 3 biological lifecycle periods
    menopause_status = case_when(
      `_AGEG5YR` %in% 1:5   ~ "Pre-Menopausal Age (18-44)",
      `_AGEG5YR` %in% 6:7   ~ "Peri/Early Post-Menopausal Age (45-54)",
      `_AGEG5YR` %in% 8:14  ~ "Post-Menopausal Age (55+)"
    ),
    menopause_status = factor(menopause_status, 
                              levels = c("Pre-Menopausal Age (18-44)", 
                                         "Peri/Early Post-Menopausal Age (45-54)", 
                                         "Post-Menopausal Age (55+)")),
    
    # Correct the BRFSS 2-decimal implicit multiplier for BMI
    bmi = `_BMI5` / 100
  ) %>%
  # Retain only final structural variables for the analytic file
  select(`_LLCPWT`, `_STATE`, `_PSU`, osteoarthritis, menopause_status, bmi)

# 6. Applying Clear Variable Metadata Labels
attr(brfss_clean$osteoarthritis, "label")   <- "Doctor-Diagnosed Osteoarthritis Status"
attr(brfss_clean$menopause_status, "label") <- "Biological Menopausal Age Categories"
attr(brfss_clean$bmi, "label")              <- "Calculated Body Mass Index (kg/m²)"


# ==============================================================================
# SECTION 7: STUDY POPULATION CHARACTERISTICS (DEMOGRAPHICS & COVARIATES)
# ==============================================================================

# 1. Summary of Categorical Variables (Frequencies & Percentages)
cat("\n--- Osteoarthritis Frequencies ---\n")
## 
## --- Osteoarthritis Frequencies ---
oa_table <- table(brfss_clean$osteoarthritis)
print(oa_table)
## 
##     No    Yes 
## 125226  84300
print(prop.table(oa_table) * 100)
## 
##       No      Yes 
## 59.76633 40.23367
cat("\n--- Menopausal Age Categories Frequencies ---\n")
## 
## --- Menopausal Age Categories Frequencies ---
meno_table <- table(brfss_clean$menopause_status)
print(meno_table)
## 
##             Pre-Menopausal Age (18-44) Peri/Early Post-Menopausal Age (45-54) 
##                                  57543                                  28387 
##              Post-Menopausal Age (55+) 
##                                 123596
print(prop.table(meno_table) * 100)
## 
##             Pre-Menopausal Age (18-44) Peri/Early Post-Menopausal Age (45-54) 
##                               27.46342                               13.54820 
##              Post-Menopausal Age (55+) 
##                               58.98838
# 2. Summary of Continuous Variable (BMI Mean & Standard Deviation)
cat("\n--- Continuous BMI Summary ---\n")
## 
## --- Continuous BMI Summary ---
bmi_summary <- brfss_clean %>%
  summarise(
    Mean_BMI = mean(bmi, na.rm = TRUE),
    SD_BMI = sd(bmi, na.rm = TRUE),
    Median_BMI = median(bmi, na.rm = TRUE),
    IQR_BMI = IQR(bmi, na.rm = TRUE)
  )
print(bmi_summary)
## # A tibble: 1 × 4
##   Mean_BMI SD_BMI Median_BMI IQR_BMI
##      <dbl>  <dbl>      <dbl>   <dbl>
## 1     28.5   7.13       27.4    8.63
# 3. Visualizations: Bar Chart for Osteoarthritis by Menopausal Status
ggplot(brfss_clean, aes(x = menopause_status, fill = osteoarthritis)) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent) +
  theme_minimal() +
  labs(
    title = "Osteoarthritis Prevalence Across Menopausal Age Groups",
    x = "Biological Menopausal Age Categories",
    y = "Percentage (%)",
    fill = "Osteoarthritis"
  ) +
  theme(axis.text.x = element_text(angle = 15, hjust = 1))

# 4. Visualizations: Histogram of Body Mass Index (BMI)
ggplot(brfss_clean, aes(x = bmi)) +
  geom_histogram(binwidth = 1, fill = "steelblue", color = "white", alpha = 0.8) +
  geom_vline(aes(xintercept = mean(bmi, na.rm = TRUE)), color = "red", linetype = "dashed", size = 1) +
  xlim(15, 50) +  # Limit x-axis to standard ranges for clean viewing
  theme_minimal() +
  labs(
    title = "Distribution of Body Mass Index (BMI) in Study Population",
    subtitle = "Red dashed line indicates the population mean",
    x = "Calculated BMI (kg/m²)",
    y = "Count"
  )