Becca Ch. 3 analysis 9.22.26

Author

R.Luttinen

Ch. 3

1. Load and perform data cleaning/ recoding

#load data
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.2.1     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   4.0.3     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ 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(janitor)


Attaching package: 'janitor'

The following objects are masked from 'package:stats':

    chisq.test, fisher.test

library(ipumsr)

options(survey.lonely.psu = "adjust")

#person-period dataset

pmauglong <- read_ipums_micro(
  ddi = "C:/Users/Rebecca/Downloads/pma_00032.xml",
  data = "C:/Users/Rebecca/Downloads/pma_00032.dat.gz")

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 data

library(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 data

UG_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"))

#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-group

pmaug2022wgps<-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 use

pmaug2022wgps<-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

#education

pmaug2022wgps<-pmaug2022wgps%>%
  mutate(EDUCATTGEN=as.factor(EDUCATTGEN))%>%
  mutate(educationlevel= recode(EDUCATTGEN, '1' = "none" ,'2' = "primary/middle school", '3'= "secondary/post-primary", '4'= 'tertiary/ post-secondary',.default = NA_character_))

tabyl(pmaug2022wgps, educationlevel)

           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 variable

pmaug2022wgps<-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 choice

library(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)

 pmaug2022wgps$startstop    n    percent valid_percent
                       1  598 0.04219290    0.04291045
                       2 1951 0.13765611    0.13999713
                       3  559 0.03944119    0.04011194
                       4 7395 0.52176674    0.53064007
                       5 3433 0.24222112    0.24634041
                    <NA>  237 0.01672194            NA

pmaug2022wgps <- pmaug2022wgps %>%
  mutate(havekid = case_when(
    STARTKIDDECWILL == 1 | NEXTKIDDEC== 1 ~ '1',
    STARTKIDDECWILL == 2 | NEXTKIDDEC== 2 ~ '2',
    STARTKIDDECWILL == 3 | NEXTKIDDEC== 3 ~ '3',
    STARTKIDDECWILL == 4 | NEXTKIDDEC== 4~ '4',
    STARTKIDDECWILL == 5 | NEXTKIDDEC== 5 ~ '5'
  ))


tabyl(pmaug2022wgps$havekid)

 pmaug2022wgps$havekid    n    percent valid_percent
                     1  343 0.02420095    0.02441976
                     2 1267 0.08939533    0.09020362
                     3  446 0.03146828    0.03175281
                     4 7311 0.51583998    0.52050406
                     5 4679 0.33013476    0.33311975
                  <NA>  127 0.00896070            NA

pmaug2022wgps <- pmaug2022wgps %>%
  mutate(negotiate = case_when(
    PTRDISCKIDSTOP == 1 | PTRDISCKIDSTART== 1 ~ '1',
    PTRDISCKIDSTOP == 2 | PTRDISCKIDSTART== 2 ~ '2',
    PTRDISCKIDSTOP == 3 | PTRDISCKIDSTART== 3 ~ '3',
    PTRDISCKIDSTOP == 4 | PTRDISCKIDSTART== 4~ '4',
    PTRDISCKIDSTOP == 5 | PTRDISCKIDSTART== 5 ~ '5'
  ))

tabyl(pmaug2022wgps$negotiate)

 pmaug2022wgps$negotiate    n    percent valid_percent
                       1  575 0.04057010    0.04133420
                       2 1341 0.09461652    0.09639853
                       3  591 0.04169901    0.04248436
                       4 7393 0.52162563    0.53144993
                       5 4011 0.28300289    0.28833297
                    <NA>  262 0.01848585            NA

#make index existence of choice

pmaug2022wgps <- pmaug2022wgps %>%
  mutate(restbtw = case_when(
    PREGREST == 1 ~ '1',
    PREGREST == 2 ~ '2',
    PREGREST == 3 ~ '3',
    PREGREST == 4 ~ '4',
    PREGREST == 5~ '5'
  ))

tabyl(pmaug2022wgps$restbtw)

 pmaug2022wgps$restbtw    n     percent valid_percent
                     1  127 0.008960700   0.009021809
                     2  221 0.015593029   0.015699368
                     3  150 0.010583504   0.010655679
                     4 6855 0.483666126   0.486964552
                     5 6724 0.474423199   0.477658592
                  <NA>   96 0.006773442            NA

pmaug2022wgps <- pmaug2022wgps %>%
  mutate(edufirst = case_when(
    EDUCPREKID == 1 | EDUCPREKIDWILL== 1 ~ '1',
    EDUCPREKID == 2 | EDUCPREKIDWILL== 2 ~ '2',
    EDUCPREKID == 3 | EDUCPREKIDWILL== 3 ~ '3',
    EDUCPREKID == 4 | EDUCPREKIDWILL== 4~ '4',
    EDUCPREKID == 5 | EDUCPREKIDWILL== 5 ~ '5'
  ))

tabyl(pmaug2022wgps$edufirst)

 pmaug2022wgps$edufirst    n     percent valid_percent
                      1  202 0.014252452    0.03886858
                      2  396 0.027940450    0.07619781
                      3  122 0.008607916    0.02347508
                      4 1684 0.118817470    0.32403310
                      5 2793 0.197064842    0.53742544
                   <NA> 8976 0.633316870            NA

#don't use because it isn't validated



#drop na's

pmauglongfilter2 <- pmaug2022wgps %>%
  filter(!is.na(startstop),
         !is.na(havekid),
         !is.na(negotiate))



pmauglongfilter2$startstopnum<-as.numeric(pmauglongfilter2$startstop)
pmauglongfilter2$havekidnum<-as.numeric(pmauglongfilter2$havekid)
pmauglongfilter2$negotiatenum<-as.numeric(pmauglongfilter2$negotiate)

#create summary variables: sum

pmauglongfilter2 <- pmauglongfilter2 %>%
  mutate(summary_score = (startstopnum + havekidnum + negotiatenum) )


tabyl(pmauglongfilter2$summary_score)

 pmauglongfilter2$summary_score    n     percent
                              3   67 0.004862472
                              4   68 0.004935046
                              5   98 0.007112272
                              6  342 0.024820379
                              7  273 0.019812759
                              8  632 0.045866899
                              9  514 0.037303142
                             10 1357 0.098483199
                             11  959 0.069598665
                             12 4526 0.328470861
                             13 1734 0.125843675
                             14 1330 0.096523695
                             15 1879 0.136366935

hist(pmauglongfilter2$summary_score)

#average

pmauglongfilter2 <- pmauglongfilter2 %>%
  mutate(summary_score_av = (startstopnum + havekidnum + negotiatenum)/3 )

tabyl(pmauglongfilter2$summary_score_av)

 pmauglongfilter2$summary_score_av    n     percent
                          1.000000   67 0.004862472
                          1.333333   68 0.004935046
                          1.666667   98 0.007112272
                          2.000000  342 0.024820379
                          2.333333  273 0.019812759
                          2.666667  632 0.045866899
                          3.000000  514 0.037303142
                          3.333333 1357 0.098483199
                          3.666667  959 0.069598665
                          4.000000 4526 0.328470861
                          4.333333 1734 0.125843675
                          4.666667 1330 0.096523695
                          5.000000 1879 0.136366935

hist(pmauglongfilter2$summary_score_av)

#convert to a spatial frame

pmaug2022wgps_sf<- st_as_sf(pmauglongfilter2, coords = c("GPSLONG", "GPSLAT"), crs=4326)

2. Descriptive table stratified by urban/ rural status

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(haven)
library(srvyr)


Attaching package: 'srvyr'

The following object is masked from 'package:stats':

    filter

desc_design <- pmaug2022wgps_sf %>%
  st_drop_geometry() %>%
  filter(!is.na(FQWEIGHT)) %>%
  mutate(
    residence        = factor(URBAN, levels = c(0, 1), labels = c("Rural", "Urban")),
    summary_score_av = as.numeric(summary_score_av),
    BIRTHEVENT       = if_else(BIRTHEVENT >= 90, NA_real_, as.numeric(BIRTHEVENT)),
    educationlevel   = as_factor(educationlevel),
    agegroup         = as_factor(agegroup),
    maritalcombined  = as_factor(maritalcombined),
    MCP              = na_if(as.numeric(as.character(MCP)), 98),
    WEALTHQ          = droplevels(as_factor(replace(WEALTHQ, !(as.numeric(WEALTHQ) %in% 1:5), NA))),
    YEAR             = factor(YEAR)
  ) %>%
  as_survey_design(ids = EAID, strata = STRATA, weights = FQWEIGHT, nest = TRUE) %>%
  filter(!is.na(summary_score_av), !is.na(residence),
         !is.na(MCP), !is.na(educationlevel), !is.na(maritalcombined))

#should equal the regression sample (12,150)
nrow(desc_design)

[1] 12150

table1 <- desc_design %>%
  tbl_svysummary(
    by = residence,
    include = c(summary_score_av, BIRTHEVENT, educationlevel, agegroup,
                maritalcombined, MCP, WEALTHQ, YEAR),
    type = list(MCP ~ "dichotomous"),
    value = list(MCP ~ 1),
    statistic = list(
      all_continuous()  ~ "{mean} ({sd})",
      all_categorical() ~ "{p}%"
    ),
    label = list(
      summary_score_av ~ "Pregnancy exercise of choice sub-scale",
      BIRTHEVENT       ~ "Number of births",
      educationlevel   ~ "Educational attainment",
      agegroup         ~ "Age group",
      maritalcombined  ~ "Marital status",
      MCP              ~ "Current modern contraceptive use",
      WEALTHQ          ~ "Wealth quintile",
      YEAR             ~ "Survey year"
    ),
    missing = "no"
  ) %>%
  add_overall() %>%
  add_p() %>%
  modify_header(all_stat_cols() ~ "**{level}**  \nN = {n_unweighted}") %>%
  bold_labels()

table1

Characteristic	Overall N = 12150¹	Rural N = 7807¹	Urban N = 4343¹	p-value²
Pregnancy exercise of choice sub-scale	3.95 (0.78)	3.93 (0.80)	4.03 (0.72)	0.14
Number of births	2.91 (2.83)	3.22 (2.99)	2.08 (2.16)	<0.001
Educational attainment				<0.001
none	5.4%	6.4%	2.7%
primary/middle school	55%	63%	36%
secondary/post-primary	31%	26%	43%
tertiary/ post-secondary	8.5%	4.8%	19%
Age group				<0.001
25-29	17%	16%	21%
35-39	12%	12%	11%
15-19	22%	24%	17%
40-44	8.0%	8.8%	6.0%
45-49	6.8%	7.2%	5.7%
20-24	20%	19%	23%
30-34	14%	14%	15%
Marital status				<0.001
not in a union	39%	37%	46%
in a union	61%	63%	54%
Current modern contraceptive use	33%	32%	37%	<0.001
Wealth quintile				<0.001
Lowest quintile	19%	24%	4.4%
Lower quintile	19%	24%	6.4%
Middle quintile	20%	23%	11%
Higher quintile	20%	18%	26%
Highest quintile	23%	12%	52%
Survey year				0.073
2020	31%	31%	33%
2021	35%	37%	28%
2022	34%	32%	38%
¹ Mean (SD); %
² Design-based KruskalWallis test; Pearson’s X^2: Rao & Scott adjustment

3. Graphics of distribution of outcome by urban status/ wealth

library(srvyr)
library(ggplot2)
library(dplyr)

options(survey.lonely.psu = "adjust")

plotdat <- pmaug2022wgps_sf %>%
  filter(!WEALTHQ %in% c(98, 99)) %>%
  mutate(
    WEALTHQ    = factor(WEALTHQ, levels = 1:5,
                        labels = c("Lowest", "Second", "Middle", "Fourth", "Highest")),
    urbanrural = factor(URBAN, levels = c(0, 1), labels = c("Rural", "Urban"))
  ) %>%
  as_survey_design(ids = EAID, strata = STRATA, weights = FQWEIGHT, nest = TRUE) %>%
  group_by(WEALTHQ, urbanrural) %>%
  summarise(mean_score = survey_mean(summary_score_av, vartype = "ci"), .groups = "drop")

fig <- ggplot(plotdat, aes(x = WEALTHQ, y = mean_score, fill = urbanrural)) +
  geom_col(position = position_dodge(0.8), width = 0.7) +
  geom_errorbar(aes(ymin = mean_score_low, ymax = mean_score_upp),
                position = position_dodge(0.8), width = 0.15, linewidth = 0.4) +
  scale_fill_manual(values = c("Rural" = "#4C6A92", "Urban" = "#C4703F")) +
  coord_cartesian(ylim = c(3, 4.5)) +
  labs(
    x = "Household wealth quintile",
    y = "Mean pregnancy decision-making score",
    fill = "Residence",
    title = "Pregnancy decision-making score by wealth quintile and residence"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    panel.grid.major.x = element_blank(),
    plot.caption = element_text(hjust = 0, color = "grey30"),
    legend.position = "top"
  )

fig

#ggsave("figure_wealth_urban.png", fig, width = 8, height = 5, dpi = 300, bg = "white")

4. Create a new PCA-derived measure separately for those in urban and rural areas; show biplots

#floor/wall/roof materials

recode_material <- function(x) {
  cat <- x %/% 100
  case_when(
    cat == 1 ~ "natural",
    cat == 2 ~ "rudimentary",
    cat == 3 ~ "finished",
    cat == 4 ~ "other",
    cat >= 9 ~ NA_character_,
    TRUE ~ NA_character_
  )
}

pmaug2022wgps <- pmaug2022wgps %>%
  mutate(
    floor_cat = recode_material(FLOOR),
    wall_cat  = recode_material(WALLS),
    roof_cat  = recode_material(ROOF)
  )

tabyl(pmaug2022wgps, floor_cat)

   floor_cat    n      percent valid_percent
    finished 6317 0.4457066253  0.4460528174
     natural 7830 0.5524589007  0.5528880102
       other   11 0.0007761236  0.0007767265
 rudimentary    4 0.0002822268  0.0002824460
        <NA>   11 0.0007761236            NA

tabyl(pmaug2022wgps, wall_cat)

    wall_cat    n      percent valid_percent
    finished 7553 0.5329146970    0.53329097
     natural 3699 0.2609892048    0.26117348
       other  144 0.0101601637    0.01016734
 rudimentary 2767 0.1952303676    0.19536821
        <NA>   10 0.0007055669            NA

tabyl(pmaug2022wgps, roof_cat)

    roof_cat     n      percent valid_percent
    finished   100 0.0070556692   0.007060651
     natural  3181 0.2244408382   0.224599308
       other    17 0.0011994638   0.001200311
 rudimentary 10865 0.7665984619   0.767139730
        <NA>    10 0.0007055669            NA

pmaug2022wgps <- pmaug2022wgps %>%
  mutate(
    floor_cat_natural     = as.integer(floor_cat == "natural"),
    floor_cat_rudimentary = as.integer(floor_cat == "rudimentary"),
    floor_cat_finished    = as.integer(floor_cat == "finished"),
    floor_cat_other       = as.integer(floor_cat == "other"),

    wall_cat_natural      = as.integer(wall_cat == "natural"),
    wall_cat_rudimentary  = as.integer(wall_cat == "rudimentary"),
    wall_cat_finished     = as.integer(wall_cat == "finished"),
    wall_cat_other        = as.integer(wall_cat == "other"),

    roof_cat_natural      = as.integer(roof_cat == "natural"),
    roof_cat_rudimentary  = as.integer(roof_cat == "rudimentary"),
    roof_cat_finished     = as.integer(roof_cat == "finished"),
    roof_cat_other        = as.integer(roof_cat == "other")
  )

#toilet type

pmaug2022wgps <- pmaug2022wgps %>%
  mutate(
    toilet_cat = case_when(
      TOILETTYPE %in% c(1, 2, 3)        ~ "flush",
      TOILETTYPE %in% c(6, 7, 9)        ~ "improved_pit",
      TOILETTYPE %in% c(8, 11)          ~ "unimproved_pit",
      TOILETTYPE %in% c(20, 22)         ~ "none",
      TOILETTYPE %in% c(4, 5, 12)       ~ "other",
      TOILETTYPE %in% c(95, 96, 98, 99) ~ NA_character_,
      TRUE ~ NA_character_
    )
  )

tabyl(pmaug2022wgps, toilet_cat)

     toilet_cat    n      percent valid_percent
          flush  675 0.0476257673   0.047659394
   improved_pit 5270 0.3718337684   0.372096307
           none  807 0.0569392507   0.056979454
          other   88 0.0062089889   0.006213373
 unimproved_pit 7323 0.5166866577   0.517051472
           <NA>   10 0.0007055669            NA

pmaug2022wgps <- pmaug2022wgps %>%
  mutate(
    toilet_cat_flush          = as.integer(toilet_cat == "flush"),
    toilet_cat_improved_pit   = as.integer(toilet_cat == "improved_pit"),
    toilet_cat_unimproved_pit = as.integer(toilet_cat == "unimproved_pit"),
    toilet_cat_none           = as.integer(toilet_cat == "none"),
    toilet_cat_other          = as.integer(toilet_cat == "other")
  )

#drinking water source

pmaug2022wgps <- pmaug2022wgps %>%
  mutate(
    water_cat = case_when(
      WATERDRINKMAIN %in% c(1, 2)       ~ "piped",
      WATERDRINKMAIN == 3               ~ "public_tap",
      WATERDRINKMAIN %in% c(4, 5, 7, 9) ~ "improved_nonpiped",
      WATERDRINKMAIN %in% c(13, 14, 15) ~ "improved_packaged",
      WATERDRINKMAIN == 16              ~ "improved_nonpiped",
      WATERDRINKMAIN %in% c(6, 8, 12)   ~ "unimproved",
      WATERDRINKMAIN %in% c(10, 11)     ~ "unimproved_vendor",
      WATERDRINKMAIN == 17              ~ "other",
      WATERDRINKMAIN %in% c(96, 98, 99) ~ NA_character_,
      TRUE ~ NA_character_
    )
  )

tabyl(pmaug2022wgps, water_cat)

         water_cat    n      percent valid_percent
 improved_nonpiped 8026 0.5662880124  0.5666878486
 improved_packaged  104 0.0073378960  0.0073430770
             other    6 0.0004233402  0.0004236391
             piped 1260 0.0889014323  0.0889642025
        public_tap 2678 0.1889508220  0.1890842336
        unimproved 2074 0.1463345798  0.1464379016
 unimproved_vendor   15 0.0010583504  0.0010590976
              <NA>   10 0.0007055669            NA

pmaug2022wgps <- pmaug2022wgps %>%
  mutate(
    water_cat_piped              = as.integer(water_cat == "piped"),
    water_cat_public_tap         = as.integer(water_cat == "public_tap"),
    water_cat_improved_nonpiped  = as.integer(water_cat == "improved_nonpiped"),
    water_cat_improved_packaged  = as.integer(water_cat == "improved_packaged"),
    water_cat_unimproved         = as.integer(water_cat == "unimproved"),
    water_cat_unimproved_vendor  = as.integer(water_cat == "unimproved_vendor"),
    water_cat_other              = as.integer(water_cat == "other")
  )

#household assets
#restricted to what IPUMS's SCORE description names: radios, televisions,
#bicycles, furniture. mobphone/motorcycle/electric/car/fridge/hhphone/
#cassette/clock/drawncart are still created here (harmless, may be useful
#elsewhere) but excluded from indicator_vars below since they aren't named
#in the description.

pmaug2022wgps <- pmaug2022wgps %>%
  mutate(
    radio      = as.integer(RADIO == 1),
    tv         = as.integer(TV == 1),
    mobphone   = as.integer(MOBPHONE == 1),
    motorcycle = as.integer(MOTORCYCL == 1),
    electric   = as.integer(ELECTRC == 1),
    bike       = as.integer(BIKE == 1),
    bed        = as.integer(BED == 1),
    cabinet    = as.integer(CABINET == 1),
    chair      = as.integer(CHAIR == 1),
    sofa       = as.integer(SOFA == 1),
    table_own  = as.integer(TABLE == 1),
    car        = as.integer(CAR == 1),
    fridge     = as.integer(FRIDGE == 1),
    hhphone    = as.integer(HHPHONE == 1),
    cassette   = as.integer(CASSETTE == 1),
    clock      = as.integer(CLOCK == 1),
    drawncart  = as.integer(DRAWNCART == 1),
    livestockown = as.integer(LIVESTOCKOWN == 1)
  )

tabyl(pmaug2022wgps, electric)

 electric    n   percent
        0 7404 0.5224017
        1 6769 0.4775983

tabyl(pmaug2022wgps, bike)

 bike    n   percent
    0 9920 0.6999224
    1 4253 0.3000776

#livestock counts
#restricted to what IPUMS's SCORE description names: chickens, sheep,
#goats, horses/camels, cows/bulls, local cattle. PIGNUM and
#OTHERLIVESTOCKNUM are excluded — not named in the description.

livestock_vars <- c("CHICKENNUM", "SHEEPNUM", "GOATNUM", "HORSENUM",
                     "EXOCATTLENUM", "LOCALCATTLENUM")

pmaug2022wgps <- pmaug2022wgps %>%
  mutate(across(all_of(livestock_vars),
                ~ if_else(.x >= 99990, NA_real_, as.numeric(.x))))

#households that own no livestock (LIVESTOCKOWN == 0) were not asked the
#counts and are coded NIU; set their counts to 0 rather than imputing them

pmaug2022wgps <- pmaug2022wgps %>%
  mutate(across(all_of(livestock_vars),
                ~ if_else(is.na(.x) & LIVESTOCKOWN == 0, 0, .x)))

#household-level PCA, split by urban/rural

df_hh <- pmaug2022wgps %>%
  distinct(HHID, YEAR, .keep_all = TRUE) %>%
  mutate(urban_f = as_factor(URBAN)) %>%
  filter(URBAN %in% c(0, 1))

#indicator_vars restricted to exactly what IPUMS's SCORE description
#names: floor/wall/roof materials, water source, toilet facility, radio,
#television, bicycle, furniture (bed/cabinet/chair/sofa/table), and the
#six named livestock categories.

indicator_vars <- c(
  grep("^floor_cat_[a-z]|^wall_cat_[a-z]|^roof_cat_[a-z]|^water_cat_[a-z]|^toilet_cat_[a-z]",
       names(df_hh), value = TRUE),
  "radio", "tv", "bike",
  "bed", "cabinet", "chair", "sofa", "table_own",
  livestock_vars
)

#check remaining missingness before stratum-mean imputation
colSums(is.na(df_hh[indicator_vars]))

          floor_cat_natural       floor_cat_rudimentary 
                          1                           1 
         floor_cat_finished             floor_cat_other 
                          1                           1 
           wall_cat_natural        wall_cat_rudimentary 
                          1                           1 
          wall_cat_finished              wall_cat_other 
                          1                           1 
           roof_cat_natural        roof_cat_rudimentary 
                          1                           1 
          roof_cat_finished              roof_cat_other 
                          1                           1 
           toilet_cat_flush     toilet_cat_improved_pit 
                          1                           1 
  toilet_cat_unimproved_pit             toilet_cat_none 
                          1                           1 
           toilet_cat_other             water_cat_piped 
                          1                           2 
       water_cat_public_tap water_cat_improved_nonpiped 
                          2                           2 
water_cat_improved_packaged        water_cat_unimproved 
                          2                           2 
water_cat_unimproved_vendor             water_cat_other 
                          2                           2 
                      radio                          tv 
                          0                           0 
                       bike                         bed 
                          0                           0 
                    cabinet                       chair 
                          0                           0 
                       sofa                   table_own 
                          0                           0 
                 CHICKENNUM                    SHEEPNUM 
                          5                           2 
                    GOATNUM                    HORSENUM 
                          3                           1 
               EXOCATTLENUM              LOCALCATTLENUM 
                          2                           7

impute_stratum_mean <- function(data, vars) {
  data %>%
    mutate(across(all_of(vars), ~ if_else(is.na(.x), mean(.x, na.rm = TRUE), .x)))
}

df_hh_imputed <- df_hh %>%
  group_by(urban_f) %>%
  group_modify(~ impute_stratum_mean(.x, indicator_vars)) %>%
  ungroup()

run_wealth_pca_imputed <- function(data, vars) {
  mat <- data %>% dplyr::select(dplyr::all_of(vars)) %>% as.matrix()

  var_check <- apply(mat, 2, var, na.rm = TRUE)
  zero_var <- names(var_check)[var_check == 0 | is.na(var_check)]
  if (length(zero_var) > 0) {
    message("Dropping constant/zero-variance columns: ", paste(zero_var, collapse = ", "))
    mat <- mat[, !colnames(mat) %in% zero_var, drop = FALSE]
  }

  pca <- prcomp(mat, center = TRUE, scale. = TRUE)
  list(
    scores = pca$x[, 1],
    loadings = pca$rotation[, 1],
    vars_used = colnames(mat),
    var_explained = summary(pca)$importance[2, 1]
  )
}


urban_hh <- filter(df_hh_imputed, urban_f == "1")   # 1 = Urban
rural_hh <- filter(df_hh_imputed, urban_f == "0")   # 0 = Rural

nrow(urban_hh)

[1] 3348

nrow(rural_hh)

[1] 6033

urban_pca_imputed <- run_wealth_pca_imputed(urban_hh, indicator_vars)

Dropping constant/zero-variance columns: floor_cat_rudimentary

rural_pca_imputed <- run_wealth_pca_imputed(rural_hh, indicator_vars)

#check PCA sign orientation matches across strata before combining
#(using electric here as the check variable, even though electric is
#no longer part of indicator_vars, since it's still a reliable,
#unambiguous wealth marker for validating orientation)

cor_urban_check <- cor(urban_hh$electric, urban_pca_imputed$scores)
cor_rural_check <- cor(rural_hh$electric, rural_pca_imputed$scores)

cor_urban_check

[1] 0.6367513

cor_rural_check

[1] 0.4820995

if (sign(cor_urban_check) != sign(cor_rural_check)) {
  rural_pca_imputed$scores   <- -1 * rural_pca_imputed$scores
  rural_pca_imputed$loadings <- -1 * rural_pca_imputed$loadings
}

loading_compare_imputed <- full_join(
  tibble(indicator = urban_pca_imputed$vars_used, urban_loading = urban_pca_imputed$loadings),
  tibble(indicator = rural_pca_imputed$vars_used, rural_loading = rural_pca_imputed$loadings),
  by = "indicator"
) %>%
  arrange(rural_loading)

print(loading_compare_imputed, n = Inf)

# A tibble: 38 × 3
   indicator                   urban_loading rural_loading
   <chr>                               <dbl>         <dbl>
 1 floor_cat_natural               -0.359         -0.346  
 2 roof_cat_natural                -0.285         -0.302  
 3 wall_cat_natural                -0.239         -0.183  
 4 toilet_cat_unimproved_pit       -0.222         -0.162  
 5 wall_cat_rudimentary            -0.144         -0.159  
 6 toilet_cat_none                 -0.137         -0.146  
 7 water_cat_improved_nonpiped     -0.220         -0.0480 
 8 HORSENUM                        -0.00280       -0.0189 
 9 water_cat_unimproved            -0.0751        -0.0140 
10 floor_cat_rudimentary           NA             -0.00938
11 toilet_cat_other                 0.0227        -0.00611
12 roof_cat_other                  -0.00475        0.00368
13 water_cat_other                  0.00306        0.0114 
14 floor_cat_other                 -0.00340        0.0146 
15 water_cat_unimproved_vendor      0.00105        0.0154 
16 SHEEPNUM                         0.00202        0.0185 
17 water_cat_improved_packaged      0.0365         0.0271 
18 wall_cat_other                  -0.000649       0.0335 
19 LOCALCATTLENUM                   0.0259         0.0347 
20 water_cat_public_tap             0.121          0.0412 
21 roof_cat_finished                0.0229         0.0447 
22 bike                            -0.0644         0.0466 
23 CHICKENNUM                       0.0199         0.0479 
24 GOATNUM                          0.0374         0.0579 
25 EXOCATTLENUM                     0.0370         0.0589 
26 toilet_cat_flush                 0.126          0.0709 
27 water_cat_piped                  0.175          0.0917 
28 chair                            0.0218         0.107  
29 table_own                        0.0650         0.164  
30 radio                            0.114          0.174  
31 bed                              0.141          0.192  
32 cabinet                          0.207          0.229  
33 sofa                             0.207          0.238  
34 toilet_cat_improved_pit          0.176          0.255  
35 tv                               0.281          0.262  
36 roof_cat_rudimentary             0.264          0.295  
37 wall_cat_finished                0.295          0.304  
38 floor_cat_finished               0.359          0.346

select <- dplyr::select
filter <- dplyr::filter

#standardize and form quintiles within stratum and survey year,
#matching IPUMS's within-country-year construction of WEALTHQ

wealth_stratified_imputed <- bind_rows(
  tibble(HHID = urban_hh$HHID, YEAR = urban_hh$YEAR, score_raw = urban_pca_imputed$scores),
  tibble(HHID = rural_hh$HHID, YEAR = rural_hh$YEAR, score_raw = rural_pca_imputed$scores)
) %>%
  left_join(df_hh %>% select(HHID, YEAR, urban_f), by = c("HHID", "YEAR")) %>%
  group_by(urban_f, YEAR) %>%
  mutate(score_z = as.numeric(scale(score_raw)),
         wealthq_stratified_imputed = ntile(score_z, 5)) %>%
  ungroup() %>%
  mutate(wealthq_stratified_imputed = factor(wealthq_stratified_imputed,
                                             labels = c("Lowest","Second","Middle","Fourth","Highest")))

tabyl(wealth_stratified_imputed, YEAR, wealthq_stratified_imputed)

 YEAR Lowest Second Middle Fourth Highest
 2020    598    598    598    596     596
 2021    642    642    642    640     640
 2022    639    639    637    637     637

pmaug2022wgps <- pmaug2022wgps %>%
  left_join(wealth_stratified_imputed %>% select(HHID, YEAR, score_z, wealthq_stratified_imputed),
            by = c("HHID", "YEAR"))


print(loading_compare_imputed, n = Inf)

# A tibble: 38 × 3
   indicator                   urban_loading rural_loading
   <chr>                               <dbl>         <dbl>
 1 floor_cat_natural               -0.359         -0.346  
 2 roof_cat_natural                -0.285         -0.302  
 3 wall_cat_natural                -0.239         -0.183  
 4 toilet_cat_unimproved_pit       -0.222         -0.162  
 5 wall_cat_rudimentary            -0.144         -0.159  
 6 toilet_cat_none                 -0.137         -0.146  
 7 water_cat_improved_nonpiped     -0.220         -0.0480 
 8 HORSENUM                        -0.00280       -0.0189 
 9 water_cat_unimproved            -0.0751        -0.0140 
10 floor_cat_rudimentary           NA             -0.00938
11 toilet_cat_other                 0.0227        -0.00611
12 roof_cat_other                  -0.00475        0.00368
13 water_cat_other                  0.00306        0.0114 
14 floor_cat_other                 -0.00340        0.0146 
15 water_cat_unimproved_vendor      0.00105        0.0154 
16 SHEEPNUM                         0.00202        0.0185 
17 water_cat_improved_packaged      0.0365         0.0271 
18 wall_cat_other                  -0.000649       0.0335 
19 LOCALCATTLENUM                   0.0259         0.0347 
20 water_cat_public_tap             0.121          0.0412 
21 roof_cat_finished                0.0229         0.0447 
22 bike                            -0.0644         0.0466 
23 CHICKENNUM                       0.0199         0.0479 
24 GOATNUM                          0.0374         0.0579 
25 EXOCATTLENUM                     0.0370         0.0589 
26 toilet_cat_flush                 0.126          0.0709 
27 water_cat_piped                  0.175          0.0917 
28 chair                            0.0218         0.107  
29 table_own                        0.0650         0.164  
30 radio                            0.114          0.174  
31 bed                              0.141          0.192  
32 cabinet                          0.207          0.229  
33 sofa                             0.207          0.238  
34 toilet_cat_improved_pit          0.176          0.255  
35 tv                               0.281          0.262  
36 roof_cat_rudimentary             0.264          0.295  
37 wall_cat_finished                0.295          0.304  
38 floor_cat_finished               0.359          0.346

library(factoextra)

Welcome to factoextra!

Want to learn more? See two factoextra-related books at https://www.datanovia.com/library/principal-component-methods

urban_mat <- urban_hh %>% dplyr::select(dplyr::all_of(urban_pca_imputed$vars_used)) %>% as.matrix()
rural_mat <- rural_hh %>% dplyr::select(dplyr::all_of(rural_pca_imputed$vars_used)) %>% as.matrix()

urban_pca_obj <- prcomp(urban_mat, center = TRUE, scale. = TRUE)
rural_pca_obj <- prcomp(rural_mat, center = TRUE, scale. = TRUE)

#urban biplot

fviz_pca_biplot(
  urban_pca_obj,
  select.var = list(contrib = 15),
  geom.ind = "none",
  label = "var",
  col.var = "steelblue",
  repel = TRUE,
  title = "PCA Biplot: Urban Household Wealth Indicators"
)

#rural biplot

fviz_pca_biplot(
  rural_pca_obj,
  select.var = list(contrib = 15),
  geom.ind = "none",
  label = "var",
  col.var = "steelblue",
  repel = TRUE,
  title = "PCA Biplot: Rural Household Wealth Indicators"
)

5. Map out cluster level pregnancy decision-making, specifying which clusters are urban

#map out the scale by EAID

UG<-st_transform(UG_boundary_4, crs=4326)

UG <- st_make_valid(UG)

joined<-st_join(UG, pmaug2022wgps_sf, join=st_intersects)

# Set the global option
options(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)

pmaug2022groupedweighted2 <- joined  %>%
  as_survey_design(ids=EAID, strata = STRATA, weights = FQWEIGHT, nest=TRUE) %>%
  group_by(EAID, YEAR, URBAN) %>%
  summarise(
    summary_score_av_mean = survey_mean(summary_score_av, na.rm = TRUE)
  )



#make map of cluster level pregnancy decision-making

UG_PMA_GPS_2022 <- UG_PMA_GPS_2022 %>%
  mutate(EAID = as.numeric(EAID))

mapdatareg<-left_join(pmaug2022groupedweighted2, UG_PMA_GPS_2022)

Joining with `by = join_by(EAID)`

#convert to a spatial frame

mapdatareg<- st_as_sf(mapdatareg, coords = c("GPSLONG", "GPSLAT"), crs=4326)


mapdatareg <- mapdatareg %>%
  mutate(residence = factor(URBAN, levels = c(0, 1), labels = c("Rural", "Urban")))

#all on one map

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(aes(color = summary_score_av_mean, shape = residence), size = 1.6, alpha = 0.8) +
  scale_color_viridis_c(option = "magma", direction = -1, name = "Pregnancy exercise of choice\nsub-scale mean") +
  scale_shape_manual(values = c(Rural = 15, Urban = 17), name = NULL) +
  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"),
    axis.text    = element_blank(),
    panel.grid   = element_blank()
  )

