Show the code
library(haven)
library(tidyverse)
library(dplyr)
library(survival)
library(survminer)
library(ggsurvfit)Examining Multigenerational Households and Fertility Using Survival Analysis
Abstract
Household structure plays an important role for individuals and families, and there are many factors that determine the living arrangements of household members. Geography, demography, and specific social and cultural characteristics also play an important role in the living arrangements of individuals. Lastly, economic uncertainty also plays an important role in determining living arrangements, specifically in multigenerational households. Acknowledging that these differences exist, the question then that this paper proposes to explore relates to the benefits, if any, of multigenerational households on increased fertility. Specifically, the research question to be answered is: do multigenerational households increase the likelihood of second births; and does this difference vary by rural or non-rural classification? To answer this question, a semi-parametric survival analysis is conducted utilizing data from the DHS Program for the country of Bangladesh for the years 1993-1994 and 2017-2018.
Literature
Like Many Countries, Bangladesh has experienced fertility declines. In 1993, which is one of the observation years for this analysis, the fertility rate of Bangladesh was 4 children per women. In 2017, that number dropped to 2.1. Studies (Roy & Hossain, 2017, Bora et.al, 2020, ) have examined fertility differentials of Bangladeshi women and found that education is a strong predictor of decreased fertility, which is not surprising as similar studies have determined the same relationship between education and fertility across different parts of the world. It is important to note that Bangladesh is listed by the United Nations as a least developed country. Economic challenges, while extremely significant, are not factored into this analysis directly (more on future research later). In addition to the economic challenges that Bangladesh faces, the country is also notable for their dramatic reduction in fertility rates. In 1994, the World Bank published “The Determinants of Reproductive Change in Bangladesh: Success in a Challenging Environment”, which credited family planning programs as the main mechanism for Bangladesh’s dramatic fertility decline since the 1960s. The demographer John C. Caldwell and his colleagues, however, disagreed with the World Bank’s analysis in their article, “The Bangladesh Fertility Decline: An Interpretation” and pronounced that family planning programs alone would not be sufficient to cause such a fertility decline. Instead, through their study, which consisted of an experimental anthropological-demographic approach, involved interview data of Bangladeshi individuals and family units in rural areas. It should be noted that some of these families were multigenerational households, but, that was not factorized into their analysis. The authors concluded that from the familial and individual perspective, Bangladeshis insinuated that socioeconomic conditions also led to their decision to have less children, therefore providing an alternative interpretation to the World Bank’s description of what led to the decline in fertility in that country.
Other studies on fertility trends in Bangladesh have examined the impact of other covariates, such as a mother’s age, whether contraception has ever been used, the death of a child at any time, whether the woman has ever worked, religion, region of residence, and female independence (Kahn & Raeside, 1977). A woman’s status, as a position of power, has also been measured to determine if there is an effect on fertility decision for women in rural Bangladesh (Balk, 1994). Researchers have also examined women’s empowerment through micro-crediting programs and their effect on contraceptive use and fertility decisions (Steele, et.al., 1994). Studies continue to analyze fertility rates in Bangledesh, assessing the impact of individual and community level variables. However, research on the impact of household structure, and in particular, multigenerational households and fertility are scarce. The analysis presented in this paper will contribute to the body of literature that has, since the 1960s, focused on the decline of fertility in Bangladesh and sought to find explanations as to what contributed to their dramatic decline. A survival analysis approach will differentiate this study from much of the research that has been published examining fertility in Bangladesh
Data and Methodology
For this study, a semi-parametric model will be created. The Cox proportional hazards regression model will provide a hazard ratio for second births based on two covariates: multigenerational households and rural or non-rural status. The following formula represents the model:
Two models will be constructed for the two different sets of years using the covariates mentioned earlier (multigenerational households and rural vs. not rural). This will allow us to make observations over two decades. Data for this analysis was obtained from the DHS Program. Data for Balgnedesh for the years, 1993-1994 and 2017-2018 were used. Two DHS data sets were joined to be able to bring in the multigenerational variable and the second birth variable from the household and individual files, respectively.
To successfully execute a cox proportional hazards regression model, the following assumptions must be met: a constant hazard ratio over time, a multiplicative effect or relationship between the covariates and the hazard, and independence of survival times between individuals. Therefore, testing for the proportionality assumption was conducted.
First, the cox.zpg function was used, which scaled the Schoenfeld residuals for “rural2” and “new_mghh”. The results showed a global chi-square that is not significant, as it is close to 0. Generally, we want to remain closer to zero, which would indicate that there has been no violation of proportional hazards. Secondly, the results of the test were plotted, and we can observe a flat pattern for both covariates, which is a good visual indicator of proportional hazards. Figure 1 below shows the residuals for the multigenerational household variable. A similar pattern was found for the rural variable. Assumptions for performing a Cox Proportional Hazards Regression were therefore met.
Results & Discussion
The results from the two cox regression models were combined into one table displaying the hazard ratios for each of the models.
Holding the other variable constant, the hazard ratio of having a second birth for rural households in the 1993-94 cohort was 1.15, indicating a strong relationship between being a rural household and having a second birth, compared to the observation group of non-rural or urban households (95% confidence interval: 1.057 –1.2447). The difference between the two groups is significant, with a resulting p-value of 0.000989. Interestingly, the hazard of having a second birth for multigenerational households are 15% lower (HR=by compared to the observation group (not a multi-generational household). The difference between both groups is significant, with a resulting p-value of 1.64e-10 . The confidence intervals at 95% are 0.8038 and 0.8906 for the multigenerational group results.
Similar results are observed for the 2017-18 data for the rural covariate, where the hazard of having a second birth are higher for rural household types. However, the results between both year cohorts differ for the multigenerational households covariate. While the hazard of having a second birth for multigenerational households are still lower, it is no longer significant with a p-value of 0.055. No explanation for the difference is offered, but it is noteworthy that studies have looked at increases of labor force participation and poverty reduction programs as important influences for turning the tide on Bangledeshi poverty.
An additional analysis was conducted as a supplement to the Cox model. The Aalen’s additive regression model provides plots that are informative regarding the effect of covariates over time. That is to say, the Aalen’s additive model allows for time varying covariate effects. The additive hazards were plotted for the 1993-1994 cohort data. The results show that over time, household type (rural vs. nonrural) has an increasing effect, whereas multigenerational status household has a decreasing effect. Therefore, whether it is a mutligenrerationl household or not, becomes less influential on the hazard of having a second birth.
A similar pattern is observed for the 2017-18 dataset, however the effect of rural vs. non-rural increase sharply before leveling out.
There are theoretical perspectives, especially from the field of anthropology, that have sought to frame the role of grandmothers in evolution. The “grandmother hypothesis” posits that grandmothers allow women to have more babies by offering help with care giving. A recent study by Anthropologist Kristen Hawkes revisited the grandmother hypothesis in an attempt to locate reproductive benefits, if any, of living with grandmothers, or by proxy, a multigenerational household. While this study did not directly challenge that hypothesis, the results from the Bangledesh dataset suggest that the grandmother hypothesis does not apply as there is no increase in hazard of having a second birth for women in multigenerational households. Additionally, this analysis opened up a dimension of understanding fertility by looking at variables outside of those commonly examined, including women’s age, education, empowerment, economic access and more.
Limitations and Conclusion
In this analysis, the impact of living in a multigenerational hazard on second births was examined for Bangladesh using two data sets separated by over two decades from 1993-94 and 2017-18. Theoretically, this paper sought to challenge existing explanations for the decline in fertility in Bangledesh that have primarily credited institutional actions (via family planning programs). In addition, an anthropological theory that posits a positive relationship between grandmothers and fertility was briefly discussed to provide additional context as to why we might consider the impact of multigenerational households on fertility viz-a-viz second births. Results from a cox proportional regression model indicate that, contrary to the grandmother hypothesis, multigenerational households do not increase the probability of women having second births. This finding was true for both year cohorts, but significant only for the 1933-1994 cohort. Adding to the existing literature that examines the many factors that have led to fertility declines in Bangladesh since 1960, this analysis provides enough evidence to further examine other household level factors that might influence fertility decisions. For example, while this analysis took into consideration multigenerational households, it did not account specifically for grandmothers in the household versus grandfathers. Fertility decisions are not made in a vacuum. Macro, meso, and micro level factors play an important role. While some relationships, such as education and decreased fertility, have been reaffirmed over time, others haven’t been sufficiently studied. Multigenerational households, and further examining the role grandparents are playing in fertility decisions that are leading to differentials across demographic groups will be important to study, to better understand the composition of the future population; and therefore determine what policies and supports need to be in place for all individuals and chidlren to thrive.
R-Analysis, Code and Steps Taken
Load Libraries
Bangladesh, 1993-94
Reading in the data for Bangladesh: Standard DHS, 1993-94.
Now, I bring in the Individual Recode file for Bangladesh and join them together.
Show the code
ind <- read_dta("C:/Users/codar/OneDrive - University of Texas at San Antonio/R 2022/BDIR31FL.DTA")
ind$hh <- substr(ind$caseid, 1,12)
bang93<- left_join(ind, bgh2, by=c("hh" ="hhid"))
# names(mdat)Show the code
ba93 <- bang93 %>%
filter(bidx_01==1&b0_01==0)%>%
transmute(CASEID=caseid,
int.cmc=v008,
fbir.cmc=b3_01,
sbir.cmc=b3_02,
marr.cmc=v509,
rural=v025,
educ=v106,
age = v012,
agec=cut(v012, breaks = seq(15,50,5), include.lowest=T),
partneredu=v701,
partnerage=v730,
weight=v005/1000000,
psu=v021,
strata=v022,
mgenhh = mgenhh,
hh = hh) %>%
select(CASEID, int.cmc, fbir.cmc, sbir.cmc, marr.cmc, rural, educ, age, agec, partneredu, partnerage, weight, psu, strata, mgenhh, hh, rural)%>%
mutate(agefb = (age - (int.cmc - fbir.cmc)/12))
#Calculating the birth intervals, observed and censored and the event indicator (did women have a second birth?)
ba93a<-ba93%>%
mutate(secbi = ifelse(is.na(sbir.cmc)==T,
int.cmc - fbir.cmc,
fbir.cmc - sbir.cmc),
b2event = ifelse(is.na(sbir.cmc)==T,0,1))Organizing Covariates: Multigenerational HHs and Urban/Rural
Show the code
#Recoding mgenhh from char to numeric
lookup <- c("notmgen" = 0, "mgenhh" = 1)
ba93a$new_mghh <- lookup[ba93a$mgenhh]
#recoding rural to 0/1 binary
ba93a$rural2 <- ifelse(ba93a$rural==2, 1, 0)
library(survey)
#survey design
options(survey.lonely.psu = "adjust")
des<-svydesign(ids=~psu, strata=~strata,
data=ba93a, weight=~weight)Fit the Model and Run Tests
Show the code
library(gtsummary)
# library(splines)
cox.ba <-svycoxph(Surv(secbi, b2event)~rural2+new_mghh, design=des)%>% tbl_regression(exponentiate = TRUE)Stratified 1 - level Cluster Sampling design (with replacement)
With (301) clusters.
svydesign(ids = ~psu, strata = ~strata, data = ba93a, weight = ~weight)Show the code
cox.ba1 <-svycoxph(Surv(secbi, b2event)~rural2+new_mghh, design=des)
summary(cox.ba1)Stratified 1 - level Cluster Sampling design (with replacement)
With (301) clusters.
svydesign(ids = ~psu, strata = ~strata, data = ba93a, weight = ~weight)Call:
svycoxph(formula = Surv(secbi, b2event) ~ rural2 + new_mghh,
design = des)
n= 8478, number of events= 6943
coef exp(coef) se(coef) robust se z Pr(>|z|)
rural2 0.13723 1.14709 0.03796 0.04167 3.294 0.000989 ***
new_mghh -0.16710 0.84611 0.02440 0.02614 -6.392 1.64e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
rural2 1.1471 0.8718 1.0571 1.2447
new_mghh 0.8461 1.1819 0.8038 0.8906
Concordance= 0.521 (se = 0.004 )
Likelihood ratio test= NA on 2 df, p=NA
Wald test = 53.69 on 2 df, p=2e-12
Score (logrank) test = NA on 2 df, p=NA
(Note: the likelihood ratio and score tests assume independence of
observations within a cluster, the Wald and robust score tests do not).The Wald Test provides insight into the overall significance of the model. Concordance: We want something closer to 1.
Testing to see if proportionality assumption is met/Testing for proportional hazards
Show the code
cox.zph(cox.ba1) chisq df p
rural2 0.000474 1 0.98
new_mghh 0.000278 1 0.99
GLOBAL 0.000694 2 1.00Show the code
plot11 <-plot(cox.zph(cox.ba1)) 
Show the code
plotfunction (x, y, ...)
UseMethod("plot")
<bytecode: 0x000002ca74ebd848>
<environment: namespace:base>“An assumption of CPH regression is that the hazard (think risk) associated with a particular variable does not change over time.” I will test for proportional hazards in two steps. First, I will use the cox.zpg function, which will scale the Schoenfeld residuals for “rural2” and “new_mghh”. The results show a global chi-square that is not significant, as it is close to 0. Generally, we want to remain closer to zero, which would indicate that there has been no violation of proportional hazards. Secondly, I plot the results of the test, and we can observe a flat pattern, which is a good visual indicator of proportional hazards. Assumptions for performing a Cox Proportional Hazards Regression are therefore met.
Show the code
fitz <- aareg(Surv(secbi, b2event) ~rural+new_mghh, ba93a, weights = weight)
summary(fitz) slope coef se(coef) z p
Intercept 0.03510 0.000215 1.65e-05 13.00 1.17e-38
rural 0.00406 0.000032 8.49e-06 3.77 1.61e-04
new_mghh -0.00628 -0.000039 5.98e-06 -6.53 6.73e-11
Chisq=62.87 on 2 df, p=2.23e-14; test weights=aalenShow the code
library(ggfortify)
autoplot(fitz)Warning: `gather_()` was deprecated in tidyr 1.2.0.
ℹ Please use `gather()` instead.
ℹ The deprecated feature was likely used in the ggfortify package.
Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.
Over time, household type (rural vs. nonrural) has an increasing effect, whereas multigenerational status household has a decreasing effect. So, over time, wether it is a mutligenrerationl household or not, becomes less influential on the hazard of having a second birth.
Bangladesh, 2017-18
Aalens Additive Hazards
Show the code
fitz2 <- aareg(Surv(secbi, b2event) ~rural0+mghh, b18a, weights = weight)
summary(fitz2) slope coef se(coef) z p
Intercept 0.024000 1.05e-04 2.05e-06 51.30 0.00e+00
rural0 0.004590 1.99e-05 2.24e-06 8.89 5.90e-19
mghh -0.000641 -4.34e-06 2.17e-06 -2.00 4.51e-02
Chisq=82.36 on 2 df, p=<2e-16; test weights=aalenShow the code
library(ggfortify)
autoplot(fitz2)
Stratified 1 - level Cluster Sampling design (with replacement)
With (672) clusters.
svydesign(ids = ~psu, strata = ~strata, data = b18a, weight = ~weight) Length Class Mode
table_body 29 broom.helpers list
table_styling 7 -none- list
N 1 -none- numeric
n 1 -none- numeric
N_event 1 -none- numeric
model_obj 26 svycoxph list
inputs 10 -none- list
call_list 15 -none- list Stratified 1 - level Cluster Sampling design (with replacement)
With (672) clusters.
svydesign(ids = ~psu, strata = ~strata, data = b18a, weight = ~weight)Call:
svycoxph(formula = Surv(secbi, b2event) ~ rural0 + mghh, design = des1)
n= 17940, number of events= 13612
coef exp(coef) se(coef) robust se z Pr(>|z|)
rural0 0.17992 1.19712 0.01940 0.02604 6.909 4.88e-12 ***
mghh -0.03874 0.96200 0.01722 0.02017 -1.921 0.0547 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
rural0 1.197 0.8353 1.1376 1.260
mghh 0.962 1.0395 0.9247 1.001
Concordance= 0.519 (se = 0.003 )
Likelihood ratio test= NA on 2 df, p=NA
Wald test = 53.46 on 2 df, p=2e-12
Score (logrank) test = NA on 2 df, p=NA
(Note: the likelihood ratio and score tests assume independence of
observations within a cluster, the Wald and robust score tests do not).| Characteristic | 1993-94 | 2017-18 | ||||
|---|---|---|---|---|---|---|
| HR1 | 95% CI1 | p-value | HR1 | 95% CI1 | p-value | |
| rural2 | 1.15 | 1.06, 1.24 | <0.001 | |||
| new_mghh | 0.85 | 0.80, 0.89 | <0.001 | |||
| rural0 | 1.20 | 1.14, 1.26 | <0.001 | |||
| mghh | 0.96 | 0.92, 1.00 | 0.055 | |||
| 1 HR = Hazard Ratio, CI = Confidence Interval | ||||||
Future studies could look into more current rate of fertility based on economic developments ?? Regression analysis to control for income?
Sources
Deborah Balk (1994) Individual and Community Aspects of Women’s Status and Fertility in Rural Bangladesh, Population Studies, 48:1, 21-45, DOI: 10.1080/0032472031000147456
Caldwell, J. C., Barkat-e-Khuda, Caldwell, B., Pieris, I., & Caldwell, P. (1999). The Bangladesh Fertility Decline: An Interpretation. Population and Development Review, 25(1), 67–84. http://www.jstor.org/stable/172372
H.T.Abdullah Khan, Robert Raeside (1997)Factors affecting the most recent fertility rates in urban-rural Bangladesh,Social Science & Medicine,Volume 44, Issue 3,1997,Pages 279-289,ISSN 0277-9536,https://doi.org/10.1016/S0277-9536(96)00076-7.
Jayanta Kumar Bora, Nandita Saikia, Endale Birhanu Kebede & Wolfgang Lutz (2022) Revisiting the causes of fertility decline in Bangladesh: the relative importance of female education and family planning programs, Asian Population Studies, DOI: 10.1080/17441730.2022.2028253
Roy, S., Hossain, S.M.I. Fertility differential of women in Bangladesh demographic and health survey 2014. Fertil Res and Pract 3, 16 (2017). https://doi.org/10.1186/s40738-017-0043-z
Steele, Fiona, Sajeda Amin, and Ruchira Tabassum Naved. 1998. “The impact of an integrated micro-credit programme on women’s empowerment and fertility behavior in rural Bangladesh,” Policy Research Division Working Paper no. 115. New York: Population Council