# =========================================================
# R Script – Predictive tasks examples using the synthetic obstetric dataset (CSV format)
# =========================================================

# Install necessary packages
install.packages(c("tidyverse", "caret", "pROC", "broom", "dplyr", "gtsummary"))

library(tidyverse)
library(caret)
library(pROC)
library(broom)
library(dplyr)
library(gtsummary)

# =========================================================
# 1. Load CSV dataset
# =========================================================

# Update path as needed
setwd("C:/Users/...")
data <- read.csv("Pregnant patients - Synthetic Dataset.csv")

# Inspect dataset global structure
str(data)

# =========================================================
# 2. Variable selection
# =========================================================

# Manual selection of clinically relevant predictors

model_data <- data %>%
  select(
    maternal_age,
    bmi,
    fasting_glucose,
    blood_pressure_week16_systolic,
    blood_pressure_week16_diastolic,
    blood_pressure_week29_systolic,
    blood_pressure_week29_diastolic,
    hemoglobin_first_trimester,
    tobacco_use,
    alcohol_use,
    pregestational_diabetes,
    gestational_diabetes = referred_to_gestational_diabetes_clinic,
    preeclampsia = preeclampsia,
    gestational_hypertension,
    macrosomia,
    preterm_birth
  )

# =========================================================
# 3. Basic data pre-processing
# =========================================================

# Convert binary variables to factors (categorical)
binary_vars <- c(
  "tobacco_use",
  "alcohol_use",
  "pregestational_diabetes",
  "gestational_diabetes",
  "preeclampsia",
  "gestational_hypertension",
  "preterm_birth"
)

model_data[binary_vars] <- lapply(model_data[binary_vars], factor)

# Remove missing values for simplicity
model_data <- na.omit(model_data)

# Create maternal age group categories, covert to categorical and define a reference level (age 25 or lower)
model_data <- model_data %>%
  mutate(
    maternal_age_group = case_when(
      maternal_age < 25 ~ "<25",
      maternal_age >= 25 & maternal_age <= 29 ~ "25-29",
      maternal_age >= 30 & maternal_age <= 34 ~ "30-34",
      maternal_age >= 35 ~ "35+"
    )
  )

model_data$maternal_age_group <- factor(
  model_data$maternal_age_group,
  levels = c("<25", "25-29", "30-34", "35+")
)

# Create BMI categories, covert to categorical and define a reference level (age 25 or lower)
model_data <- model_data %>%
  mutate(
    bmi_group = case_when(
      bmi < 18.5 ~ "Underweight",
      bmi >= 18.5 & bmi < 25 ~ "Normal",
      bmi >= 25 & bmi < 30 ~ "Overweight",
      bmi >= 30 ~ "Obesity"
    )
  )

model_data$bmi_group <- factor(
  model_data$bmi_group,
  levels = c("Normal", "Underweight", "Overweight", "Obesity")
)

# =========================================================
# 4. Descriptive Analysis
# =========================================================

# Overall descriptive statistics table
table_desc <- model_data %>%
  tbl_summary()
table_desc

# Descriptive statistics table by outcome - Gestational Diabetes
table_desc_GDM <- model_data %>%
  tbl_summary(by = gestational_diabetes) %>% add_p()
table_desc_GDM

# Descriptive statistics table by outcome - Gestational Hypertension
table_desc_GHTN <- model_data %>%
  tbl_summary(by = gestational_hypertension) %>% add_p()
table_desc_GHTN

# Descriptive statistics table by outcome - Preterm birth
table_desc_PTB <- model_data %>%
  tbl_summary(by = preterm_birth) %>% add_p()
table_desc_PTB

# Descriptive statistics table by outcome - Macrosomia
table_desc_Macros <- model_data %>%
  tbl_summary(by = macros) %>% add_p()
table_desc_Macros

# Descriptive statistics table by outcome - Preeclampsia
table_desc_PE <- model_data %>%
  tbl_summary(by = preeclampsia) %>% add_p()
table_desc_PE


# =========================================================
# 5a. Logistic Regression – Gestational Diabetes
# =========================================================
set.seed(123) #Ensures reproducibility

# Train-Test Split
train_index <- createDataPartition(
  model_data$gestational_diabetes,
  p = 0.8,
  list = FALSE
)

train_data <- model_data[train_index, ]
test_data  <- model_data[-train_index, ]

# Model estimation (using Generalized Logistic Regression)
gdm_model <- glm(
  gestational_diabetes ~ maternal_age_group + bmi_group + fasting_glucose +
    alcohol_use + pregestational_diabetes,
  data = train_data,
  family = binomial()
)

#Overview of the models' coefficients, confidence intervals and p-values
summary(gdm_model)

# Predictions
gdm_probs <- predict(gdm_model, newdata = test_data, type = "response")

# ROC-AUC
gdm_roc <- roc(test_data$gestational_diabetes, gdm_probs)
print(gdm_roc)
print(auc(gdm_roc))

# Odds ratios
options(scipen = 999)
exp(cbind(
  OddsRatio = coef(gdm_model),
  confint(gdm_model)
))

# =========================================================
# 5b. Logistic Regression – Gestational Hypertension
# =========================================================
set.seed(123) #Ensures reproducibility

