Use of data from IPUMS PMA is subject to conditions including that users should cite the data appropriately. Use command `ipums_conditions()` for more details.
#read in spatial datalibrary(sf)
Linking to GEOS 3.13.1, GDAL 3.10.2, PROJ 9.5.1; sf_use_s2() is TRUE
UG_PMA_GPS<-read_sf("C:/Users/Rebecca/Downloads/PMA_UG_GPS_v2_06July2022/UGANDA/PMA_UG_GPS_v2_06July2022.csv")#read in boundary dataUG_boundary<-read_sf("C:/Users/Rebecca/Downloads/geoug/geoug.shp")UG_boundary_4<-read_sf("C:/Users/Rebecca/Downloads/geouggen/geouggen.shp")UG_sub_boundary<-read_sf("C:/Users/Rebecca/Downloads/uganda-adm-boundaries/UGANDA BOUNDARIES SHAPEFILES AS OF 17 08 2018/COUNTIES_2018_UTM_36N.shp")pmaug2022<-pmauglong#%>%#filter(YEAR==2021)UG_PMA_GPS_2022<-UG_PMA_GPS#%>%#filter(PMAYEAR==2021)UG_PMA_GPS_2022<-UG_PMA_GPS%>%rename('EAID'='EA_ID')pmaug2022wgps<-merge(pmaug2022, UG_PMA_GPS_2022, by=c("EAID"))
Step 2. Recoding/ building summary scale, data cleaning
#marital status pmaug2022wgps<-pmaug2022wgps%>%mutate(MARSTAT=as.factor(MARSTAT))%>%mutate(maritalstatus=recode(MARSTAT, '10'="never married" ,'20'="married or living together", '21'="currently married", '22'="currently living with partner", '31'="formerly in union",'32'='widow or widower',.default =NA_character_))tabyl(pmaug2022wgps, maritalstatus)
maritalstatus n percent valid_percent
never married 3506 0.2472845253 0.24735431
currently married 3416 0.2409366624 0.24100466
currently living with partner 4980 0.3512484130 0.35134754
formerly in union 1947 0.1373254338 0.13736419
widow or widower 325 0.0229228382 0.02292931
<NA> 4 0.0002821272 NA
#age-grouppmaug2022wgps<-pmaug2022wgps %>%mutate(agegroup =case_when( AGE >=15& AGE <=19~"15-19", AGE >=20& AGE <=24~"20-24", AGE >=25& AGE <=29~"25-29", AGE >=30& AGE <=34~"30-34", AGE >=35& AGE <=39~"35-39", AGE >=40& AGE <=44~"40-44", AGE >=45& AGE <=49~"45-49" ))pmaug2022wgps$age_z_score <-scale(pmaug2022wgps$AGE)pmaug2022wgps<-pmaug2022wgps%>%filter(BIRTHEVENT<90)pmaug2022wgps<-pmaug2022wgps %>%mutate(AGE=as.numeric(AGE))pmaug2022wgps<-pmaug2022wgps %>%mutate(URBAN=as.factor(URBAN))#contraceptive usepmaug2022wgps<-pmaug2022wgps%>%mutate(MCP=as.factor(MCP))%>%mutate(usingmoderncon=recode(MCP, '1'="yes" ,'0'="no",.default =NA_character_))tabyl(pmaug2022wgps, usingmoderncon)
usingmoderncon n percent valid_percent
no 9304 0.65645947 0.6662847
yes 4660 0.32879419 0.3337153
<NA> 209 0.01474635 NA
educationlevel n percent valid_percent
none 842 0.0594087349 0.05941712
primary/middle school 7902 0.5575389826 0.55761767
secondary/post-primary 4268 0.3011359627 0.30117846
tertiary/ post-secondary 1159 0.0817752064 0.08178675
<NA> 2 0.0001411134 NA
#collapse marital status variablepmaug2022wgps<-pmaug2022wgps%>%mutate(maritalcombined=recode(maritalstatus, 'formerly in union'='not in a union' ,'never married'="not in a union", 'widow or widower'="not in a union", 'currently married'="in a union", 'currently living with partner'="in a union"))#make index: exercise of choicelibrary(janitor)pmaug2022wgps <- pmaug2022wgps %>%mutate(startstop =case_when( STARTKIDDEC ==1| PTRDISCKIDSTOPWILL==1~'1', STARTKIDDEC ==2| PTRDISCKIDSTOPWILL==2~'2', STARTKIDDEC ==3| PTRDISCKIDSTOPWILL==3~'3', STARTKIDDEC ==4| PTRDISCKIDSTOPWILL==4~'4', STARTKIDDEC ==5| PTRDISCKIDSTOPWILL==5~'5' ))tabyl(pmaug2022wgps$startstop)
#convert to a spatial framepmaug2022wgps_sf<-st_as_sf(pmauglongfilter2, coords =c("GPSLONG", "GPSLAT"), crs=4326)
Step 3: Visualize the average summary score by PMA cluster
#map out the scale by EAIDUG<-st_transform(UG_boundary_4, crs=4326)UG <-st_make_valid(UG)joined<-st_join(UG, pmaug2022wgps_sf, join=st_intersects)# Set the global optionoptions(survey.lonely.psu ="adjust")library(dplyr)joined<-joined%>%filter(summary_score_av!='NA')joined$fertautonomyscale <-as.numeric(joined$summary_score_av)library(dplyr)library(srvyr)
Attaching package: 'srvyr'
The following object is masked from 'package:stats':
filter
#convert to a spatial framemapdatareg<-st_as_sf(mapdatareg, coords =c("GPSLONG", "GPSLAT"), crs=4326)#all three yearslibrary(ggplot2)ggplot(data = mapdatareg) +geom_sf(data = UG_boundary_4, fill =NA, color ="lightgray", linewidth =0.2, inherit.aes =FALSE) +geom_sf(data = UG_sub_boundary, fill =NA, color ="lightgray", linewidth =0.2, inherit.aes =FALSE) +geom_sf(data = mapdatareg, aes(color = summary_score_av_mean), size =1.2, alpha =0.8) +scale_color_viridis_c(option ="magma", direction =-1, name ="Pregnancy exercise of choice\nsub-scale mean") +facet_wrap(~ YEAR) +theme_minimal(base_size =12) +theme(legend.title =element_text(size =10),legend.text =element_text(size =9),strip.text =element_text(size =12, face ="bold") )
Step 4: check for spatial clustering
#check for spatial clusteringlibrary(spdep)
Loading required package: spData
To access larger datasets in this package, install the spDataLarge
package with: `install.packages('spDataLarge',
repos='https://nowosad.github.io/drat/', type='source')`
Moran I test under randomisation
data: mapdata_proj$summary_score_av_mean
weights: lw
n reduced by no-neighbour observations
Moran I statistic standard deviate = 3.3426, p-value = 0.0004149
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.225683686 -0.008403361 0.004904301
Moran I test under randomisation
data: mapdata_proj$summary_score_av_mean
weights: lw
n reduced by no-neighbour observations
Moran I statistic standard deviate = 3.3277, p-value = 0.0004378
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.191137691 -0.007246377 0.003554025
Moran I test under randomisation
data: mapdata_proj$summary_score_av_mean
weights: lw
n reduced by no-neighbour observations
Moran I statistic standard deviate = 2.11, p-value = 0.01743
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.118176049 -0.007246377 0.003533431
library(spdep)library(sf)# Transform to UTM Zone 36Nmapdata_proj <-st_transform(mapdatareg, crs =32636)# Use centroids for coordinatescoords_proj <-st_coordinates(st_centroid(mapdata_proj))# Create neighbors within 100 kmnb <-dnearneigh(coords_proj, 0, 50000)
Warning in dnearneigh(coords_proj, 0, 50000): neighbour object has 13
sub-graphs
lw <-nb2listw(nb, style ="W", zero.policy =TRUE)# Check column existence and run Moran's I (univariate example)moran_result <-moran.mc(mapdata_proj$summary_score_av_mean, lw, nsim =499)print(moran_result)
Monte-Carlo simulation of Moran I
data: mapdata_proj$summary_score_av_mean
weights: lw
number of simulations + 1: 500
statistic = 0.28762, observed rank = 500, p-value = 0.002
alternative hypothesis: greater
#bivariate relationship with urban status and summary scoremodelbiv <-lm(summary_score_av ~factor(URBAN), data = individualandcluster)summary(modelbiv)
Call:
lm(formula = summary_score_av ~ factor(URBAN), data = individualandcluster)
Residuals:
Min 1Q Median 3Q Max
-2.97222 -0.30555 0.09765 0.43098 1.09765
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.902353 0.008945 436.237 < 2e-16 ***
factor(URBAN)1 0.069863 0.014962 4.669 3.05e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7904 on 12148 degrees of freedom
Multiple R-squared: 0.001792, Adjusted R-squared: 0.001709
F-statistic: 21.8 on 1 and 12148 DF, p-value: 3.055e-06
#investigate interactions between year and urban statusmodelmulti_interact <-lmer(summary_score_av ~ BIRTHEVENT + educationlevel + agegroup + maritalcombined +factor(URBAN)+factor(YEAR)+factor(MCP) +factor(URBAN) *factor(YEAR) +factor(WEALTHQ) + (1|FQINSTID),data = individualandcluster)modelsummary(modelmulti_interact, stars =TRUE)
(1)
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept)
3.683***
(0.039)
BIRTHEVENT
-0.014**
(0.005)
educationlevelprimary/middle school
-0.008
(0.033)
educationlevelsecondary/post-primary
0.208***
(0.036)
educationleveltertiary/ post-secondary
0.262***
(0.044)
agegroup20-24
-0.017
(0.024)
agegroup25-29
0.006
(0.028)
agegroup30-34
0.049
(0.033)
agegroup35-39
0.001
(0.038)
agegroup40-44
-0.031
(0.042)
agegroup45-49
-0.028
(0.046)
maritalcombinedin a union
-0.053**
(0.018)
factor(URBAN)1
-0.020
(0.027)
factor(YEAR)2021
0.101***
(0.019)
factor(YEAR)2022
0.061**
(0.020)
factor(MCP)1
0.105***
(0.016)
factor(WEALTHQ)2
0.135***
(0.023)
factor(WEALTHQ)3
0.220***
(0.024)
factor(WEALTHQ)4
0.241***
(0.026)
factor(WEALTHQ)5
0.301***
(0.028)
factor(URBAN)1 × factor(YEAR)2021
-0.127***
(0.032)
factor(URBAN)1 × factor(YEAR)2022
-0.099**
(0.033)
SD (Intercept FQINSTID)
0.391
SD (Observations)
0.658
Num.Obs.
12150
R2 Marg.
0.065
R2 Cond.
0.309
AIC
27731.1
BIC
27908.8
ICC
0.3
RMSE
0.58
models3list <-list("Model 1: Base Model"= modelmulti,"Model 2: Year x Urban"= modelmulti_interact)modelsummary( models3list,slope ="b",stars =TRUE,title ="Multivariate relationship between urban/ rural status and pregnancy decision-making",gof_map =NA)
Multivariate relationship between urban/ rural status and pregnancy decision-making
Model 1: Base Model
Model 2: Year x Urban
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept)
3.708***
3.683***
(0.038)
(0.039)
BIRTHEVENT
-0.014**
-0.014**
(0.005)
(0.005)
educationlevelprimary/middle school
-0.007
-0.008
(0.033)
(0.033)
educationlevelsecondary/post-primary
0.210***
0.208***
(0.036)
(0.036)
educationleveltertiary/ post-secondary
0.263***
0.262***
(0.044)
(0.044)
agegroup20-24
-0.017
-0.017
(0.024)
(0.024)
agegroup25-29
0.005
0.006
(0.028)
(0.028)
agegroup30-34
0.047
0.049
(0.033)
(0.033)
agegroup35-39
-0.001
0.001
(0.038)
(0.038)
agegroup40-44
-0.032
-0.031
(0.043)
(0.042)
agegroup45-49
-0.029
-0.028
(0.046)
(0.046)
maritalcombinedin a union
-0.052**
-0.053**
(0.018)
(0.018)
factor(MCP)1
0.105***
0.105***
(0.016)
(0.016)
factor(URBAN)1
-0.096***
-0.020
(0.018)
(0.027)
factor(WEALTHQ)2
0.136***
0.135***
(0.023)
(0.023)
factor(WEALTHQ)3
0.222***
0.220***
(0.024)
(0.024)
factor(WEALTHQ)4
0.242***
0.241***
(0.026)
(0.026)
factor(WEALTHQ)5
0.299***
0.301***
(0.028)
(0.028)
factor(YEAR)2021
0.056***
0.101***
(0.016)
(0.019)
factor(YEAR)2022
0.026
0.061**
(0.016)
(0.020)
factor(URBAN)1 × factor(YEAR)2021
-0.127***
(0.032)
factor(URBAN)1 × factor(YEAR)2022
-0.099**
(0.033)
SD (Intercept FQINSTID)
0.392
0.391
SD (Observations)
0.658
0.658
Urban status is significant and positive at the bivariate level. At the multivariate level it is also significant but it shows a negative relationship with pregnancy decision-making. When year is interacted with urban status we see that the urban effect is concentrated in 2021 and 2022 of the panel.
modelsummary( models,slope ="b",stars =TRUE,title ="Multivariate relationship between urban/ rural status and pregnancy decision-making, stratified by year",gof_map =NA)
Multivariate relationship between urban/ rural status and pregnancy decision-making, stratified by year
2020
2021
2022
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept)
3.638***
3.840***
3.774***
(0.059)
(0.058)
(0.062)
BIRTHEVENT
-0.011
-0.021**
-0.007
(0.009)
(0.008)
(0.008)
educationlevelprimary/middle school
0.025
-0.022
-0.065
(0.052)
(0.050)
(0.055)
educationlevelsecondary/post-primary
0.251***
0.167**
0.163**
(0.059)
(0.056)
(0.060)
educationleveltertiary/ post-secondary
0.330***
0.221***
0.183*
(0.072)
(0.067)
(0.074)
agegroup20-24
-0.009
-0.047
-0.016
(0.042)
(0.038)
(0.040)
agegroup25-29
0.033
-0.062
0.029
(0.048)
(0.044)
(0.045)
agegroup30-34
0.042
0.017
0.067
(0.055)
(0.051)
(0.051)
agegroup35-39
0.052
-0.024
-0.062
(0.064)
(0.057)
(0.059)
agegroup40-44
-0.056
-0.067
-0.026
(0.071)
(0.065)
(0.065)
agegroup45-49
-0.033
-0.010
-0.031
(0.077)
(0.069)
(0.072)
maritalcombinedin a union
-0.068*
-0.016
-0.073*
(0.032)
(0.028)
(0.028)
factor(MCP)1
0.043
0.147***
0.148***
(0.029)
(0.025)
(0.026)
factor(URBAN)1
-0.034
-0.118***
-0.120***
(0.030)
(0.028)
(0.028)
factor(WEALTHQ)2
0.134**
0.091*
0.159***
(0.042)
(0.038)
(0.036)
factor(WEALTHQ)3
0.225***
0.243***
0.189***
(0.042)
(0.038)
(0.037)
factor(WEALTHQ)4
0.269***
0.214***
0.232***
(0.044)
(0.039)
(0.041)
factor(WEALTHQ)5
0.385***
0.220***
0.317***
(0.047)
(0.042)
(0.045)
Urban status is significant in 2021 and 2022 but not 2020