Warning: Removed 243 rows containing missing values or values outside the scale range
(`geom_sf()`).

#separated by panel

mapdatareg %>%
  filter(!is.na(residence)) %>%
  ggplot() +
  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(aes(color = summary_score_av_mean), size = 1.8, alpha = 0.8) +
  scale_color_viridis_c(option = "magma", direction = -1, name = "Pregnancy exercise of choice\nsub-scale mean") +
  facet_grid(residence ~ 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"),
    axis.text    = element_blank(),
    panel.grid   = element_blank()
  )

6. Map out original and udpated measures of wealth

wealth_clusters <- pmaug2022wgps %>%
  as_survey_design(ids = EAID, strata = STRATA, weights = FQWEIGHT, nest = TRUE) %>%
  group_by(EAID, YEAR) %>%
  summarise(
    wealth_orig_mean = survey_mean(SCORE,   na.rm = TRUE),
    wealth_pca_mean  = survey_mean(score_z, na.rm = TRUE)
  ) %>%
  left_join(UG_PMA_GPS_2022, by = "EAID") %>%
  st_as_sf(coords = c("GPSLONG", "GPSLAT"), crs = 4326)


map_wealth <- function(var, legend_title) {
  ggplot(wealth_clusters) +
    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(aes(color = .data[[var]]), size = 1.8, alpha = 0.8) +
    scale_color_gradient2(low = "#b2182b", mid = "grey85", high = "#2166ac",
                          midpoint = 0, name = legend_title) +
    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"),
      axis.text    = element_blank(),
      panel.grid   = element_blank()
    )
}

