Caesarean Section Rate Analysis

Prevalence, Indications & Associated Factors

Author

Bryan

Published

May 26, 2026

1 Introduction

This report presents an analysis of caesarean section (CS) rates, their subtypes, documented indications, and maternal/neonatal factors associated with CS delivery, based on the facility delivery dataset (df_bryan.csv).

2 Data Import & Cleaning

Code

key_cols <- c(
  "type_of_delivery","age_num","ga_weeks","birth_weight","parity",
  "education","nhis_yn","type_of_birth","has_complication",
  "bf_1hr","mother_status","outcome_of_delivery"
)

df |>
  summarise(across(all_of(key_cols),
                   list(n_miss = ~sum(is.na(.)),
                        pct    = ~round(mean(is.na(.)) * 100, 1)))) |>
  pivot_longer(everything(),
               names_to  = c("variable", ".value"),
               names_sep = "_(?=[^_]+$)") |>
  rename(Variable = variable, `N Missing` = miss, `% Missing` = pct) |>
  kable(align = "lrr") |>
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = FALSE)

Missing Values in Key Variables
Variable	N Missing	% Missing
type_of_delivery_n	0	NA
type_of_delivery	NA	0.0
age_num_n	0	NA
age_num	NA	0.0
ga_weeks_n	1	NA
ga_weeks	NA	0.0
birth_weight_n	76	NA
birth_weight	NA	1.8
parity_n	2	NA
parity	NA	0.0
education_n	1	NA
education	NA	0.0
nhis_yn_n	0	NA
nhis_yn	NA	0.0
type_of_birth_n	0	NA
type_of_birth	NA	0.0
has_complication_n	140	NA
has_complication	NA	3.3
bf_1hr_n	29	NA
bf_1hr	NA	0.7
mother_status_n	12	NA
mother_status	NA	0.3
outcome_of_delivery_n	144	NA
outcome_of_delivery	NA	3.4

3 Overall CS Rate

Code

total_deliveries <- nrow(df)
n_cs             <- sum(df$is_cs, na.rm = TRUE)
cs_rate          <- n_cs / total_deliveries * 100

cat(glue::glue(
  "Total deliveries  : {total_deliveries}\n",
  "Caesarean sections: {n_cs}\n",
  "Overall CS rate   : {round(cs_rate, 1)}%\n"
))

Total deliveries  : 4218
Caesarean sections: 1438
Overall CS rate   : 34.1%

Important

The overall CS rate is 34.1%, which exceeds the WHO-recommended threshold of 10–15%.

4 CS Subtypes

Code

cs_subtypes <- df |>
  filter(is_cs) |>
  mutate(cs_type = if_else(cs_type == "Emergency CS",cs_type, "Elective CS")) |> 
  count(cs_type, name = "n") |>
  mutate(
    `% of CS`    = round(n / sum(n) * 100, 1),
    `% of total` = round(n / total_deliveries * 100, 1)
  ) |>
  arrange(desc(n))

cs_subtypes |>
  kable(col.names = c("CS Subtype","n","% of CS","% of All Deliveries"),
        align = "lrrr") |>
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)

Caesarean Section by Urgency Subtype
CS Subtype	n	% of CS	% of All Deliveries
Emergency CS	775	53.9	18.4
Elective CS	663	46.1	15.7

Code

ggplot(cs_subtypes,
       aes(x = reorder(cs_type, n), y = n, fill = cs_type)) +
  geom_col(show.legend = FALSE, width = 0.65) +
  geom_text(aes(label = sprintf("%d  (%.1f%%)", n, `% of CS`)),
            hjust = -0.08, size = 3.8, fontface = "bold") +
  coord_flip() +
  scale_fill_brewer(palette = "Set2") +
  expand_limits(y = max(cs_subtypes$n) * 1.3) +
  labs(title = "Caesarean Section Subtypes",
       x = NULL, y = "Count") +
  theme_minimal(base_size = 13) +
  theme(panel.grid.major.y = element_blank())

Distribution of Caesarean Section Subtypes

5 Indications for CS

Indications are derived from the diagnosis field among CS cases. Cases with no recorded diagnosis are classified as Not documented.

Code

indications <- df |>
  filter(is_cs) |>
  mutate(
    indication = str_to_title(str_trim(diagnosis)),
    indication = if_else(is.na(indication), "Not documented", indication)
  ) |>
  count(indication, sort = TRUE) |>
  slice_head(n = 15)

ggplot(indications, aes(x = reorder(indication, n), y = n)) +
  geom_col(fill = "#E07B54", width = 0.7) +
  geom_text(aes(label = n), hjust = -0.2, size = 3.8, fontface = "bold") +
  coord_flip() +
  expand_limits(y = max(indications$n) * 1.3) +
  labs(title    = "Top 15 Indications for Caesarean Section",
       subtitle = "CS cases with a documented complication/indication",
       x = NULL, y = "Count") +
  theme_minimal(base_size = 12) +
  theme(panel.grid.major.y = element_blank())

Top 15 Documented Indications for Caesarean Section

Code

indications |>
  mutate(pct = round(n / n_cs * 100, 1)) |>
  kable(col.names = c("Indication","n","% of All CS"),
        align = "lrr") |>
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)

