Caesarean Section Rate Analysis

Prevalence, Indications & Associated Factors

Author

Bryan

Published

April 4, 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 
# ── 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() |>
  # ── Filter rows with no delivery type ──────────────────────────
  filter(!is.na(type_of_delivery)) |>
  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)
    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, "") |>
      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)
Rows: 4,392
Columns: 56
$ patient_name              <chr> "AWINEBOYA MARTINA", "HOWARD BETTY", "ASANTE…
$ occupation                <chr> "SECRETARY", "TRADING", "TRADING", "TRADING"…
$ gestational_age           <chr> "38+5", "40+1", "36+6", "40+4", "38+6", "41"…
$ admission_date            <dttm> 2023-01-01, 2022-12-31, 2023-01-01, 2023-01…
$ delivery_date             <dttm> 2023-01-01, 2023-01-01, 2023-01-01, 2023-01…
$ type_of_birth             <fct> Single, Single, Twin, Single, Single, Single…
$ breastfeeding_in_one_hour <chr> "Yes", "Yes", "No", "Yes", "Yes", "Yes", "No…
$ cord_care                 <chr> "M. Spirit", "M. Spirit", "M. Spirit", "M. S…
$ female_child              <dbl> 1, 1, 2, NA, 1, 1, 0, NA, 1, 1, NA, 0, 1, 1,…
$ male_child                <dbl> 0, 0, 0, NA, 0, 0, 1, NA, 0, 0, NA, NA, 0, 0…
$ female_alive              <dbl> 1, 1, 2, NA, 1, 1, 0, NA, 1, 1, NA, 0, 1, 1,…
$ male_alive                <dbl> 0, 0, 0, 1, 0, 0, 1, NA, 0, 0, NA, NA, 0, 0,…
$ female_dead               <dbl> 0, 0, 0, 0, 0, 0, 0, NA, 0, 0, NA, 0, 0, 0, …
$ male_dead                 <dbl> 0, 0, 0, NA, 0, 0, 0, NA, 0, 0, NA, NA, 0, 0…
$ female_macerated          <dbl> 0, 0, 0, 0, 0, 0, 0, NA, 0, 0, NA, 0, 0, 0, …
$ male_macerated            <dbl> 0, 0, 0, NA, 0, 0, 0, NA, 0, 0, NA, NA, 0, 0…
$ female_fresh              <dbl> 0, 0, 0, 0, 0, 0, 0, NA, 0, 0, NA, 0, 0, 0, …
$ male_fresh                <dbl> 0, 0, 0, NA, 0, 0, 0, NA, 0, 0, NA, NA, 0, 0…
$ address                   <chr> NA, "APRADE", "PLOT7BLOCKB", "STADIUM", "AWO…
$ education                 <ord> TERTIARY, JHS, SHS, SHS, TERTIARY, TERTIARY,…
$ discharge_date            <dttm> 2023-01-02, 2023-01-02, 2023-01-04, 2023-01…
$ mother_status             <fct> Alive, Alive, Alive, Alive, Alive, Alive, Al…
$ pp_vita                   <chr> "No", NA, NA, "No", NA, NA, "No", NA, NA, NA…
$ ppfp                      <chr> "No", NA, NA, "No", NA, NA, "No", NA, NA, NA…
$ estimated_blood_loss      <dbl> 100, 200, 500, 150, 200, 200, 400, 150, 400,…
$ symp_fundal_height        <dbl> 38, NA, 14, 37, NA, 40, 36, 37, 14, 14, 31, …
$ age                       <chr> "38 Year(s)", "43 Year(s)", "33 Year(s)", "3…
$ complication_of_pregnancy <chr> "none", "none", "nil", "none", "none", "none…
$ folder_no                 <chr> "0034/23", "1870/22", "1455/22", "0570/22", …
$ parity                    <dbl> 1, 4, 0, 1, 0, 1, 5, 3, 0, 0, 1, 2, 2, 0, 1,…
$ nhis                      <chr> "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Y…
$ type_of_delivery          <chr> "Spontaneous Vaginal Delivery without Episio…
$ outcome_of_delivery       <chr> "A", "A", "AA", "A", "A", "A", "A", NA, "A",…
$ no_of_children            <dbl> 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1,…
$ birth_weight              <dbl> 2.9, 3.6, 2.2, 2.9, 2.9, 3.8, 2.5, 2.9, 2.9,…
$ eye_care                  <chr> "Tetracycline", "Tetracycline", "Tetracyclin…
$ reg_date                  <dttm> NA, NA, NA, NA, 2023-11-27, NA, NA, NA, 202…
$ gravida                   <dbl> NA, NA, NA, NA, 2, NA, NA, NA, 3, 3, NA, NA,…
$ g_age_wks                 <dbl> NA, NA, NA, NA, 13, NA, NA, NA, 22, 38, NA, …
$ t_code                    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ a_start_time              <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ a_end_time                <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ diagnosis                 <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ surgeon                   <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ is_cs                     <lgl> FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, TRU…
$ cs_type                   <chr> "Vaginal Delivery", "Vaginal Delivery", "Eme…
$ delivery_label            <chr> "Vaginal Delivery", "Vaginal Delivery", "Cae…
$ age_num                   <dbl> 38, 43, 33, 38, 30, 28, 46, 42, 30, 30, 28, …
$ ga_weeks                  <dbl> 38, 40, 36, 40, 38, 41, 40, 39, 40, 40, 36, …
$ parity_grp                <fct> Multiparous (1–3), Grand Multiparous (≥4), N…
$ complication_clean        <chr> "none", "none", "nil", "none", "none", "none…
$ has_complication          <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FA…
$ nhis_yn                   <fct> Yes, Yes, Yes, Yes, Yes, Yes, Yes, Yes, Yes,…
$ bf_1hr                    <fct> Yes, Yes, No, Yes, Yes, Yes, No, Yes, Yes, Y…
$ twin                      <lgl> FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, FAL…
$ nhis_yes                  <lgl> TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TR…
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.7
parity_n 2 NA
parity NA 0.0
education_n 2 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.2
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.3

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  : 4392
Caesarean sections: 1438
Overall CS rate   : 32.7%
Important

