Healthcare accessibility among women remains a major public health and development concern in Bangladesh, where socioeconomic, geographic, educational, and informational inequalities can shape access to care. Using data from the Bangladesh Demographic and Health Survey (BDHS) 2017–18, this analysis examines healthcare access problems among women by wealth index, rural-urban residence, education level, marital status, and media exposure.

This project connects women’s health access to poverty and inequality research by exploring whether disadvantaged groups face greater barriers to healthcare. The analysis uses R to clean and prepare DHS microdata, recode key variables, apply survey weights, and generate descriptive and statistical outputs. By examining healthcare access through socioeconomic and demographic factors, the project provides policy-relevant evidence for understanding health equity, rural service gaps, and poverty-sensitive public health planning in Bangladesh.

The original DHS microdata are not included due to DHS data-use restrictions. Researchers can request access from The DHS Program and place individual.DTA in the local data/raw folder before running the analysis.

#all necessary libraries for data cleaning, analysis, and visualization
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.1     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.3     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(readxl)
library(janitor)

## 
## Attaching package: 'janitor'
## 
## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test

library(haven)
library(labelled)
library(survey)

## Loading required package: grid
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## Loading required package: survival
## 
## Attaching package: 'survey'
## 
## The following object is masked from 'package:graphics':
## 
##     dotchart

library(srvyr)

## 
## Attaching package: 'srvyr'
## 
## The following object is masked from 'package:stats':
## 
##     filter

library(skimr)
library(gtsummary)
library(flextable)

## 
## Attaching package: 'flextable'
## 
## The following object is masked from 'package:gtsummary':
## 
##     continuous_summary
## 
## The following object is masked from 'package:purrr':
## 
##     compose

library(ggplot2)
library(gt)
library(dplyr)

Cleaning and preparing data for analysis from main dataset, as per the research question and variables of interest. Used dhs_codebook.pdf for variable selection and coding. It can be found in the same link as the dataset().

main_dataset <- read_dta("D:/DHS_research/dataraw/bangladesh/individual.DTA")
analysis_data <- main_dataset %>%
select(v001,v002,v005,v021,v022, v012, v024, v025, v106, v130, v190, v502,v394,v467b,v467c,v467d,v467f,v481,v157,v158,v159)
#demographics,healthcare, media exposure
analysis_data <- analysis_data %>% #turning into factors for analysis
mutate(across(where(is.labelled),to_factor))
#adding weight for survey desing(as per dhs rules)
analysis_data <-analysis_data %>%
mutate(weight=v005/1000000)
#Checking the na in the data
table(analysis_data$v467b, useNA = "always")

## 
##        no problem       big problem not a big problem              <NA> 
##                 0              2173             17952                 2

#access problem and binary indicator for women's healthcare access barriers
analysis_data <- analysis_data %>%
mutate(access_problem = if_else(v467b == "big problem"|v467c =="big problem"|v467d=="big problem"|v467f == "big problem","Has Access Problem","No Access Problem")) #
analysis_data <- analysis_data %>% 
mutate(access_bin = if_else(access_problem == "Has Access Problem",1,0)) 
#survey design for analysis
dhs_design <- svydesign(ids = ~v021,strata = ~v022,weights = ~weight,data = analysis_data,nest = TRUE)
dhs_design

## Stratified 1 - level Cluster Sampling design (with replacement)
## With (672) clusters.
## svydesign(ids = ~v021, strata = ~v022, weights = ~weight, data = analysis_data, 
##     nest = TRUE)

women’s healthcare access problems (%) by wealth group

#grouped wealth variable
analysis_data <- analysis_data %>%
mutate(wealth_group = case_when(  v190 %in% c("poorest", "poorer") ~ "Poor",  v190 == "middle" ~ "Middle Income",  v190 %in% c("richer", "richest") ~ "Rich"))
#order setting for better visualization
analysis_data$wealth_group <- factor(
analysis_data$wealth_group,
levels = c("Poor", "Middle Income", "Rich"))
# recreate survey design after adding wealth_group
dhs_design <- svydesign(ids = ~v021,strata = ~v022,weights = ~weight,data = analysis_data,nest = TRUE)
# weighted results by wealth group
wealth_group_results <- svyby(~access_bin,~wealth_group,dhs_design,svymean,na.rm = TRUE)
# convert to percentage
wealth_group_results <- wealth_group_results %>%
  mutate(  access_bin = round(access_bin * 100, 2),  se = round(se * 100, 2))