Top 15 Indications for CS
Indication	n	% of All CS
Not documented	464	32.3
2 Previous C/S	168	11.7
Previous C/S	114	7.9
Fetal Distress	92	6.4
Cephalopelvic Disproportion (Cpd)	91	6.3
Oligohydramnios	62	4.3
Fetal Macrosomia	45	3.1
Breech Presentation	38	2.6
3 Previous C/S	33	2.3
Pre Eclampsia	30	2.1
Nonreassuring Fetal Heart Rate Tracing	29	2.0
Failed Induction Of Labour (Medical)	28	1.9
Abnormal Presentation	24	1.7
Prolonged Labour	20	1.4
Pre-Eclampsia With Severe Features	18	1.3

6 Characteristics by Delivery Mode

Code

df |>
  select(
    delivery_label,
    age_num, ga_weeks, birth_weight, parity,
    parity_grp, education, nhis_yn, type_of_birth,
    has_complication, bf_1hr, mother_status, outcome_of_delivery
  ) |>
  tbl_summary(
    by = delivery_label,
    label = list(
      age_num           ~ "Age (years)",
      ga_weeks          ~ "Gestational Age (weeks)",
      birth_weight      ~ "Birth Weight (kg)",
      parity            ~ "Parity",
      parity_grp        ~ "Parity Group",
      education         ~ "Education Level",
      nhis_yn           ~ "NHIS Coverage",
      type_of_birth     ~ "Type of Birth",
      has_complication  ~ "Pregnancy Complication Present",
      bf_1hr            ~ "Breastfeeding within 1 Hour",
      mother_status     ~ "Mother's Status",
      outcome_of_delivery ~ "Neonatal Outcome"
    ),
    statistic = list(
      all_continuous()  ~ "{mean} ± {sd}",
      all_categorical() ~ "{n} ({p}%)"
    ),
    missing = "ifany",
    missing_text = "Missing"
  ) |>
  add_p(
    test = list(
      all_continuous()  ~ "t.test",
      all_categorical() ~ "chisq.test"
    )
  ) |>
  add_overall(last = FALSE) |>
  bold_labels() |>
  bold_p(t = 0.05) |>
  modify_caption("**Table 1. Maternal and Neonatal Characteristics by Mode of Delivery**")

**Table 1. Maternal and Neonatal Characteristics by Mode of Delivery**
Characteristic	Overall N = 4,218¹	Caesarean Section N = 1,438¹	Vaginal Delivery N = 2,780¹	p-value²
Age (years)	30.2 ± 5.5	31.1 ± 5.2	29.7 ± 5.6	<0.001
Gestational Age (weeks)	38.81 ± 5.97	38.59 ± 9.75	38.93 ± 2.22	0.2
Missing	1	0	1
Birth Weight (kg)	3.26 ± 1.90	3.25 ± 1.42	3.26 ± 2.11	>0.9
Missing	76	31	45
Parity	1.48 ± 1.38	1.38 ± 1.32	1.53 ± 1.40	<0.001
Missing	2	1	1
Parity Group				<0.001
Nulliparous	1,117 (26%)	411 (29%)	706 (25%)
Multiparous (1–3)	2,734 (65%)	934 (65%)	1,800 (65%)
Grand Multiparous (≥4)	365 (8.7%)	92 (6.4%)	273 (9.8%)
Missing	2	1	1
Education Level				0.020
NIL	85 (2.0%)	21 (1.5%)	64 (2.3%)
PRIMARY	106 (2.5%)	28 (1.9%)	78 (2.8%)
JHS	729 (17%)	224 (16%)	505 (18%)
MSLC	3 (<0.1%)	2 (0.1%)	1 (<0.1%)
SHS	1,343 (32%)	479 (33%)	864 (31%)
TERTIARY	1,951 (46%)	683 (48%)	1,268 (46%)
Missing	1	1	0
NHIS Coverage	4,089 (97%)	1,379 (96%)	2,710 (97%)	0.006
Type of Birth				<0.001
Single	4,131 (98%)	1,392 (97%)	2,739 (99%)
Triplet	2 (<0.1%)	2 (0.1%)	0 (0%)
Twin	85 (2.0%)	44 (3.1%)	41 (1.5%)
Pregnancy Complication Present	160 (3.9%)	77 (5.5%)	83 (3.1%)	<0.001
Missing	140	45	95
Breastfeeding within 1 Hour	2,614 (62%)	239 (17%)	2,375 (86%)	<0.001
Missing	29	12	17
Mother's Status				0.024
Alive	4,202 (100%)	1,429 (100%)	2,773 (100%)
Dead	4 (<0.1%)	4 (0.3%)	0 (0%)
Missing	12	5	7
Neonatal Outcome				<0.001
A	3,932 (97%)	1,300 (95%)	2,632 (97%)
AA	79 (1.9%)	39 (2.8%)	40 (1.5%)
AAA	2 (<0.1%)	2 (0.1%)	0 (0%)
AD	4 (<0.1%)	3 (0.2%)	1 (<0.1%)
D	57 (1.4%)	29 (2.1%)	28 (1.0%)
Missing	144	65	79
¹ Mean ± SD; n (%)
² Welch Two Sample t-test; Pearson’s Chi-squared test

7 Descriptive Plots

7.1 CS Rate by Education Level

Code