The overall CS rate is 32.7%, which exceeds the WHO-recommended threshold of 10–15%, consistent with a high-volume tertiary/referral facility.

4 CS Subtypes

Code 
cs_subtypes <- df |>
  filter(is_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 17.6
Elective CS 645 44.9 14.7
CS on Maternal Request 16 1.1 0.4
CS (other) 2 0.1 0.0
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 complication_of_pregnancy field among CS cases. Entries recorded as “nil / none / no / NAD” are classified as no complication documented.

Code 
indications <- df |>
  filter(is_cs, !is.na(has_complication), has_complication == TRUE) |>
  mutate(indication = str_to_title(str_trim(complication_of_pregnancy))) |>
  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
Pre Eclampsia 7 0.5
Prom 6 0.4
Gdm 5 0.3
Hpt 5 0.3
Pih 5 0.3
Abruptio 3 0.2
Aph 3 0.2
Dm 3 0.2
Oligohydramnios 3 0.2
N/A 2 0.1
Pre-Eclampsia 2 0.1
Abruptio Placentae 1 0.1
Anaemia 1 0.1
Anaemia In Pregn 1 0.1
Aneamia 1 0.1

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,3921		 Caesarean Section
N = 1,4381		 Vaginal Delivery
N = 2,9541		 p-value2
Age (years) 30.2 ± 5.5 31.1 ± 5.2 29.8 ± 5.6 <0.001
Gestational Age (weeks) 38.78 ± 5.89 38.59 ± 9.75 38.87 ± 2.29 0.3
    Missing 1 0 1
Birth Weight (kg) 3.47 ± 15.06 3.25 ± 1.42 3.58 ± 18.32 0.3
    Missing 76 31 45
Parity 1.48 ± 1.37 1.38 ± 1.32 1.53 ± 1.39 <0.001
    Missing 2 1 1
Parity Group


<0.001
    Nulliparous 1,146 (26%) 411 (29%) 735 (25%)
    Multiparous (1–3) 2,868 (65%) 934 (65%) 1,934 (65%)
    Grand Multiparous (≥4) 376 (8.6%) 92 (6.4%) 284 (9.6%)
    Missing 2 1 1
Education Level


0.023
    NIL 87 (2.0%) 21 (1.5%) 66 (2.2%)
    PRIMARY 117 (2.7%) 28 (1.9%) 89 (3.0%)
    JHS 747 (17%) 224 (16%) 523 (18%)
    MSLC 3 (<0.1%) 2 (0.1%) 1 (<0.1%)
    SHS 1,411 (32%) 479 (33%) 932 (32%)
    TERTIARY 2,025 (46%) 683 (48%) 1,342 (45%)
    Missing 2 1 1
NHIS Coverage 4,252 (97%) 1,379 (96%) 2,873 (97%) 0.020
Type of Birth


0.001
    Single 4,297 (98%) 1,392 (97%) 2,905 (98%)
    Triplet 2 (<0.1%) 2 (0.1%) 0 (0%)
    Twin 93 (2.1%) 44 (3.1%) 49 (1.7%)
Pregnancy Complication Present 176 (4.1%) 77 (5.5%) 99 (3.5%) 0.002
    Missing 140 45 95
Breastfeeding within 1 Hour 2,657 (61%) 239 (17%) 2,418 (82%) <0.001
    Missing 29 12 17
Mother's Status


0.019
    Alive 4,376 (100%) 1,429 (100%) 2,947 (100%)
    Dead 4 (<0.1%) 4 (0.3%) 0 (0%)
    Missing 12 5 7
Neonatal Outcome


<0.001
    A 4,094 (96%) 1,300 (95%) 2,794 (97%)
    AA 87 (2.0%) 39 (2.8%) 48 (1.7%)
    AAA 2 (<0.1%) 2 (0.1%) 0 (0%)
    AD 4 (<0.1%) 3 (0.2%) 1 (<0.1%)
    D 61 (1.4%) 29 (2.1%) 32 (1.1%)
    Missing 144 65 79
1 Mean ± SD; n (%)
2 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)