wealth_group_results

##                wealth_group access_bin   se
## Poor                   Poor      77.74 0.88
## Middle Income Middle Income      68.04 1.09
## Rich                   Rich      56.38 0.98

healthcare access problems by wealth group: figure

wealth_group_results$wealth_group <- factor(wealth_group_results$wealth_group,
levels = c("Poor", "Middle Income", "Rich"))
ggplot(wealth_group_results,aes(x = wealth_group, y = access_bin)) + geom_col(width = 0.7)+
geom_text(aes(label = access_bin),vjust = -0.5,size = 6) +
labs(title= "Healthcare Access Problems by Wealth Group",subtitle = "Bangladesh DHS Women Dataset",x = "Wealth Group",y = "Women With Healthcare Access Problems (%)",
     caption= "Weighted DHS estimates")+ ylim(0, 100)+theme_minimal(base_size = 14)

women’s access to healthcare access(%) by education

education_results <- svyby(~access_bin,~v106,dhs_design,svymean,na.rm = TRUE)
education_results <- education_results %>%
  mutate(access_bin = round(access_bin * 100, 2),se = round(se * 100, 2))
education_results

##                      v106 access_bin   se
## no education no education      77.40 1.01
## primary           primary      73.27 0.85
## secondary       secondary      64.14 0.89
## higher             higher      45.97 1.35

Healthcare access by education level: figure

ggplot(education_results,aes(x = v106, y = access_bin)) +
geom_col(width = 0.7) +geom_text(aes(label = access_bin),vjust = -0.5,size = 8) +
labs(title = "Healthcare Access Problems by Education Level",subtitle = "Bangladesh DHS Women Dataset",x = "Education Level",y = "Women With Healthcare Access Problems (%)",caption = "Weighted DHS estimates") +ylim(0, 100) +theme_minimal(base_size = 14)

women’s access to healthcare access(%) by urban vs rural

residence_results <- svyby(
~access_bin,~v025,dhs_design,svymean,na.rm = TRUE)
residence_results <- residence_results %>%
mutate(
access_bin = round(access_bin * 100, 2),
se = round(se * 100, 2))
residence_results

##        v025 access_bin   se
## urban urban      57.92 1.28
## rural rural      70.47 0.86

urban vs rural healthcare access figure

ggplot(residence_results,aes(x = v025, y = access_bin)) +
geom_col(width = 0.6) +geom_text(aes(label = access_bin),vjust = -0.5,size = 8) +
labs(title = "Healthcare Access Problems by Residence",subtitle = "Bangladesh DHS Women Dataset",x = "Residence",y = "Women With Healthcare Access Problems (%)",
caption ="Weighted DHS estimates")+ ylim(0, 100) + theme_minimal(base_size = 14)

women’s access to healthcare access(%) by marital status

analysis_data <- analysis_data %>%
mutate(marital_group = case_when(v502 == "currently in union/living with a man" ~ "Currently Married",v502 == "formerly in union/living with a man" ~ "Formerly Married",
TRUE ~ NA_character_))
#survey design to include marital status
dhs_design <- svydesign(ids = ~v021,strata = ~v022,weights = ~weight,data = analysis_data,nest = TRUE)
marital_results <- svyby(~access_bin,~marital_group,dhs_design,svymean,na.rm = TRUE)

marital_results

##                       marital_group access_bin          se
## Currently Married Currently Married  0.6663951 0.007278476
## Formerly Married   Formerly Married  0.7127945 0.016542411

marital_results <- marital_results %>%
mutate(access_bin = round(access_bin * 100, 2),se = round(se * 100, 2))
marital_results

##                       marital_group access_bin   se
## Currently Married Currently Married      66.64 0.73
## Formerly Married   Formerly Married      71.28 1.65

women’s access to healthcare by media exposure

#coding numeric for statistical analysis
analysis_data <- analysis_data %>%
mutate(media_exposure = if_else(v157 != "not at all" |v158 != "not at all" |v159 != "not at all",1,0))
table(analysis_data$media_exposure)