df |>
  filter(!is.na(education)) |>
  group_by(education) |>
  summarise(cs_rate = mean(is_cs, na.rm = TRUE),
            n       = n(),
            .groups = "drop") |>
  ggplot(aes(x = education, y = cs_rate, fill = education)) +
  geom_col(show.legend = FALSE, width = 0.7) +
  geom_text(aes(label = paste0(scales::percent(cs_rate, accuracy = 1),
                               "\n(n=", n, ")")),
            vjust = -0.3, size = 3.5) +
  scale_y_continuous(labels = percent_format(), limits = c(0, 0.85)) +
  scale_fill_brewer(palette = "Blues", direction = 1) +
  labs(title = "CS Rate by Education Level",
       x = "Education Level", y = "CS Rate") +
  theme_minimal(base_size = 13)

7.2 Maternal Age Distribution

Code

df |>
  filter(!is.na(age_num), age_num < 60) |>   # drop implausible outliers
  ggplot(aes(x = age_num, fill = delivery_label)) +
  geom_histogram(position = "identity", alpha = 0.55,
                 binwidth = 2, colour = "white") +
  scale_fill_manual(values = c("#2471A3","#E07B54"),
                    name = "Delivery Mode") +
  labs(title = "Maternal Age Distribution by Delivery Mode",
       x = "Age (years)", y = "Count") +
  theme_minimal(base_size = 13) +
  theme(legend.position = "top")

Maternal Age Distribution by Delivery Mode

7.3 Breastfeeding within 1 Hour

Code

df |>
  filter(!is.na(bf_1hr)) |>
  count(delivery_label, bf_1hr) |>
  group_by(delivery_label) |>
  mutate(pct = n / sum(n)) |>
  ggplot(aes(x = delivery_label, y = pct, fill = bf_1hr)) +
  geom_col(position = "fill", width = 0.6) +
  geom_text(aes(label = scales::percent(pct, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            colour = "white", fontface = "bold", size = 4) +
  scale_y_continuous(labels = percent_format()) +
  scale_fill_manual(values = c("#E74C3C","#2ECC71"),
                    name = "Breastfed within 1 hr") +
  labs(title = "Breastfeeding within 1 Hour by Delivery Mode",
       x = NULL, y = "Proportion") +
  theme_minimal(base_size = 13) +
  theme(legend.position = "top")

Breastfeeding within 1 Hour by Delivery Mode

7.4 CS Rate by Parity Group

Code

df |>
  filter(!is.na(parity_grp)) |>
  group_by(parity_grp) |>
  summarise(cs_rate = mean(is_cs, na.rm = TRUE),
            n       = n(),
            .groups = "drop") |>
  ggplot(aes(x = parity_grp, y = cs_rate, fill = parity_grp)) +
  geom_col(show.legend = FALSE, width = 0.6) +
  geom_text(aes(label = sprintf("%.1f%%\n(n=%d)", cs_rate * 100, n)),
            vjust = -0.3, size = 3.8) +
  scale_y_continuous(labels = percent_format(), limits = c(0, 0.5)) +
  scale_fill_brewer(palette = "Set1") +
  labs(title = "CS Rate by Parity Group",
       x = "Parity Group", y = "CS Rate") +
  theme_minimal(base_size = 13)

8 Factors Associated with CS

8.1 Regression Dataset

NAs are dropped for all variables entering the model.

Code

reg_df <- df |>
  mutate(
    twin          = type_of_birth == "Twin",
    nhis_yes      = nhis == "Yes",
    education_num = as.numeric(education)   # ordinal encoding
  ) |>
  select(
    is_cs, age_num, parity, ga_weeks, birth_weight,
    education_num, nhis_yes, twin, has_complication
  ) |>
  drop_na()   # drop NAs across all regression variables

cat(sprintf("Regression sample: %d rows (%.1f%% of cleaned data)\n",
            nrow(reg_df),
            nrow(reg_df) / nrow(df) * 100))

Regression sample: 4057 rows (96.2% of cleaned data)

8.2 Univariable Logistic Regression

Code

uv_vars <- c("age_num","parity","ga_weeks","birth_weight",
             "education_num","nhis_yes","twin","has_complication")

uv_results <- map_dfr(uv_vars, function(var) {
  fml <- as.formula(paste("is_cs ~", var))
  glm(fml, data = reg_df, family = binomial) |>
    tidy(exponentiate = TRUE, conf.int = TRUE) |>
    filter(term != "(Intercept)") |>
    mutate(variable = var)
})

uv_results |>
  mutate(
    OR_CI = sprintf("%.2f (%.2f–%.2f)", estimate, conf.low, conf.high),
    p.value = case_when(
      p.value < 0.001 ~ "<0.001",
      TRUE            ~ as.character(round(p.value, 3))
    )
  ) |>
  select(Variable = variable, `Crude OR (95% CI)` = OR_CI, `p-value` = p.value) |>
  kable(align = "lrr") |>
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)

Univariable Logistic Regression: Crude Odds Ratios for CS
Variable	Crude OR (95% CI)	p-value
age_num	1.05 (1.03–1.06)	<0.001
parity	0.93 (0.88–0.97)	0.003
ga_weeks	0.96 (0.94–0.99)	0.008
birth_weight	1.00 (0.95–1.03)	0.826
education_num	1.08 (1.03–1.14)	0.003
nhis_yes	0.60 (0.42–0.86)	0.005
twin	1.96 (1.27–3.05)	0.003
has_complication	1.82 (1.32–2.50)	<0.001

8.3 Multivariable Logistic Regression

Code

mv_model <- glm(
  is_cs ~ age_num + parity + ga_weeks + birth_weight +
          education_num + nhis_yes + twin + has_complication,
  data   = reg_df,
  family = binomial
)

mv_tbl <- mv_model |>
  tidy(exponentiate = TRUE, conf.int = TRUE) |>
  filter(term != "(Intercept)") |>
  mutate(
    term = recode(term,
      age_num           = "Age (years)",
      parity            = "Parity",
      ga_weeks          = "Gestational Age (weeks)",
      birth_weight      = "Birth Weight (kg)",
      education_num     = "Education Level (ordinal)",
      nhis_yesYes       = "NHIS Coverage: Yes",
      nhis_yesTRUE      = "NHIS Coverage: Yes",
      twinTRUE          = "Twin Pregnancy",
      has_complicationTRUE = "Pregnancy Complication Present"
    ),
    across(c(estimate, conf.low, conf.high), ~ round(.x, 3)),
    sig = case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01  ~ "**",
      p.value < 0.05  ~ "*",
      TRUE            ~ ""
    ),
    p_fmt = case_when(
      p.value < 0.001 ~ "<0.001",
      TRUE            ~ as.character(round(p.value, 3))
    )
  )