p_wealth_orig <- map_wealth("wealth_orig_mean", "Original wealth score\ncluster mean")
p_wealth_pca  <- map_wealth("wealth_pca_mean",  "PCA-derived wealth score\ncluster mean")

p_wealth_orig

p_wealth_pca

7. Regression analysis

#drop missing values



individualandcluster<-left_join(joined, pmaug2022groupedweighted2)

Joining with `by = join_by(EAID, YEAR, URBAN)`

#drop missing values

individualandcluster<-individualandcluster%>%
  filter(MCP!='98')


individualandcluster<-individualandcluster%>%
  filter(educationlevel!='NA')



individualandcluster<-individualandcluster%>%
  filter(maritalcombined!='NA')



individualandcluster<-individualandcluster%>%
  filter(URBAN!='98')

library(modelsummary)



library(lme4)

Loading required package: Matrix


Attaching package: 'Matrix'

The following objects are masked from 'package:tidyr':

    expand, pack, unpack

Registered S3 method overwritten by 'lme4':
  method           from
  na.action.merMod car

# relationship with urban status and summary score: controlling for year and a random effect for individual 

modelbiv <- lmer(summary_score_av ~factor(URBAN)+ factor(YEAR)+ (1|FQINSTID), 
              data = individualandcluster)

modelsummary(
  modelbiv,
  slope = "b",
  stars = TRUE,
  title = "Bivariate relationship between urban/ rural status and pregnancy decision-makiong",
  gof_map = list(list(raw = "nobs", clean = "N", fmt = 0))
)