#Analyse considering 5-unit variation in BMI
model_data$bmi_5unit <- model_data$bmi / 5

# Train-Test Split
train_index <- createDataPartition(
  model_data$gestational_hypertension,
  p = 0.8,
  list = FALSE
)

train_data <- model_data[train_index, ]
test_data  <- model_data[-train_index, ]

# Model estimation (using Generalized Logistic Regression)
ghtn_model <- glm(
  gestational_hypertension ~ maternal_age_group + bmi_5unit + pregestational_diabetes,
  data = train_data,
  family = binomial()
)

#Overview of the models' coefficients, confidence intervals and p-values
summary(ghtn_model)

# Predictions
ghtn_probs <- predict(ghtn_model, newdata = test_data, type = "response")

# ROC-AUC
ghtn_roc <- roc(test_data$gestational_hypertension, ghtn_probs)
print(ghtn_roc)
print(auc(ghtn_roc))

# Odds ratios
options(scipen = 999)
exp(cbind(
  OddsRatio = coef(ghtn_model),
  confint(ghtn_model)
))

# =========================================================
# 5c. Logistic Regression – Preterm Births
# =========================================================
set.seed(123) #Ensures reproducibility

# Train-Test Split
train_index <- createDataPartition(
  model_data$preterm_birth,
  p = 0.8,
  list = FALSE
)

train_data <- model_data[train_index, ]
test_data  <- model_data[-train_index, ]

# Model estimation (using Generalized Logistic Regression)
ptb_model <- glm(
  preterm_birth ~ maternal_age_group + bmi_group + preeclampsia + pregestational_diabetes,
  data = train_data,
  family = binomial()
)

#Overview of the models' coefficients, confidence intervals and p-values
summary(ptb_model)

# Predictions
ptb_probs <- predict(ptb_model, newdata = test_data, type = "response")

# ROC-AUC
ptb_roc <- roc(test_data$preterm_birth, ptb_probs)
print(ptb_roc)
print(auc(ptb_roc))

# Odds ratios
options(scipen = 999)
exp(cbind(
  OddsRatio = coef(ptb_model),
  confint(ptb_model)
))


# =========================================================
# 5d. Logistic Regression – Macrosomia
# =========================================================
set.seed(123) #Ensures reproducibility

# Train-Test Split
train_index <- createDataPartition(
  model_data$macrosomia,
  p = 0.8,
  list = FALSE
)

train_data <- model_data[train_index, ]
test_data  <- model_data[-train_index, ]

# Model estimation (using Generalized Logistic Regression)
macros_model <- glm(
  macrosomia ~ maternal_age_group + bmi_group + fasting_glucose + gestational_diabetes,
  data = train_data,
  family = binomial()
)

#Overview of the models' coefficients, confidence intervals and p-values
summary(macros_model)

# Predictions
macros_probs <- predict(macros_model, newdata = test_data, type = "response")

# ROC-AUC
macros_roc <- roc(test_data$macrosomia, macros_probs)
print(macros_roc)
print(auc(macros_roc))

# Odds ratios
options(scipen = 999)
exp(cbind(
  OddsRatio = coef(macros_model),
  confint(macros_model)
))


# =========================================================
# 5e. Logistic Regression – Preeclampsia
# =========================================================
set.seed(123) #Ensures reproducibility

#Analyse considering 5-unit variation in BMI
model_data$bmi_5unit <- model_data$bmi / 5

# Train-Test Split
train_index <- createDataPartition(
  model_data$preeclampsia,
  p = 0.8,
  list = FALSE
)

train_data <- model_data[train_index, ]
test_data  <- model_data[-train_index, ]

# Model estimation (using Generalized Logistic Regression)
pe_model <- glm(
  preeclampsia ~ maternal_age_group + bmi_5unit + gestational_diabetes,
  data = train_data,
  family = binomial()
)

#Overview of the models' coefficients, confidence intervals and p-values
summary(pe_model)

# Predictions
pe_probs <- predict(pe_model, newdata = test_data, type = "response")

# ROC-AUC
pe_roc <- roc(test_data$preeclampsia, pe_probs)
print(pe_roc)
print(auc(pe_roc))

# Odds ratios
options(scipen = 999)
exp(cbind(
  OddsRatio = coef(pe_model),
  confint(pe_model)
))

# =========================================================
# 6. Variable importance overview
# =========================================================

cat("\n==============================\n")
cat("Most relevant predictors\n")
cat("==============================\n")

cat("\nGestational Diabetes Model\n")
print(tidy(gdm_model))

cat("
Gestational Hypertension Model
")
print(tidy(ghtn_model))

cat("
Preterm birth Model
")
print(tidy(ptb_model))

cat("
Macrosomia Model
")
print(tidy(macros_model))

cat("
Preeclampsia Model
")
print(tidy(pe_model))



# =========================================================
# 7. Optional ROC plots
# =========================================================

plot(gdm_roc, main = "ROC Curve – Gestational Diabetes")
plot(ghtn_roc, main = "ROC Curve – Gestational Hypertension")
plot(ptb_roc, main = "ROC Curve – Preterm birth")
plot(macros_roc, main = "ROC Curve – Macrosomia")
plot(pe_roc, main = "ROC Curve – Preeclampsia")