Risk Factors Associated with Osteoporosis: A Statistical Analysis of Demographic, Lifestyle, and Medical Predictors

# 1.1 Firstly we load the appropriate libraries
library(ggplot2)
library(dplyr)

library(patchwork)
library(pROC)

library(janitor)

# 1.2 Load data
setwd("/Users/warshaautar/Desktop/Statistics")
getwd()

data <- read.csv("osteoporosis.csv")

#1.3 Descriptive labels
str(data)

#1.2.1 Gender Count
tabyl(data$Gender)

#1.2.2 To begin, what are the mean and median ages in the data
summary(data$Age)

#1.2.3 Shapiro Sample test
shapiro.test(sample(data$Age, 500))

#1.2.4 Wilcoxon test
wilcox.test(Age ~ Osteoporosis, data = data)

#1.2.5 Plotting Figure 1
ggplot(data, aes(x = Age)) +
  geom_histogram(bins = 25, fill = "salmon", color = "white") +
  theme_minimal() +
  labs(title = "Figure 1: Visual Inspection of Age Distribution")

#1.2.6 Age to age category conversion
data <- data %>%
  mutate(
    Age_Group = cut(
      Age,
      breaks = c(0, 24, 39, 49, 59, 69, 79, 100),
      labels = c("<25", "26-40", "40-49", "50-59", "60-69", "70-79", "80+")
    )
  )

# 2.1 Listing all the variables
variables <- c("Age_Group","Smoking", "Hormonal.Changes", "Family.History", "Body.Weight",
               "Calcium.Intake", "Vitamin.D.Intake", "Physical.Activity",
               "Alcohol.Consumption", "Medical.Conditions", "Medications", "Prior.Fractures")

# 2.2 Running the Chi-Square test for each one automatically
results_list <- lapply(variables, function(var) {
  test <- chisq.test(data[[var]], data$Osteoporosis)

  data.frame(
    Factor = var,
    X_Squared = round(test$statistic, 3),
    df = test$parameter,
    P_Value = round(test$p.value, 4),
    Significant = ifelse(test$p.value < 0.05, "Yes (p < 0.05)", "No")
  )
})

# 2.3 Combine them all into one clean table
chi_square_table <- do.call(rbind, results_list)

chi_square_table

#Fig.2: Osteoporosis by Gender
demo1 <- ggplot(data, aes(x = Gender, fill = factor(Osteoporosis))) +
  geom_bar(position = "fill", alpha = 0.9) +
  scale_fill_discrete(labels = c("0" = "No", "1" = "Yes")) +
  theme_minimal() +
  labs(title = "Fig.2: Osteoporosis by Gender", x = "Gender", y = "Proportion", fill = "Osteoporosis")
scale_fill_discrete(labels = c("0" = "No", "1" = "Yes"))

#Fig.3:Osteoporosis by Age Group
demo2 <- ggplot(data, aes(x = Age_Group, fill = factor(Osteoporosis))) +
  geom_bar(position = "fill", alpha = 0.9) +
  scale_fill_discrete(labels = c("0" = "No", "1" = "Yes")) +
  theme_minimal() +
  labs(title = "Fig.3:Osteoporosis by Age Group", x = "Age Group", y = "Proportion", fill = "Osteoporosis")
scale_fill_discrete(labels = c("0" = "No", "1" = "Yes"))

#Fig.4:Osteoporosis by Race/Ethnicity
demo3 <- ggplot(data, aes(x = `Race.Ethnicity`, fill = factor(Osteoporosis))) +
  geom_bar(position = "fill", alpha = 0.9) +
  scale_fill_discrete(labels = c("0" = "No", "1" = "Yes")) +
  theme_minimal() +
  labs(title = "Fig.4:Osteoporosis by Race/Ethnicity", x = "Race / Ethnicity", y = "Proportion", fill = "Osteoporosis")


# Render Page 1
demo1 / demo2 / demo3 +
  plot_layout(guides = "collect") +
  plot_annotation(title = "Figures 2-4: Demographic Profiles & Osteoporosis Status")