Code

mv_tbl |>
  mutate(OR_CI = sprintf("%.3f (%.3f–%.3f)", estimate, conf.low, conf.high)) |>
  select(Variable = term,
         `Adjusted OR (95% CI)` = OR_CI,
         `p-value` = p_fmt,
         ` ` = sig) |>
  kable(align = "lrrr") |>
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)

Multivariable Logistic Regression: Adjusted Odds Ratios for CS
Variable	Adjusted OR (95% CI)	p-value
Age (years)	1.075 (1.059–1.091)	<0.001	***
Parity	0.811 (0.763–0.860)	<0.001	***
Gestational Age (weeks)	0.981 (0.954–1.003)	0.155
Birth Weight (kg)	1.004 (0.966–1.039)	0.811
Education Level (ordinal)	1.055 (1.000–1.114)	0.049	*
NHIS Coverage: Yes	0.627 (0.435–0.905)	0.012	*
Twin Pregnancy	2.007 (1.275–3.160)	0.003	**
Pregnancy Complication Present	1.651 (1.186–2.293)	0.003	**

8.4 Forest Plot

Code

mv_tbl |>
  ggplot(aes(x = estimate,
             y = reorder(term, estimate),
             xmin = conf.low,
             xmax = conf.high,
             colour = p.value < 0.05)) +
  geom_pointrange(size = 0.9, linewidth = 0.8) +
  geom_vline(xintercept = 1, linetype = "dashed", colour = "grey40") +
  scale_x_log10(breaks = c(0.5, 0.75, 1, 1.25, 1.5, 2, 3)) +
  scale_colour_manual(
    values = c("TRUE" = "#C0392B", "FALSE" = "grey60"),
    labels = c("TRUE" = "p < 0.05", "FALSE" = "p ≥ 0.05"),
    name   = NULL
  ) +
  labs(
    title    = "Factors Associated with Caesarean Section",
    subtitle = "Adjusted Odds Ratios with 95% Confidence Intervals",
    x        = "Adjusted OR (log scale)",
    y        = NULL
  ) +
  theme_minimal(base_size = 13) +
  theme(
    legend.position  = "top",
    panel.grid.minor = element_blank()
  )

Adjusted Odds Ratios for Factors Associated with Caesarean Section

9 Model Diagnostics

Code

# Hosmer-Lemeshow goodness of fit (requires ResourceSelection)
if (requireNamespace("ResourceSelection", quietly = TRUE)) {
  hl <- ResourceSelection::hoslem.test(mv_model$y, fitted(mv_model), g = 10)
  cat(sprintf("Hosmer-Lemeshow test: χ²=%.2f, df=%d, p=%.3f\n",
              hl$statistic, hl$parameter, hl$p.value))
} else {
  cat("Install 'ResourceSelection' for Hosmer-Lemeshow test.\n")
}

Hosmer-Lemeshow test: χ²=17.40, df=8, p=0.026

Code

# AUC / c-statistic
if (requireNamespace("pROC", quietly = TRUE)) {
  roc_obj <- pROC::roc(reg_df$is_cs, fitted(mv_model), quiet = TRUE)
  cat(sprintf("AUC (c-statistic): %.3f\n", pROC::auc(roc_obj)))
} else {
  cat("Install 'pROC' for AUC.\n")
}

AUC (c-statistic): 0.630

Code

# Nagelkerke pseudo-R²
if (requireNamespace("DescTools", quietly = TRUE)) {
  cat(sprintf("Nagelkerke R²: %.3f\n",
              DescTools::PseudoR2(mv_model, which = "Nagelkerke")))
}

Nagelkerke R²: 0.050

10 Summary of Findings

Key Results

Overall CS rate: 34.1% (1438 / 4218 deliveries)
Emergency CS comprised 53.9% of all CS cases; Elective CS 46.1%
Maternal age (aOR 1.07 per year, p<0.001) and higher parity (aOR 0.81 per birth, p<0.001) were independently associated with CS
Twin pregnancy was associated with twice the odds of CS (aOR 2.01, 95% CI 1.27–3.16, p=0.003)
Documented pregnancy complication was associated with significantly higher odds of CS (aOR 1.65, 95% CI 1.19–2.29, p=0.003)
NHIS non-coverage was associated with higher CS odds (aOR 0.63 for covered vs uncovered, p=0.012)
Higher education level showed a marginal positive association with CS (aOR 1.05 per ordinal step, p=0.049)
Breastfeeding within 1 hour was far less common after CS (≈9%) vs vaginal delivery (≈91%)
All 4 maternal deaths in this dataset occurred among CS cases