Bivariate relationship between urban/ rural status and pregnancy decision-makiong
	(1)
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001
(Intercept)	3.879***
	(0.014)
factor(URBAN)1	0.070***
	(0.017)
factor(YEAR)2021	0.061***
	(0.016)
factor(YEAR)2022	0.015
	(0.016)
SD (Intercept FQINSTID)	0.436
SD (Observations)	0.659
N	12150

#add a random effect for cluster


modelbivwithcluster <- lmer(summary_score_av ~ factor(URBAN)+ factor(YEAR)+ (1|EAID)+ (1|FQINSTID), 
              data = individualandcluster)

modelsummary(
  modelbivwithcluster,
  slope = "b",
  stars = TRUE,
  title = "Relationship between urban/ rural status and pregnancy decision-making",
  gof_map = list(list(raw = "nobs", clean = "N", fmt = 0))
)

Relationship between urban/ rural status and pregnancy decision-making
	(1)
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001
(Intercept)	3.896***
	(0.038)
factor(URBAN)1	0.089
	(0.060)
factor(YEAR)2021	0.040**
	(0.015)
factor(YEAR)2022	-0.003
	(0.016)
SD (Intercept FQINSTID)	0.287
SD (Intercept EAID)	0.331
SD (Observations)	0.658
N	12150