#Fig.5: Physical Activity Level
life1 <- ggplot(data, aes(x = `Physical.Activity`, fill = factor(Osteoporosis))) +
  geom_bar(position = "fill", alpha = 0.9) +
  scale_fill_discrete(labels = c("0" = "No", "1" = "Yes")) +
  theme_minimal() +
  labs(title = "Fig.5: Physical Activity Level", x = "Activity Level", y = "Proportion", fill = "Osteoporosis")

#Fig.6: Vitamin D Intake
life2 <- ggplot(data, aes(x = `Vitamin.D.Intake`, fill = factor(Osteoporosis))) +
  geom_bar(position = "fill", alpha = 0.9) +
  scale_fill_discrete(labels = c("0" = "No", "1" = "Yes")) +
  theme_minimal() +
  labs(title = "Fig.6: Vitamin D Intake", x = "Vitamin D", y = "Proportion", fill = "Osteoporosis")

#Fig.7: Calcium Intake
life3 <- ggplot(data, aes(x = `Calcium.Intake`, fill = factor(Osteoporosis))) +
  geom_bar(position = "fill", alpha = 0.9) +
  scale_fill_discrete(labels = c("0" = "No", "1" = "Yes")) +
  theme_minimal() +
  labs(title = "Fig.7: Calcium Intake", x = "Calcium", y = "Proportion", fill = "Osteoporosis")

#Fig.8: Smoking Status
life4 <- ggplot(data , aes(x = `Smoking`, fill = factor(Osteoporosis))) +
  geom_bar(position = "fill", alpha = 0.9) +
  scale_fill_discrete(labels = c("0" = "No", "1" = "Yes")) +
  theme_minimal() +
  labs(title = "Fig.8: Smoking Status", x = "Smoking (Yes only)", y = "Proportion", fill = "Osteoporosis")

#Fig.9: Alcohol Consumption
life5 <- ggplot(data, aes(x = `Alcohol.Consumption`, fill = factor(Osteoporosis))) +
  geom_bar(position = "fill", alpha = 0.9) +
  scale_fill_discrete(labels = c("0" = "No", "1" = "Yes")) +
  theme_minimal() +
  labs(title = "Fig.9: Alcohol Consumption", x = "Alcohol Level", y = "Proportion", fill = "Osteoporosis")

#Render Page 2
(life1 + life2 + life3 + life4 + life5) +
  plot_layout(ncol = 2, guides = "collect") +
  plot_annotation(title = "Figures 5-9: Behavioral & Diet Risk Factors")

#Fig.10: Hormonal Status
med1 <- ggplot(data, aes(x = `Hormonal.Changes`, fill = factor(Osteoporosis))) +
  geom_bar(position = "fill", alpha = 0.9) +
  scale_fill_discrete(labels = c("0" = "No", "1" = "Yes")) +
  theme_minimal() +
  labs(title = "Fig.10: Hormonal Status", x = "Hormonal Changes", y = "Proportion", fill = "Osteoporosis")

#Fig.11: Family History
med2 <- ggplot(data, aes(x = `Family.History`, fill = factor(Osteoporosis))) +
  geom_bar(position = "fill", alpha = 0.9) +
  scale_fill_discrete(labels = c("0" = "No", "1" = "Yes")) +
  theme_minimal() +
  labs(title = "Fig.11: Family History", x = "Family History", y = "Proportion", fill = "Osteoporosis")

#Fig.12: Prior Fractures
med3 <- ggplot(data, aes(x = `Prior.Fractures`, fill = factor(Osteoporosis))) +
  geom_bar(position = "fill", alpha = 0.9) +
  scale_fill_discrete(labels = c("0" = "No", "1" = "Yes")) +
  theme_minimal() +
  labs(title = "Fig.12: Prior Fractures", x = "Prior Fractures", y = "Proportion", fill = "Osteoporosis")