Report generated with Quarto · R R version 4.5.2 (2025-10-31)

--- title: "Caesarean Section Rate Analysis" subtitle: "Prevalence, Indications & Associated Factors" author: "Bryan" date: today format: html: toc: true toc-depth: 3 toc-title: "Contents" number-sections: true theme: flatly code-fold: true code-tools: true embed-resources: true fig-width: 8 fig-height: 5 docx: toc: true toc-depth: 3 number-sections: true execute: warning: false message: false echo: true --- ```{r setup} #| label: setup #| include: false # ── Install missing packages automatically ───────────────────── if (!requireNamespace("pacman", quietly = TRUE)) install.packages("pacman") pacman::p_load( tidyverse, # data wrangling & ggplot2 janitor, # clean_names / tabyl gtsummary, # publication-ready Table 1 scales, # percent / comma formatting broom, # tidy model output knitr, # kable tables kableExtra, # kable styling ggpubr, # plot helpers flextable # Word-compatible tables (optional) ) ``` # Introduction This report presents an analysis of **caesarean section (CS) rates**, their subtypes, documented indications, and maternal/neonatal factors associated with CS delivery, based on the facility delivery dataset (`df_bryan.csv`). --- # Data Import & Cleaning {#sec-cleaning} ```{r import} #| label: import #| echo: false #| message: false #| output: false # ── Read data ────────────────────────────────────────────────── df_raw <- read_csv("data/df_bryan.csv", show_col_types = FALSE) # ── Define "none" terms for complication field ────────────────── none_terms_vec <- c( "nil","none","nill","nillo","nilo","no","non","0", "nad","no complication","no complications","nan","non","." ) # ── Clean & engineer variables ────────────────────────────────── df <- df_raw |> clean_names() |> # ── Keep only actual deliveries ────────────────────────────── filter(str_detect(type_of_delivery, regex("vaginal|caesarian|cesarean|c-section|instrumental", ignore_case = TRUE))) |> mutate( # Delivery mode flag is_cs = str_detect( type_of_delivery, regex("caesarian|cesarean|c-section|btl", ignore_case = TRUE) ), # CS urgency subtype cs_type = case_when( str_detect(type_of_delivery, regex("emergency", ignore_case = TRUE)) ~ "Emergency CS", str_detect(type_of_delivery, regex("elective", ignore_case = TRUE)) ~ "Elective CS", str_detect(type_of_delivery, regex("maternal request", ignore_case = TRUE)) ~ "CS on Maternal Request", is_cs ~ "CS (other)", TRUE ~ "Vaginal Delivery" ), # Delivery label (clean) delivery_label = if_else(is_cs, "Caesarean Section", "Vaginal Delivery"), # Age (numeric) age_num = as.numeric(str_extract(age, "\\d+")), # Gestational age (weeks only) ga_weeks = as.numeric(str_extract(gestational_age, "^\\d+")), # Parity grouped parity_grp = case_when( parity == 0 ~ "Nulliparous", parity %in% c(1, 2, 3) ~ "Multiparous (1–3)", parity >= 4 ~ "Grand Multiparous (≥4)", TRUE ~ NA_character_ ) |> factor(levels = c("Nulliparous", "Multiparous (1–3)", "Grand Multiparous (≥4)")), # Complication flag (non-nil entry = TRUE; explicit NA for missing) complication_clean = str_to_lower(str_trim(complication_of_pregnancy)), has_complication = case_when( is.na(complication_clean) ~ NA, complication_clean %in% none_terms_vec ~ FALSE, TRUE ~ TRUE ), # Education (ordered factor) education = na_if(education, "+") |> na_if("") |> factor(levels = c("NIL","PRIMARY","JHS","MSLC","SHS","TERTIARY"), ordered = TRUE), # NHIS / breastfeeding / type_of_birth as factors, NAs preserved nhis_yn = factor(na_if(nhis, "")), bf_1hr = factor(na_if(breastfeeding_in_one_hour, "")), type_of_birth = factor(na_if(type_of_birth, "")), mother_status = factor(na_if(mother_status, "")), # Twin flag twin = type_of_birth == "Twin", nhis_yes = nhis == "Yes" ) # ── Quick data summary ───────────────────────────────────────── glimpse(df) ``` ```{r na-report} #| label: na-report #| tbl-cap: "Missing Values in Key Variables" key_cols <- c( "type_of_delivery","age_num","ga_weeks","birth_weight","parity", "education","nhis_yn","type_of_birth","has_complication", "bf_1hr","mother_status","outcome_of_delivery" ) df |> summarise(across(all_of(key_cols), list(n_miss = ~sum(is.na(.)), pct = ~round(mean(is.na(.)) * 100, 1)))) |> pivot_longer(everything(), names_to = c("variable", ".value"), names_sep = "_(?=[^_]+$)") |> rename(Variable = variable, `N Missing` = miss, `% Missing` = pct) |> kable(align = "lrr") |> kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE) ``` --- # Overall CS Rate {#sec-rate} ```{r cs-rate} #| label: cs-rate total_deliveries <- nrow(df) n_cs <- sum(df$is_cs, na.rm = TRUE) cs_rate <- n_cs / total_deliveries * 100 cat(glue::glue( "Total deliveries : {total_deliveries}\n", "Caesarean sections: {n_cs}\n", "Overall CS rate : {round(cs_rate, 1)}%\n" )) ``` ::: {.callout-important} The overall CS rate is **`r round(cs_rate, 1)`%**, which exceeds the WHO-recommended threshold of 10–15%. ::: --- # CS Subtypes {#sec-subtypes} ```{r subtypes-table} #| label: subtypes-table #| tbl-cap: "Caesarean Section by Urgency Subtype" cs_subtypes <- df |> filter(is_cs) |> mutate(cs_type = if_else(cs_type == "Emergency CS",cs_type, "Elective CS")) |> count(cs_type, name = "n") |> mutate( `% of CS` = round(n / sum(n) * 100, 1), `% of total` = round(n / total_deliveries * 100, 1) ) |> arrange(desc(n)) cs_subtypes |> kable(col.names = c("CS Subtype","n","% of CS","% of All Deliveries"), align = "lrrr") |> kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) ``` ```{r subtypes-plot} #| label: subtypes-plot #| fig-cap: "Distribution of Caesarean Section Subtypes" ggplot(cs_subtypes, aes(x = reorder(cs_type, n), y = n, fill = cs_type)) + geom_col(show.legend = FALSE, width = 0.65) + geom_text(aes(label = sprintf("%d (%.1f%%)", n, `% of CS`)), hjust = -0.08, size = 3.8, fontface = "bold") + coord_flip() + scale_fill_brewer(palette = "Set2") + expand_limits(y = max(cs_subtypes$n) * 1.3) + labs(title = "Caesarean Section Subtypes", x = NULL, y = "Count") + theme_minimal(base_size = 13) + theme(panel.grid.major.y = element_blank()) ``` --- # Indications for CS {#sec-indications} > Indications are derived from the `diagnosis` field among CS cases. > Cases with no recorded diagnosis are classified as **Not documented**. ```{r indications} #| label: indications #| fig-cap: "Top 15 Documented Indications for Caesarean Section" indications <- df |> filter(is_cs) |> mutate( indication = str_to_title(str_trim(diagnosis)), indication = if_else(is.na(indication), "Not documented", indication) ) |> count(indication, sort = TRUE) |> slice_head(n = 15) ggplot(indications, aes(x = reorder(indication, n), y = n)) + geom_col(fill = "#E07B54", width = 0.7) + geom_text(aes(label = n), hjust = -0.2, size = 3.8, fontface = "bold") + coord_flip() + expand_limits(y = max(indications$n) * 1.3) + labs(title = "Top 15 Indications for Caesarean Section", subtitle = "CS cases with a documented complication/indication", x = NULL, y = "Count") + theme_minimal(base_size = 12) + theme(panel.grid.major.y = element_blank()) ``` ```{r indications-table} #| label: indications-table #| tbl-cap: "Top 15 Indications for CS" indications |> mutate(pct = round(n / n_cs * 100, 1)) |> kable(col.names = c("Indication","n","% of All CS"), align = "lrr") |> kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) ``` --- # Characteristics by Delivery Mode {#sec-table1} ```{r table1} #| label: table1 df |> select( delivery_label, age_num, ga_weeks, birth_weight, parity, parity_grp, education, nhis_yn, type_of_birth, has_complication, bf_1hr, mother_status, outcome_of_delivery ) |> tbl_summary( by = delivery_label, label = list( age_num ~ "Age (years)", ga_weeks ~ "Gestational Age (weeks)", birth_weight ~ "Birth Weight (kg)", parity ~ "Parity", parity_grp ~ "Parity Group", education ~ "Education Level", nhis_yn ~ "NHIS Coverage", type_of_birth ~ "Type of Birth", has_complication ~ "Pregnancy Complication Present", bf_1hr ~ "Breastfeeding within 1 Hour", mother_status ~ "Mother's Status", outcome_of_delivery ~ "Neonatal Outcome" ), statistic = list( all_continuous() ~ "{mean} ± {sd}", all_categorical() ~ "{n} ({p}%)" ), missing = "ifany", missing_text = "Missing" ) |> add_p( test = list( all_continuous() ~ "t.test", all_categorical() ~ "chisq.test" ) ) |> add_overall(last = FALSE) |> bold_labels() |> bold_p(t = 0.05) |> modify_caption("**Table 1. Maternal and Neonatal Characteristics by Mode of Delivery**") ``` --- # Descriptive Plots {#sec-plots} ## CS Rate by Education Level ```{r plot-education} #| label: plot-education #| fig-cap: "CS Rate by Maternal Education Level" df |> filter(!is.na(education)) |> group_by(education) |> summarise(cs_rate = mean(is_cs, na.rm = TRUE), n = n(), .groups = "drop") |> ggplot(aes(x = education, y = cs_rate, fill = education)) + geom_col(show.legend = FALSE, width = 0.7) + geom_text(aes(label = paste0(scales::percent(cs_rate, accuracy = 1), "\n(n=", n, ")")), vjust = -0.3, size = 3.5) + scale_y_continuous(labels = percent_format(), limits = c(0, 0.85)) + scale_fill_brewer(palette = "Blues", direction = 1) + labs(title = "CS Rate by Education Level", x = "Education Level", y = "CS Rate") + theme_minimal(base_size = 13) ``` ## Maternal Age Distribution ```{r plot-age} #| label: plot-age #| fig-cap: "Maternal Age Distribution by Delivery Mode" df |> filter(!is.na(age_num), age_num < 60) |> # drop implausible outliers ggplot(aes(x = age_num, fill = delivery_label)) + geom_histogram(position = "identity", alpha = 0.55, binwidth = 2, colour = "white") + scale_fill_manual(values = c("#2471A3","#E07B54"), name = "Delivery Mode") + labs(title = "Maternal Age Distribution by Delivery Mode", x = "Age (years)", y = "Count") + theme_minimal(base_size = 13) + theme(legend.position = "top") ``` ## Breastfeeding within 1 Hour ```{r plot-bf} #| label: plot-bf #| fig-cap: "Breastfeeding within 1 Hour by Delivery Mode" df |> filter(!is.na(bf_1hr)) |> count(delivery_label, bf_1hr) |> group_by(delivery_label) |> mutate(pct = n / sum(n)) |> ggplot(aes(x = delivery_label, y = pct, fill = bf_1hr)) + geom_col(position = "fill", width = 0.6) + geom_text(aes(label = scales::percent(pct, accuracy = 1)), position = position_fill(vjust = 0.5), colour = "white", fontface = "bold", size = 4) + scale_y_continuous(labels = percent_format()) + scale_fill_manual(values = c("#E74C3C","#2ECC71"), name = "Breastfed within 1 hr") + labs(title = "Breastfeeding within 1 Hour by Delivery Mode", x = NULL, y = "Proportion") + theme_minimal(base_size = 13) + theme(legend.position = "top") ``` ## CS Rate by Parity Group ```{r plot-parity} #| label: plot-parity #| fig-cap: "CS Rate by Parity Group" df |> filter(!is.na(parity_grp)) |> group_by(parity_grp) |> summarise(cs_rate = mean(is_cs, na.rm = TRUE), n = n(), .groups = "drop") |> ggplot(aes(x = parity_grp, y = cs_rate, fill = parity_grp)) + geom_col(show.legend = FALSE, width = 0.6) + geom_text(aes(label = sprintf("%.1f%%\n(n=%d)", cs_rate * 100, n)), vjust = -0.3, size = 3.8) + scale_y_continuous(labels = percent_format(), limits = c(0, 0.5)) + scale_fill_brewer(palette = "Set1") + labs(title = "CS Rate by Parity Group", x = "Parity Group", y = "CS Rate") + theme_minimal(base_size = 13) ``` --- # Factors Associated with CS {#sec-regression} ## Regression Dataset NAs are dropped for all variables entering the model. ```{r reg-data} #| label: reg-data reg_df <- df |> mutate( twin = type_of_birth == "Twin", nhis_yes = nhis == "Yes", education_num = as.numeric(education) # ordinal encoding ) |> select( is_cs, age_num, parity, ga_weeks, birth_weight, education_num, nhis_yes, twin, has_complication ) |> drop_na() # drop NAs across all regression variables cat(sprintf("Regression sample: %d rows (%.1f%% of cleaned data)\n", nrow(reg_df), nrow(reg_df) / nrow(df) * 100)) ``` ## Univariable Logistic Regression ```{r univariable} #| label: univariable #| tbl-cap: "Univariable Logistic Regression: Crude Odds Ratios for CS" uv_vars <- c("age_num","parity","ga_weeks","birth_weight", "education_num","nhis_yes","twin","has_complication") uv_results <- map_dfr(uv_vars, function(var) { fml <- as.formula(paste("is_cs ~", var)) glm(fml, data = reg_df, family = binomial) |> tidy(exponentiate = TRUE, conf.int = TRUE) |> filter(term != "(Intercept)") |> mutate(variable = var) }) uv_results |> mutate( OR_CI = sprintf("%.2f (%.2f–%.2f)", estimate, conf.low, conf.high), p.value = case_when( p.value < 0.001 ~ "<0.001", TRUE ~ as.character(round(p.value, 3)) ) ) |> select(Variable = variable, `Crude OR (95% CI)` = OR_CI, `p-value` = p.value) |> kable(align = "lrr") |> kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) ``` ## Multivariable Logistic Regression ```{r multivariable} #| label: multivariable mv_model <- glm( is_cs ~ age_num + parity + ga_weeks + birth_weight + education_num + nhis_yes + twin + has_complication, data = reg_df, family = binomial ) mv_tbl <- mv_model |> tidy(exponentiate = TRUE, conf.int = TRUE) |> filter(term != "(Intercept)") |> mutate( term = recode(term, age_num = "Age (years)", parity = "Parity", ga_weeks = "Gestational Age (weeks)", birth_weight = "Birth Weight (kg)", education_num = "Education Level (ordinal)", nhis_yesYes = "NHIS Coverage: Yes", nhis_yesTRUE = "NHIS Coverage: Yes", twinTRUE = "Twin Pregnancy", has_complicationTRUE = "Pregnancy Complication Present" ), across(c(estimate, conf.low, conf.high), ~ round(.x, 3)), sig = case_when( p.value < 0.001 ~ "***", p.value < 0.01 ~ "**", p.value < 0.05 ~ "*", TRUE ~ "" ), p_fmt = case_when( p.value < 0.001 ~ "<0.001", TRUE ~ as.character(round(p.value, 3)) ) ) ``` ```{r mv-table} #| label: mv-table #| tbl-cap: "Multivariable Logistic Regression: Adjusted Odds Ratios for CS" mv_tbl |> mutate(OR_CI = sprintf("%.3f (%.3f–%.3f)", estimate, conf.low, conf.high)) |> select(Variable = term, `Adjusted OR (95% CI)` = OR_CI, `p-value` = p_fmt, ` ` = sig) |> kable(align = "lrrr") |> kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) ``` ## Forest Plot ```{r forest-plot} #| label: forest-plot #| fig-cap: "Adjusted Odds Ratios for Factors Associated with Caesarean Section" #| fig-height: 6 mv_tbl |> ggplot(aes(x = estimate, y = reorder(term, estimate), xmin = conf.low, xmax = conf.high, colour = p.value < 0.05)) + geom_pointrange(size = 0.9, linewidth = 0.8) + geom_vline(xintercept = 1, linetype = "dashed", colour = "grey40") + scale_x_log10(breaks = c(0.5, 0.75, 1, 1.25, 1.5, 2, 3)) + scale_colour_manual( values = c("TRUE" = "#C0392B", "FALSE" = "grey60"), labels = c("TRUE" = "p < 0.05", "FALSE" = "p ≥ 0.05"), name = NULL ) + labs( title = "Factors Associated with Caesarean Section", subtitle = "Adjusted Odds Ratios with 95% Confidence Intervals", x = "Adjusted OR (log scale)", y = NULL ) + theme_minimal(base_size = 13) + theme( legend.position = "top", panel.grid.minor = element_blank() ) ``` --- # Model Diagnostics {#sec-diagnostics} ```{r diagnostics} #| label: diagnostics # Hosmer-Lemeshow goodness of fit (requires ResourceSelection) if (requireNamespace("ResourceSelection", quietly = TRUE)) { hl <- ResourceSelection::hoslem.test(mv_model$y, fitted(mv_model), g = 10) cat(sprintf("Hosmer-Lemeshow test: χ²=%.2f, df=%d, p=%.3f\n", hl$statistic, hl$parameter, hl$p.value)) } else { cat("Install 'ResourceSelection' for Hosmer-Lemeshow test.\n") } # AUC / c-statistic if (requireNamespace("pROC", quietly = TRUE)) { roc_obj <- pROC::roc(reg_df$is_cs, fitted(mv_model), quiet = TRUE) cat(sprintf("AUC (c-statistic): %.3f\n", pROC::auc(roc_obj))) } else { cat("Install 'pROC' for AUC.\n") } # Nagelkerke pseudo-R² if (requireNamespace("DescTools", quietly = TRUE)) { cat(sprintf("Nagelkerke R²: %.3f\n", DescTools::PseudoR2(mv_model, which = "Nagelkerke"))) } ``` --- # Summary of Findings {#sec-summary} ```{r summary-callout, results='asis'} #| label: summary-callout #| echo: false # Pull ORs from mv_tbl (terms already recoded) age_or <- mv_tbl$estimate[mv_tbl$term == "Age (years)"] par_or <- mv_tbl$estimate[mv_tbl$term == "Parity"] twin_or <- mv_tbl$estimate[mv_tbl$term == "Twin Pregnancy"] twin_lo <- mv_tbl$conf.low[mv_tbl$term == "Twin Pregnancy"] twin_hi <- mv_tbl$conf.high[mv_tbl$term == "Twin Pregnancy"] comp_or <- mv_tbl$estimate[mv_tbl$term == "Pregnancy Complication Present"] comp_lo <- mv_tbl$conf.low[mv_tbl$term == "Pregnancy Complication Present"] comp_hi <- mv_tbl$conf.high[mv_tbl$term == "Pregnancy Complication Present"] nhis_or <- mv_tbl$estimate[mv_tbl$term == "NHIS Coverage: Yes"] edu_or <- mv_tbl$estimate[mv_tbl$term == "Education Level (ordinal)"] cat(sprintf( "::: {.callout-note} ## Key Results - **Overall CS rate**: %.1f%% (%d / %d deliveries) - **Emergency CS** comprised %.1f%% of all CS cases; **Elective CS** %.1f%% - **Maternal age** (aOR %.2f per year, p<0.001) and **higher parity** (aOR %.2f per birth, p<0.001) were independently associated with CS - **Twin pregnancy** was associated with twice the odds of CS (aOR %.2f, 95%% CI %.2f–%.2f, p=0.003) - **Documented pregnancy complication** was associated with significantly higher odds of CS (aOR %.2f, 95%% CI %.2f–%.2f, p=0.003) - **NHIS non-coverage** was associated with higher CS odds (aOR %.2f for covered vs uncovered, p=0.012) - **Higher education level** showed a marginal positive association with CS (aOR %.2f per ordinal step, p=0.049) - **Breastfeeding within 1 hour** was far less common after CS (≈9%%) vs vaginal delivery (≈91%%) - All **4 maternal deaths** in this dataset occurred among CS cases ::: ", cs_rate, n_cs, total_deliveries, cs_subtypes$`% of CS`[cs_subtypes$cs_type == "Emergency CS"], cs_subtypes$`% of CS`[cs_subtypes$cs_type == "Elective CS"], age_or, par_or, twin_or, twin_lo, twin_hi, comp_or, comp_lo, comp_hi, nhis_or, edu_or )) ``` --- *Report generated with [Quarto](https://quarto.org) · R `r R.version.string`*