Multilevel Modelling - Predictors of Client Satisfaction with Family Planning Services in Ethiopia
Michael Endale Mengesha and Henrik Holmberg
Department of Epidemiology and Global Health, Umeå UniversityBackground
This is a supplementary material for the paper titled ‘A multilevel analysis on the predictors of client satisfaction with family planning services in Ethiopia: Evidence from the Ethiopian Service Provision Assessment (ESPA) 2021/22’. It will explore the steps taken during the analysis and reports additional information that might be of interest for the reader.
The data was accessed through the demographic and health survey website.
For the paper two data sets were merged; namely the client data-set (containing client-provider observations and client exit interviews) and the facility data-set (containing facility inventory data). This are named ‘fpclient’ and’fpfacility’ respectively.
We used a left_join with inv_id as unique identifier so the unit of analysis remained as the number of clients.
alldata <- fpclient %>%
left_join(fpfacility, by='inv_id')The variables included in our analysis are described in detail as an appendix in the paper.
Composite measure for client satisfaction
We initially created a composite variable to measure client satisfaction. This composite variable was constructed using survey question number 202 (refer to the final SPA report), which inquiries about common problems clients face at health facilities.
attach(alldata)
overall_client_satisfaction <- data.frame(c502a,c502b,c502c,c502e,
c502f,c502g,c502h,c502i,c502j,c502k,
c502l,c_overall_satisfaction)
overall_client_satisfaction[overall_client_satisfaction == 8] <- NA
overall_client_satisfaction[overall_client_satisfaction == 9] <- NA
overall_client_satisfaction <- na.omit(overall_client_satisfaction)
composite_client_satisfaction <- overall_client_satisfaction %>%
select(c502a, c502b, c502c, c502e, c502f, c502g, c502h, c502i, c502j, c502k, c502l)
pca_composite <- prcomp(composite_client_satisfaction, scale. = TRUE)
pc1 <- pca_composite$x[, 1]
library(factoextra)
fviz_eig(pca_composite, addlabels = TRUE, main = "Variance Explained by Principal Components")
However, once the principal components of the composite measure were identified, a spearman correlation test was conducted to see if the composite variable is actually measuring in line with the overall client satisfaction.
single_variable <- overall_client_satisfaction %>%
pull(c_overall_satisfaction)
if (length(pc1) == length(single_variable)) {
cor_spearman <- cor(pc1, single_variable, method = "spearman")
print(paste("Spearman Rank Correlation: ", cor_spearman))
} else {
stop("Dimensions of PC1 and c_overall_satisfaction do not match!")
}## [1] "Spearman Rank Correlation: 0.269412558848388"Single variable measure for client satisfaction
The very low correlation between the composite measure and the overall client satisfaction indicates the composite variable is measuring a dimension of patient service that is not correlated with client satisfaction. As a result, a decision to use the single variable measure was made.
In our model client satisfaction is measured using a 4-level likert scale that has been dicothomized into more satisfied and less satisfied.
Weights
#Creating facility weight
weight_list <- data.frame(final_df1$f_weight, final_df1$inv_id)
weight_list <- unique(weight_list)
facility_list <- data.frame(final_df1$f_number, final_df1$inv_id)
facility_list<- unique(facility_list)
weight_facility <- weight_list %>%
inner_join(facility_list, by='final_df1.inv_id')
n1 <- nrow(weight_facility)
weight_facility$final_df1.f_weight <- n1 * (weight_facility$final_df1.f_weight /1000000) /
sum(weight_facility$final_df1.f_weight /1000000)
sum(weight_facility$final_df1.f_weight)## [1] 529# Creating client-weight
n <- nrow(final_df1)
final_df1$c_houseweight <- n * (final_df1$weight/1000000) / sum(final_df1$weight/1000000)
sum(final_df1$c_houseweight)## [1] 2071A descriptive table before and after applying weights is presented in the paper.
Multilevel modelling
Facility as second level
nullmodel_adapt <- mixed_model(fixed = c_overall_satisfaction ~1 ,
random = ~ 1 |f_number,
data = final_df1,
family = binomial("logit"),
weights = weight_facility$final_df1.f_weight)
summary(nullmodel_adapt)##
## Call:
## mixed_model(fixed = c_overall_satisfaction ~ 1, random = ~1 |
## f_number, data = final_df1, family = binomial("logit"), weights = weight_facility$final_df1.f_weight)
##
## Data Descriptives:
## Number of Observations: 2071
## Number of Groups: 529
##
## Model:
## family: binomial
## link: logit
##
## Fit statistics:
## log.Lik AIC BIC
## -963.5288 1931.058 1939.6
##
## Random effects covariance matrix:
## StdDev
## (Intercept) 2.0497
##
## Fixed effects:
## Estimate Std.Err z-value p-value
## (Intercept) 0.3539 0.1213 2.9166 0.003539
##
## Integration:
## method: adaptive Gauss-Hermite quadrature rule
## quadrature points: 11
##
## Optimization:
## method: hybrid EM and quasi-Newton
## converged: TRUEtab_model(nullmodel_adapt)| Â | c_overall_satisfaction | ||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 1.42 | 1.12 – 1.81 | 0.004 |
| Random Effects | |||
| σ2 | 3.29 | ||
| Ï„00 f_number | 4.20 | ||
| ICC | 0.56 | ||
| N f_number | 529 | ||
| Observations | 2071 | ||
| Marginal R2 / Conditional R2 | 0.000 / 0.561 | ||
#Median odds ratio
var.area.mnull <- nullmodel_adapt$D
or.mnull <- exp(0.95*sqrt(var.area.mnull))
or.mnull## (Intercept)
## (Intercept) 7.00914
## attr(,"L")
## [1] 2.0497Region as third level?
Conducting a similar test for the variable ‘region’ gives the following result.
Non significant variance attribution value for region renders it’s use as a third level irrelevant.
Level two variance visualized
library(ggplot2)
library(egg)
# ggplot butterfly
re = data.frame(condval = unlist(nullmodel_adapt$post_modes)[,1],
condsd = sqrt(unlist(nullmodel_adapt$post_vars)),
grp = 1:length(nullmodel_adapt$post_vars))
# sort for plot
dataPlot = re[order(re$condval), ]
dataPlot$rank = 1:NROW(dataPlot)
dataPlot$li = dataPlot$condval - dataPlot$condsd*qnorm(0.975)
dataPlot$ui = dataPlot$condval + dataPlot$condsd*qnorm(0.975)
ggplot(data = dataPlot, mapping = aes(x = rank, y = condval,
ymin = li,
ymax = ui))+
geom_ribbon(fill = "steelblue", color = NA)+
geom_line()+
egg::theme_article()
#less dense plot
ggplot(data = dataPlot[as.numeric(as.character(dataPlot$grp)) < 50,],
mapping = aes(x = rank,
y = condval, ymin = li, ymax = ui))+
ylab("Random effects for facility intercept")+
xlab("Facility rank")+
geom_hline(color = "pink", yintercept = 0, linetype = "dashed", alpha
= 0.5)+
geom_errorbar(color = "steelblue", width = 5)+
geom_line(color = "steelblue", alpha = 0.1)+
geom_point(color = "steelblue")+
egg::theme_article()
Client level predictors
Crude estimates
| Variable | Level | OR | Lower.CI | Upper.CI | p.value |
|---|---|---|---|---|---|
| c_fp_history_physical | c_fp_history_physical | 1.0314735 | 0.9558653 | 1.1130622 | 0.4249611 |
| c_previous_contact | c_previous_contactYes | 1.2112652 | 0.8464520 | 1.7333096 | 0.2945159 |
| c_visual_aids | c_visual_aidsYes | 1.1499795 | 0.7200107 | 1.8367127 | 0.5585737 |
| c_confidentiality | c_confidentialityYes | 1.2211428 | 0.7630783 | 1.9541767 | 0.4049390 |
| c_auditory_privacy | c_auditory_privacyYes | 1.6067342 | 1.0898727 | 2.3687121 | 0.0166394 |
| c_visual_privacy | c_visual_privacyYes | 1.6599048 | 1.0631111 | 2.5917178 | 0.0257997 |
| c_facility_closeness | c_facility_closenessYes | 1.3985798 | 0.8394186 | 2.3302145 | 0.1977675 |
| c_prepayment_plan | c_prepayment_planYes | 0.9681164 | 0.6915261 | 1.3553347 | 0.8502779 |
| c_fp_status | c_fp_statusPast User | 1.8384331 | 1.1434912 | 2.9557168 | 0.0119546 |
| c_fp_status | c_fp_statusNever Used | 1.0526124 | 0.6133096 | 1.8065798 | 0.8524022 |
| c_wait_time | c_wait_time | 0.9919341 | 0.9862190 | 0.9976823 | 0.0060126 |
| c_asked_ques | c_asked_quesYes | 0.5763606 | 0.4070327 | 0.8161299 | 0.0019035 |
| c_concerns_disscussed | c_concerns_disscussedYes | 0.6921383 | 0.4877266 | 0.9822213 | 0.0393558 |
| c_partner_attitude | c_partner_attitudeYes | 0.9381833 | 0.4988424 | 1.7644609 | 0.8430454 |
| c_followup | c_followupYes | 1.1020399 | 0.4545849 | 2.6716506 | 0.8297236 |
| c_problem_occur | c_problem_occurYes | 2.6360889 | 1.8987786 | 3.6597025 | 0.0000000 |
| c_sideeffect | c_sideeffectYes | 1.9166866 | 1.3884664 | 2.6458598 | 0.0000764 |
| c_how_to_use | c_how_to_useYes | 1.2376985 | 0.8479216 | 1.8066502 | 0.2691101 |
| c_educationlevel | c_educationlevelPrimary | 0.5224971 | 0.3678555 | 0.7421480 | 0.0002884 |
| c_educationlevel | c_educationlevelHigher | 0.6390332 | 0.4312163 | 0.9470039 | 0.0256594 |
| c_maritalstatus | c_maritalstatusMarried | 0.9940164 | 0.5220273 | 1.8927526 | 0.9854277 |
| c_age | c_age22 TO 30 | 1.3739005 | 0.9764225 | 1.9331823 | 0.0682937 |
| c_age | c_ageMore than 30 | 3.5774911 | 2.2794761 | 5.6146421 | 0.0000000 |
Adjusted estimates
Adjusted estimates along with LLR, median odds ratio and proportional change in variance are reported in the original paper.
Facility level predictors
Crude estimates
| Variable | Level | OR | Lower.CI | Upper.CI | p.value |
|---|---|---|---|---|---|
| f_training | f_trainingYes | 0.5532235 | 0.3330734 | 0.9188853 | 0.0222091 |
| f_catchment_pop | f_catchment_pop25K TO 1.5M | 1.0809172 | 0.6401130 | 1.8252745 | 0.7709833 |
| f_catchment_pop | f_catchment_popMore than 1.5M | 2.8299819 | 1.2142779 | 6.5955227 | 0.0159640 |
| f_no_of_days | f_no_of_daysAbove 20 | 1.1870360 | 0.7214108 | 1.9531928 | 0.4997951 |
| f_quality_structure | f_quality_structureQuality Unit & Committe | 2.4422251 | 1.2224772 | 4.8789975 | 0.0114408 |
| f_quality_structure | f_quality_structureQuality Committe | 3.1102251 | 1.9098011 | 5.0651874 | 0.0000051 |
| f_guideline_national | f_guideline_nationalAvailable | 1.2224549 | 0.7581554 | 1.9710947 | 0.4098922 |
| f_routine_action | f_routine_actionYes | 0.9107589 | 0.5672359 | 1.4623224 | 0.6988041 |
| f_location | f_locationRural | 0.4770849 | 0.2923226 | 0.7786263 | 0.0030641 |
| f_ownership | f_ownershipPrivate | 2.3901438 | 1.0490418 | 5.4457196 | 0.0380835 |
| f_no_trained_provider | f_no_trained_providerYes | 0.6500749 | 0.3908701 | 1.0811708 | 0.0970544 |
| f_methods_offered | f_methods_offeredMore than/equal 8 | 0.9599982 | 0.5915292 | 1.5579897 | 0.8687513 |
| fp_equipment | fp_equipment | 1.1034273 | 0.9773649 | 1.2457495 | 0.1118118 |
| f_supervisory_visit | f_supervisory_visitWithin 6 month | 1.4155655 | 0.5973812 | 3.3543499 | 0.4297978 |
| f_supervisory_visit | f_supervisory_visitMore than 6 month | 1.0326358 | 0.5342701 | 1.9958756 | 0.9239020 |
| f_routine_fees | f_routine_feesFixed | 0.3650711 | 0.1649555 | 0.8079568 | 0.0129142 |
| f_routine_fees | f_routine_feesSeparate | 2.4775930 | 1.5269827 | 4.0199979 | 0.0002386 |
| f_monthly_meetings | f_monthly_meetingsYes | 0.7123443 | 0.4053152 | 1.2519501 | 0.2384079 |
| f_3month_comm_meetings | f_3month_comm_meetingsYes | 0.8519934 | 0.5274847 | 1.3761399 | 0.5126016 |
| f_waiting_area | f_waiting_areaYes | 4.9687153 | 2.6125993 | 9.4496432 | 0.0000010 |
| f_standard_precaution | f_standard_precautionYes | 1.5990471 | 0.9738851 | 2.6255166 | 0.0635384 |
| f_performance_check | f_performance_checkClinical audit | 1.3936506 | 0.7250508 | 2.6787942 | 0.3194362 |
| f_performance_check | f_performance_checkRegular audit | 1.1589273 | 0.6638520 | 2.0232108 | 0.6038762 |
| f_performance_check | f_performance_checkFeedback | 3.1588163 | 1.4410269 | 6.9243122 | 0.0040735 |
| f_dhis2 | f_dhis2Yes | 1.7049208 | 1.0626578 | 2.7353630 | 0.0269691 |
| f_risk_assessment | f_risk_assessmentYes | 1.0449422 | 0.6019792 | 1.8138571 | 0.8758448 |
Adjusted estimates
Adjusted estimates along with LLR, median odds ratio and proportional change in variance are reported in the original paper.
Final model
A combination of client and facility level variables that were significant in the previous models are reported for adjusted estimates along with LLR, median odds ratio and proportional change in variance in the original paper.
Thank you for taking interest in our work.