Note: adding a random effect for cluster renders urban as insignificant when controlling for year in addition to a random effect for individual.

However, if we build a model that only includes a random effect for individual, not for cluster, urban has a significant and negative on pregnancy decision-making, so the suppression by wealth is happening here.

#full model without the cluster random effect 

library(lme4)

modelmulti <- lmer(summary_score_av ~BIRTHEVENT + educationlevel + agegroup+ maritalcombined+ factor(MCP)+ factor(URBAN)+ factor(WEALTHQ)+ factor(YEAR)+ (1|FQINSTID), 
              data = individualandcluster)

modelsummary(
  modelmulti,
  slope = "b",
  stars = TRUE,
  title = "Relationship between urban/ rural status and pregnancy decision-making, individual random effect only",
  gof_map = list(list(raw = "nobs", clean = "N", fmt = 0))
)

Relationship between urban/ rural status and pregnancy decision-making, individual random effect only
	(1)
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001
(Intercept)	3.708***
	(0.038)
BIRTHEVENT	-0.014**
	(0.005)
educationlevelprimary/middle school	-0.007
	(0.033)
educationlevelsecondary/post-primary	0.210***
	(0.036)
educationleveltertiary/ post-secondary	0.263***
	(0.044)
agegroup20-24	-0.017
	(0.024)