## 
##     0     1 
##  6981 13146

#labeling: convinient for interpretation
analysis_data <- analysis_data %>% #c
mutate(media_exposure = case_when(media_exposure == 1 ~ "Has Media Exposure",media_exposure == 0 ~ "No Media Exposure"))
table(analysis_data$media_exposure)

## 
## Has Media Exposure  No Media Exposure 
##              13146               6981

dhs_design <- svydesign(ids = ~v021,strata = ~v022, weights = ~weight,data = analysis_data, nest = TRUE)

media_results <- svyby(~access_bin,~media_exposure,dhs_design,svymean,na.rm = TRUE)
media_results

##                        media_exposure access_bin          se
## Has Media Exposure Has Media Exposure  0.6228659 0.008076747
## No Media Exposure   No Media Exposure  0.7585671 0.009204702

media_results <- media_results %>%
mutate(access_bin = round(access_bin * 100, 2),se = round(se * 100, 2))

media_results

##                        media_exposure access_bin   se
## Has Media Exposure Has Media Exposure      62.29 0.81
## No Media Exposure   No Media Exposure      75.86 0.92

Media exposure figure

ggplot(media_results,aes(x = media_exposure, y = access_bin)) + geom_col(width = 0.6) +
geom_text(aes(label = access_bin),vjust = -0.5,size = 10)+
labs(title = "Healthcare Access Problems by Media Exposure",subtitle = "Bangladesh DHS Women Dataset",
x = "Media Exposure",y = "Women With Healthcare Access Problems (%)",caption = "Weighted DHS estimates")+ylim(0, 100)+
theme_minimal(base_size = 14)

All results together for better understanding.

combined_results <- bind_rows(wealth_group_results %>%  mutate(variable = "Wealth", category = as.character(wealth_group)) %>% 
                                select(variable, category, access_bin, se),
education_results %>% 
  mutate(variable = "Education", category = as.character(v106)) %>% 
  select(variable, category, access_bin, se),
residence_results %>% 
  mutate(variable = "Residence", category = as.character(v025)) %>% 
  select(variable, category, access_bin, se),
marital_results %>% 
  mutate(variable = "Marital Status", category = as.character(marital_group)) %>% 
  select(variable, category, access_bin, se),
media_results %>% 
  mutate(variable = "Media Exposure", category = as.character(media_exposure)) %>% 
  select(variable, category, access_bin, se))
combined_results

##                          variable           category access_bin   se
## Poor                       Wealth               Poor      77.74 0.88
## Middle Income              Wealth      Middle Income      68.04 1.09
## Rich                       Wealth               Rich      56.38 0.98
## no education            Education       no education      77.40 1.01
## primary                 Education            primary      73.27 0.85
## secondary               Education          secondary      64.14 0.89
## higher                  Education             higher      45.97 1.35
## urban                   Residence              urban      57.92 1.28
## rural                   Residence              rural      70.47 0.86
## Currently Married  Marital Status  Currently Married      66.64 0.73
## Formerly Married   Marital Status   Formerly Married      71.28 1.65
## Has Media Exposure Media Exposure Has Media Exposure      62.29 0.81
## No Media Exposure  Media Exposure  No Media Exposure      75.86 0.92

Easily understandable table for all

combined_results %>%
  arrange(desc(access_bin)) %>%
  rename(  `Determinant` = variable,  `Category` = category,  `Healthcare Access Problems (%)` = access_bin,  `Standard Error` = se) %>%
  gt(groupname_col = "Determinant") %>%fmt_number(columns = c(`Healthcare Access Problems (%)`,    `Standard Error`  ),  decimals = 2) %>%
  tab_header(  title = "Women’s Healthcare Access Problems in Bangladesh",  subtitle = "Ranked weighted estimates from Bangladesh DHS") %>%
  cols_label(  Category = "Category",  `Healthcare Access Problems (%)` = "Access Problem (%)",  `Standard Error` = "SE") %>%
  tab_source_note(  source_note = "Note: Estimates account for DHS sampling weights, clustering, and stratification.")

