# Research Metrics Africa — public-health logistic-regression workflow
# Synthetic teaching example. No participant or client data.
# Place this script beside synthetic-public-health-logistic-regression.csv.

data <- read.csv("synthetic-public-health-logistic-regression.csv")

# Preserve source fields; derive an interpretable 10-year age unit.
data$age_10_years <- data$age_years / 10
data$sex <- factor(data$sex, levels = c("Male", "Female"))
data$current_smoker <- factor(data$current_smoker, levels = c(0, 1), labels = c("No", "Yes"))
data$low_physical_activity <- factor(data$low_physical_activity, levels = c(0, 1), labels = c("No", "Yes"))

# Reconcile fields used in the model before fitting it.
model_fields <- c("elevated_bp", "age_10_years", "sex", "current_smoker", "low_physical_activity")
stopifnot(all(data$elevated_bp %in% c(0, 1)))
print(colSums(is.na(data[model_fields])))

fit <- glm(
  elevated_bp ~ age_10_years + sex + current_smoker + low_physical_activity,
  data = data,
  family = binomial(link = "logit")
)

# Report adjusted odds ratios and Wald confidence intervals.
beta <- coef(fit)
se <- sqrt(diag(vcov(fit)))
results <- data.frame(
  term = names(beta),
  adjusted_odds_ratio = exp(beta),
  lower_95 = exp(beta - 1.96 * se),
  upper_95 = exp(beta + 1.96 * se),
  p_value = summary(fit)$coefficients[, "Pr(>|z|)"]
)
print(results, row.names = FALSE)

# Basic review outputs; a full study needs design-specific diagnostics.
data$predicted_probability <- predict(fit, type = "response")
print(summary(fit))
print(summary(data$predicted_probability))
print(influence.measures(fit))

# Evidence boundary:
# A cross-sectional association does not establish temporal order or causality.