agegroup25-29	0.005
	(0.028)
agegroup30-34	0.047
	(0.033)
agegroup35-39	-0.001
	(0.038)
agegroup40-44	-0.032
	(0.043)
agegroup45-49	-0.029
	(0.046)
maritalcombinedin a union	-0.052**
	(0.018)
factor(MCP)1	0.105***
	(0.016)
factor(URBAN)1	-0.096***
	(0.018)
factor(WEALTHQ)2	0.136***
	(0.023)
factor(WEALTHQ)3	0.222***
	(0.024)
factor(WEALTHQ)4	0.242***
	(0.026)
factor(WEALTHQ)5	0.299***
	(0.028)
factor(YEAR)2021	0.056***
	(0.016)
factor(YEAR)2022	0.026
	(0.016)
SD (Intercept FQINSTID)	0.392
SD (Observations)	0.658
N	12150

#full model with the cluster random effect


modelmulti_EAID <- lmer(summary_score_av ~ BIRTHEVENT + educationlevel + agegroup + maritalcombined +
                          factor(MCP) + factor(URBAN) + factor(WEALTHQ) + factor(YEAR) +
                          (1 | FQINSTID) + (1 | EAID),
                        data = individualandcluster)

modelsummary(
  modelmulti_EAID,
  stars = TRUE,
  title = "Relationship between urban/rural status and pregnancy decision-making (IPUMS wealth) with random effect for EAID",
 gof_map = list(list(raw = "nobs", clean = "N", fmt = 0))
)

