Chapter 2- pregnancy decision-making, age, union status and contraceptive use
Author
R. Luttinen
Pregnancy Decision-making
This chapter explores panel data from the performance monitoring for action sample in Uganda collected from 2020-2022. I employ random effects for individual to study the pregnancy exercise of choice scale over time. I pay special attention to interactions between age group and marital status as well as age group and usage of modern contraception.
#read in datalibrary(ipumsr)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
#read in long format#person-period datasetpmauglong <-read_ipums_micro(ddi ="C:/Users/Rebecca/Downloads/pma_00029.xml",data ="C:/Users/Rebecca/Downloads/pma_00029.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.
#filter for 2022 only#pmauglong<-pmauglong%>%# filter(YEAR==2021)
#recoding and filtering#check distribution of variables of interestlibrary(janitor)
Attaching package: 'janitor'
The following objects are masked from 'package:stats':
chisq.test, fisher.test
#marital status pmauglongfilter2<-pmauglongfilter2%>%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(pmauglongfilter2, maritalstatus)
maritalstatus n percent valid_percent
never married 3251 0.2358873893 0.23593875
currently married 3384 0.2455376578 0.24559112
currently living with partner 4934 0.3580031926 0.35808114
formerly in union 1897 0.1376433029 0.13767327
widow or widower 313 0.0227107822 0.02271573
<NA> 3 0.0002176752 NA
#age-grouppmauglongfilter2<-pmauglongfilter2 %>%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" ))pmauglongfilter2<-pmauglongfilter2%>%filter(BIRTHEVENT<90)pmauglongfilter2<-pmauglongfilter2 %>%mutate(AGE=as.numeric(AGE))pmauglongfilter2<-pmauglongfilter2 %>%mutate(URBAN=as.factor(URBAN))#contraceptive usepmauglongfilter2<-pmauglongfilter2%>%mutate(MCP=as.factor(MCP))%>%mutate(usingmoderncon=recode(MCP, '1'="yes" ,'0'="no",.default =NA_character_))tabyl(pmauglongfilter2, usingmoderncon)
usingmoderncon n percent valid_percent
no 8955 0.64990202 0.6594256
yes 4625 0.33565571 0.3405744
<NA> 199 0.01444227 NA
educationlevel n percent valid_percent
none 797 0.0578416431 0.05785004
primary/middle school 7666 0.5563538718 0.55643464
secondary/post-primary 4165 0.3022715727 0.30231545
tertiary/ post-secondary 1149 0.0833877640 0.08339987
<NA> 2 0.0001451484 NA
#collapse marital status variablepmauglongfilter2<-pmauglongfilter2%>%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"))
The following objects are masked from 'package:tidyr':
expand, pack, unpack
#filter out cases with a missing wealth index pmauglongfilter3<-pmauglongfilter2%>%filter(WEALTHQ!=99)#filter out any infertile or missing values for fertility preferencespmauglongfilter3<-pmauglongfilter3%>%filter(!FERTPREF %in%c(3, 97, 98, 99))#begin multivariate testing for the average score#base modellm_model <-lmer(summary_score_av ~ agegroup + maritalcombined + usingmoderncon+ educationlevel+ BIRTHEVENT+factor(URBAN)+ SCORE+factor(YEAR)+factor(FERTPREF)+ (1|FQINSTID) , data = pmauglongfilter3)summary(lm_model)
Linear mixed model fit by REML ['lmerMod']
Formula: summary_score_av ~ agegroup + maritalcombined + usingmoderncon +
educationlevel + BIRTHEVENT + factor(URBAN) + SCORE + factor(YEAR) +
factor(FERTPREF) + (1 | FQINSTID)
Data: pmauglongfilter3
REML criterion at convergence: 26130.3
Scaled residuals:
Min 1Q Median 3Q Max
-4.0836 -0.4252 0.0967 0.5978 2.7144
Random effects:
Groups Name Variance Std.Dev.
FQINSTID (Intercept) 0.1511 0.3887
Residual 0.4294 0.6553
Number of obs: 11516, groups: FQINSTID, 6930
Fixed effects:
Estimate Std. Error t value
(Intercept) 3.889354 0.038342 101.439
agegroup20-24 -0.017550 0.024465 -0.717
agegroup25-29 0.004917 0.028544 0.172
agegroup30-34 0.060063 0.033340 1.802
agegroup35-39 0.023171 0.038986 0.594
agegroup40-44 0.007372 0.044944 0.164
agegroup45-49 0.033609 0.049520 0.679
maritalcombinedin a union -0.056571 0.018614 -3.039
usingmodernconyes 0.114311 0.015783 7.242
educationlevelprimary/middle school 0.001430 0.033537 0.043
educationlevelsecondary/post-primary 0.215301 0.037116 5.801
educationleveltertiary/ post-secondary 0.236454 0.045574 5.188
BIRTHEVENT -0.014452 0.005486 -2.634
factor(URBAN)1 -0.113544 0.018709 -6.069
SCORE 0.040521 0.003627 11.173
factor(YEAR)2021 0.069118 0.015994 4.322
factor(YEAR)2022 0.014716 0.016350 0.900
factor(FERTPREF)2 -0.039629 0.022720 -1.744
Correlation matrix not shown by default, as p = 18 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need it
#interaction between union status and age grouplm_modelint <-lmer(summary_score_av ~ agegroup + maritalcombined + usingmoderncon+ educationlevel+ BIRTHEVENT+ agegroup*maritalcombined+factor(URBAN)+ SCORE+factor(YEAR)+factor(FERTPREF)+ (1|FQINSTID), data = pmauglongfilter3)summary(lm_modelint)
Linear mixed model fit by REML ['lmerMod']
Formula: summary_score_av ~ agegroup + maritalcombined + usingmoderncon +
educationlevel + BIRTHEVENT + agegroup * maritalcombined +
factor(URBAN) + SCORE + factor(YEAR) + factor(FERTPREF) +
(1 | FQINSTID)
Data: pmauglongfilter3
REML criterion at convergence: 26131.6
Scaled residuals:
Min 1Q Median 3Q Max
-4.0958 -0.4241 0.0924 0.6012 2.7227
Random effects:
Groups Name Variance Std.Dev.
FQINSTID (Intercept) 0.1499 0.3872
Residual 0.4296 0.6554
Number of obs: 11516, groups: FQINSTID, 6930
Fixed effects:
Estimate Std. Error t value
(Intercept) 3.914472 0.038989 100.399
agegroup20-24 -0.060843 0.032363 -1.880
agegroup25-29 -0.103069 0.041473 -2.485
agegroup30-34 0.013144 0.049292 0.267
agegroup35-39 0.060883 0.058314 1.044
agegroup40-44 0.018586 0.061355 0.303
agegroup45-49 -0.032120 0.067546 -0.476
maritalcombinedin a union -0.184068 0.040244 -4.574
usingmodernconyes 0.117610 0.015801 7.443
educationlevelprimary/middle school 0.001528 0.033520 0.046
educationlevelsecondary/post-primary 0.211516 0.037118 5.699
educationleveltertiary/ post-secondary 0.239650 0.045628 5.252
BIRTHEVENT -0.015686 0.005518 -2.843
factor(URBAN)1 -0.111441 0.018698 -5.960
SCORE 0.040169 0.003630 11.066
factor(YEAR)2021 0.068837 0.015994 4.304
factor(YEAR)2022 0.015245 0.016348 0.933
factor(FERTPREF)2 -0.038943 0.022707 -1.715
agegroup20-24:maritalcombinedin a union 0.156064 0.051567 3.026
agegroup25-29:maritalcombinedin a union 0.241218 0.057604 4.188
agegroup30-34:maritalcombinedin a union 0.160018 0.062918 2.543
agegroup35-39:maritalcombinedin a union 0.054901 0.068598 0.800
agegroup40-44:maritalcombinedin a union 0.083634 0.072247 1.158
agegroup45-49:maritalcombinedin a union 0.193261 0.079360 2.435
Correlation matrix not shown by default, as p = 24 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need it
library(ggeffects)#visualize the interactionpredictions <-predict_response(lm_modelint, terms =c("agegroup", "maritalcombined"))predictions %>%as.data.frame() %>%mutate(x =fct_relevel(x, "15-19", "20-24", "25-29", "30-34", "35-39", "40-44", "45-49")) %>%ggplot(aes(x = x, y = predicted, color = group, group = group)) +geom_point(position =position_dodge(0.4), size =2.5) +geom_errorbar(aes(ymin = conf.low, ymax = conf.high), position =position_dodge(0.4), width =0.2) +labs(title ="Marginal Effect of Age Group and Marital Status on Summary Score",x ="Age Group",y ="Predicted Summary Score Average",color ="Marital Status" ) +theme_minimal()
#interaction between contraceptive use and age grouplm_modelint2 <-lmer(summary_score_av ~ agegroup + maritalcombined + BIRTHEVENT+ usingmoderncon+ educationlevel+agegroup*usingmoderncon+ SCORE+factor(URBAN)+factor(FERTPREF)+factor(YEAR)+ (1|FQINSTID), data = pmauglongfilter3)summary(lm_modelint2)
#make a list to put into model summary# 3. Put your models into a named list for the table columnsmodels_list_con <-list("Model 1: Main Effects"= lm_model,"Model 2: Age × Marital"= lm_modelint,"Model 3: Age × Contraception"= lm_modelint2)
library(modelsummary)library(marginaleffects)
Attaching package: 'marginaleffects'
The following object is masked from 'package:lme4':
refit
modelsummary( models_list_con,slope ="b",stars =TRUE,title ="Regression Table: Average Marginal Effects",gof_map =NA)