#Fig.13: Medical Conditions
med4 <- ggplot(data, aes(x = `Medical.Conditions`, fill = factor(Osteoporosis))) +
  geom_bar(position = "fill", alpha = 0.9) +
  scale_fill_discrete(labels = c("0" = "No", "1" = "Yes")) +
  theme_minimal() +
  labs(title = "Fig.13: Medical Conditions", x = "Medical Conditions", y = "Proportion", fill = "Osteoporosis")

#Fig.14: Medications Use
med5 <- ggplot(data, aes(x = Medications, fill = factor(Osteoporosis))) +
  geom_bar(position = "fill", alpha = 0.9) +
  scale_fill_discrete(labels = c("0" = "No", "1" = "Yes")) +
  theme_minimal() +
  labs(title = "Fig.14: Medications Use", x = "Medications", y = "Proportion", fill = "Osteoporosis")

# Render Page 3
(med1 + med2 + med3 + med4 + med5) +
  plot_layout(ncol = 2, guides = "collect") +
  plot_annotation(title = "Figures 10-14: Medical History & Clinical Factors; Hormonal Changes, Family History, Prior Fractures")

#4.1.1 Splits data and generates a summary
data %>%
  group_by(factor(Osteoporosis)) %>%
  summarise(
    Average_Age = mean(Age),
    Median_Age = median(Age),
    Min_Age = min(Age),
    Max_Age = max(Age),
    Total_Patients = n()
  )

# 4.1.2 Draw a Boxplot to once again visually prove how different the ages are
# This graph will clearly show that the groups barely overlap!
ggplot(data, aes(x = factor(Osteoporosis), y = Age, fill = factor(Osteoporosis))) +
  geom_boxplot(alpha = 0.7, outlier.color = "red") +
  scale_fill_discrete(labels = c("0" = "No", "1" = "Yes")) +
  theme_minimal() +
  labs(
    title = "The Primary Driver: Age Distribution by Osteoporosis Status",
    subtitle = "Proving that age profiles do not overlap significantly between groups",
    x = "Has Osteoporosis? (No vs Yes)",
    y = "Age (Years)",
    fill = "Osteoporosis"
  )

# 4.1.3 Fit a model using ONLY Age as a predictor
age_only_model <- glm(Osteoporosis ~ Age,
                       data = data,
                       family = binomial)

# 4.1.4 Look at the statistical significance of Age by itself
summary(age_only_model)

# 4.1.5 Generate risk probabilities based solely on a patient's age
data$predicted_risk_age <- predict(age_only_model, type = "response")

# 4.1.6 Calculate the ROC curve data
roc_age <- roc(data$Osteoporosis, data$predicted_risk_age)

# 4.1.7 Plot the ROC curve for Age alone
plot(roc_age,
     col = "darkgreen",   # Using a different color to distinguish it
     lwd = 3,
     main = paste("Age Standalone Accuracy (AUC =", round(auc(roc_age), 2), ")"))

# 4.2 XXXX NEED TO LOOK AT THIS Summary Statistics for Categorical Variables by Osteoporosis Status
# Since everything is categorical, let's look at counts for key risk factors
summary_stats <- data %>%
  group_by(Osteoporosis = factor(Osteoporosis)) %>%
  summarise(
    Total_Patients = n(),
    # Counts for Calcium Intake levels
    Calcium_Low = sum(Calcium.Intake == "Low", na.rm = TRUE),
    Calcium_Adequate = sum(Calcium.Intake == "Adequate", na.rm = TRUE),
    # Counts for Vitamin D Intake levels
    VitD_Insufficient = sum(Vitamin.D.Intake == "Insufficient", na.rm = TRUE),
    VitD_Sufficient = sum(Vitamin.D.Intake == "Sufficient", na.rm = TRUE),
    # Count of Smokers
    Smokers = sum(Smoking == "Yes", na.rm = TRUE)
  )
print(summary_stats)