Relationship between urban/rural status and pregnancy decision-making (IPUMS wealth) with random effect for EAID
	(1)
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001
(Intercept)	3.733***
	(0.051)
BIRTHEVENT	-0.018***
	(0.005)
educationlevelprimary/middle school	0.078*
	(0.033)
educationlevelsecondary/post-primary	0.269***
	(0.036)
educationleveltertiary/ post-secondary	0.338***
	(0.043)
agegroup20-24	-0.016
	(0.022)
agegroup25-29	0.011
	(0.026)
agegroup30-34	0.057+
	(0.030)
agegroup35-39	0.026
	(0.034)
agegroup40-44	-0.010
	(0.039)
agegroup45-49	-0.004
	(0.042)
maritalcombinedin a union	-0.037*
	(0.017)
factor(MCP)1	0.109***
	(0.015)
factor(URBAN)1	-0.042
	(0.056)
factor(WEALTHQ)2	0.051*
	(0.026)
factor(WEALTHQ)3	0.098***
	(0.027)
factor(WEALTHQ)4	0.114***
	(0.029)
factor(WEALTHQ)5	0.145***
	(0.032)
factor(YEAR)2021	0.037*
	(0.015)
factor(YEAR)2022	0.004
	(0.016)
SD (Intercept FQINSTID)	0.254
SD (Intercept EAID)	0.304
SD (Observations)	0.658
N	12150

Now if we run this with the PCA-derived measure of wealth, urban is not significant.

#run this with the PCA derived measure for wealth

individualandcluster <- individualandcluster %>%
  left_join(wealth_stratified_imputed %>% select(HHID, YEAR, score_z, wealthq_stratified_imputed),
            by = c("HHID", "YEAR"))


modelmulti_pca <- lmer(summary_score_av ~ BIRTHEVENT + educationlevel + agegroup + maritalcombined +
                         factor(MCP) + factor(URBAN) + wealthq_stratified_imputed + factor(YEAR) +
                         (1 | FQINSTID),
                       data = individualandcluster)

modelsummary(
  modelmulti_pca,
  stars = TRUE,
  title = "Relationship between urban/rural status and pregnancy decision-making (PCA-derived wealth)",
 gof_map = list(list(raw = "nobs", clean = "N", fmt = 0))
)

Relationship between urban/rural status and pregnancy decision-making (PCA-derived wealth)
	(1)
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001
(Intercept)	3.690***
	(0.038)
BIRTHEVENT	-0.015**
	(0.005)
educationlevelprimary/middle school	-0.004
	(0.033)
educationlevelsecondary/post-primary	0.219***
	(0.036)
educationleveltertiary/ post-secondary	0.266***
	(0.044)
agegroup20-24	-0.017
	(0.024)
agegroup25-29	0.007
	(0.028)
agegroup30-34	0.051
	(0.033)
agegroup35-39	0.005
	(0.038)
agegroup40-44	-0.028
	(0.043)
agegroup45-49	-0.024
	(0.046)
maritalcombinedin a union	-0.052**
	(0.018)
factor(MCP)1	0.107***
	(0.016)
factor(URBAN)1	-0.017
	(0.017)
wealthq_stratified_imputedSecond	0.132***
	(0.024)
wealthq_stratified_imputedMiddle	0.195***
	(0.024)
wealthq_stratified_imputedFourth	0.226***
	(0.025)
wealthq_stratified_imputedHighest	0.267***
	(0.026)
factor(YEAR)2021	0.058***
	(0.016)
factor(YEAR)2022	0.016
	(0.016)
SD (Intercept FQINSTID)	0.395
SD (Observations)	0.658
N	12150

Now when we add a random effect for cluster, urban remains non-significant.

modelmulti_pca_EAID <- lmer(summary_score_av ~ BIRTHEVENT + educationlevel + agegroup + maritalcombined +
                         factor(MCP) + factor(URBAN) + wealthq_stratified_imputed + factor(YEAR) +
                         (1 | FQINSTID)+ (1|EAID),
                       data = individualandcluster)

modelsummary(
  modelmulti_pca_EAID,
  stars = TRUE,
  title = "Relationship between urban/rural status and pregnancy decision-making (PCA-derived wealth) with random effect for EAID",
   gof_map = list(list(raw = "nobs", clean = "N", fmt = 0))
)

Relationship between urban/rural status and pregnancy decision-making (PCA-derived wealth) with random effect for EAID
	(1)
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001
(Intercept)	3.693***
	(0.052)
BIRTHEVENT	-0.018***
	(0.005)
educationlevelprimary/middle school	0.079*
	(0.033)
educationlevelsecondary/post-primary	0.271***
	(0.036)
educationleveltertiary/ post-secondary	0.340***
	(0.042)
agegroup20-24	-0.018
	(0.022)
agegroup25-29	0.009
	(0.026)
agegroup30-34	0.054+
	(0.030)
agegroup35-39	0.024
	(0.034)
agegroup40-44	-0.013
	(0.039)
agegroup45-49	-0.007
	(0.042)
maritalcombinedin a union	-0.036*
	(0.017)
factor(MCP)1	0.110***
	(0.015)
factor(URBAN)1	-0.000
	(0.055)
wealthq_stratified_imputedSecond	0.115***
	(0.024)
wealthq_stratified_imputedMiddle	0.115***
	(0.026)
wealthq_stratified_imputedFourth	0.140***
	(0.027)
wealthq_stratified_imputedHighest	0.153***
	(0.028)
factor(YEAR)2021	0.039*
	(0.015)
factor(YEAR)2022	-0.000
	(0.016)
SD (Intercept FQINSTID)	0.254
SD (Intercept EAID)	0.305
SD (Observations)	0.658
N	12150

7. How much between-EA variance is explained urban status or wealth?

library(lme4)
library(performance)
library(dplyr)


model_null <- lmer(
  summary_score_av ~ (1|FQINSTID) + (1|EAID),
  data = individualandcluster
)

icc(model_null, by_group = TRUE)

# ICC by Group

Group    |   ICC
----------------
FQINSTID | 0.131
EAID     | 0.177

model_urban_only <- lmer(
  summary_score_av ~ factor(URBAN) + (1|FQINSTID) + (1|EAID),
  data = individualandcluster
)

icc(model_urban_only, by_group = TRUE)

# ICC by Group

Group    |   ICC
----------------
FQINSTID | 0.132
EAID     | 0.176