Category	Access Problem (%)	SE
Women’s Healthcare Access Problems in Bangladesh
Ranked weighted estimates from Bangladesh DHS
Wealth
Poor	77.74	0.88
Middle Income	68.04	1.09
Rich	56.38	0.98
Education
no education	77.40	1.01
primary	73.27	0.85
secondary	64.14	0.89
higher	45.97	1.35
Media Exposure
No Media Exposure	75.86	0.92
Has Media Exposure	62.29	0.81
Marital Status
Formerly Married	71.28	1.65
Currently Married	66.64	0.73
Residence
rural	70.47	0.86
urban	57.92	1.28
Note: Estimates account for DHS sampling weights, clustering, and stratification.

Multivariable logistic regression to identify independent predictors of healthcare access problems among

health_model <- svyglm(
access_bin ~wealth_group + v106 +v025 +marital_group + media_exposure + v012 +v130,design = dhs_design,family = quasibinomial())
summary(health_model)

## 
## Call:
## svyglm(formula = access_bin ~ wealth_group + v106 + v025 + marital_group + 
##     media_exposure + v012 + v130, design = dhs_design, family = quasibinomial())
## 
## Survey design:
## svydesign(ids = ~v021, strata = ~v022, weights = ~weight, data = analysis_data, 
##     nest = TRUE)
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      1.316227   0.124011  10.614  < 2e-16 ***
## wealth_groupMiddle Income       -0.307761   0.060641  -5.075 5.09e-07 ***
## wealth_groupRich                -0.567137   0.061890  -9.164  < 2e-16 ***
## v106primary                     -0.175082   0.060309  -2.903 0.003823 ** 
## v106secondary                   -0.495325   0.066469  -7.452 3.00e-13 ***
## v106higher                      -1.063991   0.081217 -13.101  < 2e-16 ***
## v025rural                        0.213466   0.068193   3.130 0.001826 ** 
## marital_groupFormerly Married    0.085588   0.081682   1.048 0.295121    
## media_exposureNo Media Exposure  0.194731   0.052857   3.684 0.000249 ***
## v012                            -0.004112   0.002237  -1.838 0.066468 .  
## v130hinduism                     0.211106   0.086555   2.439 0.015000 *  
## v130buddhism                     0.670836   0.385684   1.739 0.082457 .  
## v130christianity                 0.827997   0.541169   1.530 0.126509    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasibinomial family taken to be 1.000717)
## 
## Number of Fisher Scoring iterations: 4

Interpretation of the results:

1. Women from poor households have significantly higher odds of facing healthcare access barriers compared to women from middle-income households, while women from rich households have significantly lower odds of healthcare access problems.
2. Women with higher education have substantially lower odds of healthcare access problems.
3. Rural women have higher odds of healthcare access problems than urban women.
4. Women without media exposure have higher healthcare access barriers.
5. After controlling for other socioeconomic and demographic factors, marital status does not independently predict healthcare access problems.

Additional visualization for extra points:

ranking the vulnerbility

combined_ranked <- combined_results %>%
arrange(desc(access_bin)) %>%
mutate(
category_label = paste(variable, category, sep = ": "),
category_label = reorder(category_label, access_bin))

ggplot(combined_ranked, aes(x = category_label, y = access_bin)) +
geom_col(width = 0.7) +
geom_text(aes(label = paste0(access_bin, "%")), hjust = -0.1, size = 3.5) +coord_flip() +
labs(
  title = "Most Vulnerable Groups for Healthcare Access Problems",
subtitle = "Ranked weighted estimates from Bangladesh DHS women dataset", x = "",y = "Women With Healthcare Access Problems (%)",
caption = "Weighted DHS estimates") + ylim(0, 100) + theme_minimal(base_size = 13)

Overall, the findings suggest that poverty, lower educational attainment, rural residence, and limited media exposure remain important determinants of healthcare accessibility among women in Bangladesh.

RA_MSU_STAT-Social and demographic inequalities in women’s healthcare accessibility: Bangladesh 2017-18

Shahan Ahmed

2/18/2023

Multivariable logistic regression to identify independent predictors of healthcare access problems among

Interpretation of the results:

Additional visualization for extra points:

Overall, the findings suggest that poverty, lower educational attainment, rural residence, and limited media exposure remain important determinants of healthcare accessibility among women in Bangladesh.