# 4.2.XXX Fit the Multivariable Logistic Regression Model (EXCLUDING Age)
# All column names have been updated to match R's dot-notation conversion of your dataset
model_no_age <- glm(
  factor(Osteoporosis) ~ Gender + Hormonal.Changes + Family.History +
                         Race.Ethnicity + Body.Weight + Calcium.Intake +
                         Vitamin.D.Intake + Physical.Activity + Smoking +
                         Alcohol.Consumption + Medical.Conditions +
                         Medications + Prior.Fractures,
  data = data,
  family = binomial
)

# 4.2.3 Look at the statistical significance of the remaining variables
summary(model_no_age)

# 4.4.4 Generate risk probabilities based on this new model
data$predicted_risk_no_age <- predict(model_no_age, type = "response")

# 4.4.5 Calculate the ROC curve data
# Making sure response is matching the format used in glm
roc_no_age <- roc(data$Osteoporosis, data$predicted_risk_no_age)

# 4.4.6 Plot the ROC Curve for the model without Age
plot(
  roc_no_age,
  col = "darkblue",
  lwd = 3,
  main = paste("Model Accuracy Without Age (AUC =", round(auc(roc_no_age), 2), ")")
)

# 5.1.1 Baseline Model: Age Only
model_baseline <- glm(
  factor(Osteoporosis) ~ Age,
  data = data,
  family = binomial
)

# 5.1.2 Expanded Model: Age + All Clinical & Lifestyle Factors
model_expanded <- glm(
  factor(Osteoporosis) ~ Age + Gender + Hormonal.Changes + Family.History +
                         Race.Ethnicity + Body.Weight + Calcium.Intake +
                         Vitamin.D.Intake + Physical.Activity + Smoking +
                         Alcohol.Consumption + Medical.Conditions +
                         Medications + Prior.Fractures,
  data = data,
  family = binomial
)


# 5.2.2 COMPUTE INFORMATION FIT DIAGNOSTICS (AIC)
aic_baseline <- AIC(model_baseline)
aic_expanded <- AIC(model_expanded)

cat("--- Information Fit Diagnostics ---\n")

cat("Baseline Model (Age Only) AIC: ", round(aic_baseline, 2), "\n")

cat("Expanded Model (Age + All) AIC:", round(aic_expanded, 2), "\n")

cat("AIC Difference:                ", round(aic_baseline - aic_expanded, 2), "\n\n")

# Note: A lower AIC is better. If the expanded model's AIC is NOT significantly
# lower than the baseline, the extra variables are just adding noise.


#5.3.1 Generate predictions
prob_baseline <- predict(model_baseline, type = "response")
prob_expanded <- predict(model_expanded, type = "response")

#5.3.2 Calculate ROC
roc_baseline <- roc(data$Osteoporosis, prob_baseline)

roc_expanded <- roc(data$Osteoporosis, prob_expanded)

#5.3.3 Run DeLong's test to see if the two AUCs are statistically different
delong_test <- roc.test(roc_baseline, roc_expanded, method = "delong")
print(delong_test)

#5.4.1 Plot the baseline curve first
plot(
  roc_baseline,
  col = "darkgreen",
  lwd = 3,
  main = "Comparative Diagnostic Evaluation",
  xlab = "False Positive Rate (1 - Specificity)",
  ylab = "True Positive Rate (Sensitivity)"
)

#5.4.2 Add the expanded model curve to the same plot
plot(
  roc_expanded,
  col = "darkblue",
  lwd = 3,
  add = TRUE
)

#5.4.3 Add a legend to distinguish the models, including AIC and AUC metrics
legend(
  "bottomright",
  legend = c(
    paste("Age Only (AUC =", round(auc(roc_baseline), 2), ", AIC =", round(aic_baseline, 0), ")"),
    paste("Age + Clinical (AUC =", round(auc(roc_expanded), 2), ", AIC =", round(aic_expanded, 0), ")")
  ),
  col = c("darkgreen", "darkblue"),
  lwd = 3,
  bty = "n" # Removes the legend border box
)