# ---- Extract variance components from both models ----
var_null  <- as.data.frame(VarCorr(model_null))
var_urban <- as.data.frame(VarCorr(model_urban_only))

eaid_var_null   <- var_null$vcov[var_null$grp == "EAID"]
eaid_var_urban  <- var_urban$vcov[var_urban$grp == "EAID"]

fqinstid_var_null  <- var_null$vcov[var_null$grp == "FQINSTID"]
fqinstid_var_urban <- var_urban$vcov[var_urban$grp == "FQINSTID"]

residual_var_null  <- var_null$vcov[var_null$grp == "Residual"]
residual_var_urban <- var_urban$vcov[var_urban$grp == "Residual"]

# ---- Build a before/after comparison table ----
variance_decomp_table <- tibble(
  component = c("EAID (between-cluster)", "FQINSTID (within-person)", "Residual"),
  variance_before_urban = c(eaid_var_null, fqinstid_var_null, residual_var_null),
  variance_after_urban  = c(eaid_var_urban, fqinstid_var_urban, residual_var_urban)
) %>%
  mutate(
    pct_change = 100 * (variance_before_urban - variance_after_urban) / variance_before_urban
  )

print(variance_decomp_table)

# A tibble: 3 × 4
  component                variance_before_urban variance_after_urban pct_change
  <chr>                                    <dbl>                <dbl>      <dbl>
1 EAID (between-cluster)                  0.111                0.110     0.785  
2 FQINSTID (within-person)                0.0823               0.0823    0.0287 
3 Residual                                0.434                0.434    -0.00285

pct_eaid_explained_by_urban <- 100 * (eaid_var_null - eaid_var_urban) / eaid_var_null

cat(sprintf(
  "Between-EA variance before adding URBAN: %.4f\nBetween-EA variance after adding URBAN:  %.4f\nPercent of between-EA variance explained by urban/rural status: %.1f%%\n",
  eaid_var_null, eaid_var_urban, pct_eaid_explained_by_urban
))

Between-EA variance before adding URBAN: 0.1108
Between-EA variance after adding URBAN:  0.1099
Percent of between-EA variance explained by urban/rural status: 0.8%

reg_vars <- c("summary_score_av", "BIRTHEVENT", "educationlevel", "agegroup",
              "maritalcombined", "MCP", "URBAN", "WEALTHQ",
              "wealthq_stratified_imputed", "YEAR", "FQINSTID", "EAID")

icc_data <- individualandcluster %>%
  filter(if_all(all_of(reg_vars), ~ !is.na(.x)))

get_eaid_var <- function(model) {
  vc <- as.data.frame(VarCorr(model))
  vc$vcov[vc$grp == "EAID"]
}

model_null_common  <- lmer(summary_score_av ~ 1 + (1|FQINSTID) + (1|EAID), data = icc_data)
model_wealth_ipums <- lmer(summary_score_av ~ factor(WEALTHQ) + (1|FQINSTID) + (1|EAID), data = icc_data)
model_wealth_pca   <- lmer(summary_score_av ~ wealthq_stratified_imputed + (1|FQINSTID) + (1|EAID), data = icc_data)

eaid_var_null_common <- get_eaid_var(model_null_common)

wealth_var_compare <- tibble(
  model    = c("Null", "IPUMS wealth quintile", "Stratified PCA wealth quintile"),
  eaid_var = c(eaid_var_null_common,
               get_eaid_var(model_wealth_ipums),
               get_eaid_var(model_wealth_pca))
) %>%
  mutate(pct_eaid_explained = 100 * (eaid_var_null_common - eaid_var) / eaid_var_null_common)

nrow(icc_data)

[1] 12150

nobs(modelmulti)

[1] 12150

print(wealth_var_compare)

# A tibble: 3 × 3
  model                          eaid_var pct_eaid_explained
  <chr>                             <dbl>              <dbl>
1 Null                             0.111                 0  
2 IPUMS wealth quintile            0.0947               14.6
3 Stratified PCA wealth quintile   0.0975               12.0

#get ICC for other covariates

modelmulti_cluster_nourban <- lmer(summary_score_av ~ BIRTHEVENT + educationlevel + agegroup + maritalcombined +
                              factor(MCP) + (1|FQINSTID)+ (1| EAID),
                            data = individualandcluster)
library(performance)
icc(modelmulti_cluster_nourban)

# Intraclass Correlation Coefficient

    Adjusted ICC: 0.273
  Unadjusted ICC: 0.262

icc(modelmulti_cluster_nourban, by_group = TRUE)

# ICC by Group

Group    |   ICC
----------------
FQINSTID | 0.109
EAID     | 0.163

var_covariates <- as.data.frame(VarCorr(modelmulti_cluster_nourban))
eaid_var_covariates <- var_covariates$vcov[var_covariates$grp == "EAID"]

pct_eaid_explained_by_covariates <- 100 * (eaid_var_null - eaid_var_covariates) / eaid_var_null

cat(sprintf(
  "Percent of between-EA variance explained by full individual covariate set: %.1f%%\n",
  pct_eaid_explained_by_covariates
))

Percent of between-EA variance explained by full individual covariate set: 12.1%

#cumulative between-EA variance explained: covariates, + wealth, + urban

cumulative_formulas <- list(
  "Covariates"                        = summary_score_av ~ BIRTHEVENT + educationlevel + agegroup + maritalcombined + factor(MCP),
  "Covariates + IPUMS wealth"         = summary_score_av ~ BIRTHEVENT + educationlevel + agegroup + maritalcombined + factor(MCP) + factor(WEALTHQ),
  "Covariates + IPUMS wealth + urban" = summary_score_av ~ BIRTHEVENT + educationlevel + agegroup + maritalcombined + factor(MCP) + factor(WEALTHQ) + factor(URBAN),
  "Covariates + PCA wealth"           = summary_score_av ~ BIRTHEVENT + educationlevel + agegroup + maritalcombined + factor(MCP) + wealthq_stratified_imputed,
  "Covariates + PCA wealth + urban"   = summary_score_av ~ BIRTHEVENT + educationlevel + agegroup + maritalcombined + factor(MCP) + wealthq_stratified_imputed + factor(URBAN)
)

cumulative_models <- lapply(cumulative_formulas, function(f) {
  lmer(update(f, . ~ . + (1 | FQINSTID) + (1 | EAID)), data = icc_data)
})

cumulative_var <- tibble(
  model    = c("Null", names(cumulative_models)),
  eaid_var = c(eaid_var_null_common, sapply(cumulative_models, get_eaid_var))
) %>%
  mutate(pct_eaid_explained = 100 * (eaid_var_null_common - eaid_var) / eaid_var_null_common)

print(cumulative_var)

# A tibble: 6 × 3
  model                             eaid_var pct_eaid_explained
  <chr>                                <dbl>              <dbl>
1 Null                                0.111                 0  
2 Covariates                          0.0974               12.1
3 Covariates + IPUMS wealth           0.0925               16.5
4 Covariates + IPUMS wealth + urban   0.0926               16.4
5 Covariates + PCA wealth             0.0925               16.5
6 Covariates + PCA wealth + urban     0.0932               15.9