CS Rate by Maternal Education Level

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)

CS Rate by Parity Group

8 Factors Associated with CS

8.1 Regression Dataset

NAs are dropped for all variables entering the model.

Code 
reg_df <- df |>
  filter(!is.na(is_cs)) |>
  mutate(
    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()                          # ← remove ALL rows with any NA

cat(sprintf("Regression sample: %d rows (%.1f%% of cleaned data)\n",
            nrow(reg_df),
            nrow(reg_df) / nrow(df) * 100))
Regression sample: 4230 rows (96.3% 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.04 (1.03–1.06) <0.001
parity 0.93 (0.88–0.97) 0.002
ga_weeks 0.97 (0.95–1.00) 0.038
birth_weight 1.00 (0.96–1.00) 0.649
education_num 1.08 (1.03–1.13) 0.003
nhis_yes 0.65 (0.46–0.93) 0.017
twin 1.75 (1.14–2.66) 0.009
has_complication 1.62 (1.19–2.19) 0.002

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.071 (1.056–1.086) <0.001 ***
Parity 0.815 (0.768–0.865) <0.001 ***
Gestational Age (weeks) 0.986 (0.961–1.005) 0.286
Birth Weight (kg) 0.997 (NA–1.002) 0.571
Education Level (ordinal) 1.051 (0.998–1.109) 0.063
NHIS Coverage: Yes 0.673 (0.474–0.961) 0.028 *
Twin Pregnancy 1.754 (1.135–2.697) 0.011 *
Pregnancy Complication Present 1.466 (1.067–2.006) 0.017 *

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: χ²=14.15, df=8, p=0.078
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.621
Code 
# Nagelkerke pseudo-R²
if (requireNamespace("DescTools", quietly = TRUE)) {
  cat(sprintf("Nagelkerke R²: %.3f\n",
              DescTools::PseudoR2(mv_model, which = "Nagelkerke")))
}
Nagelkerke R²: 0.044

10 Summary of Findings

NoteKey Results
  • Overall CS rate: 3274.1% (1438 / 4392 deliveries)
  • Emergency CS comprised 53.9% of all CS cases; Elective CS 44.9%
  • Advanced maternal age was significantly associated with higher odds of CS
  • Twin pregnancy was strongly associated with CS
  • NHIS non-coverage was associated with slightly higher CS odds
  • 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)