We used secondary data that was extracted from the Bangladesh Demographic and Health Survey (BDHS) 2017–2018. It was the latest cross-sectional household survey throughout the country. We used 7235 Bangladeshi adult women aged 18–49 years and their last-born under-five children as samples. The BDHS collected household, socio-demographic, lifestyle, and health-related information of mothers and their under-five children from October 2017 to March 2018. Moreover, BDHS 2017–18 measured the height and weight of the selected women and their under-five children. The study population, sample, study design, questionnaire, instruments, data collection procedure, and data reliability were described elsewhere [28].

Bangladeshi non-pregnant married women, living in Bangladesh, aged 18–49 years, and having at least one under-five child living with mothers (who were eligible for height and weight measurements) were considered as samples for the analysis.

BDHS 2017–18 used two-stage stratified cluster sampling for selecting households from Bangladesh. In the first stage, 675 enumeration areas (EAs) (250 in urban and 425 in rural areas) were selected by stratified sampling with proportional allocation. In the second stage, 30 households were selected from each selected EA using systematic sampling. BDHS 2017–18 eliminated three EAs due to communication problems and finally considered 672 EAs and 20,160 households for the survey. They mentioned that the sampling weights were not expected to lead to any significant differences in the overall survey indicators [28]. For the present study, we first considered 8,653 women with at least one under-five child who were eligible for height and weight measurements. BDHS 2017–18 considered one child if a woman had twin babies. According to our exclusion criteria, we excluded some women and their under-five children. Data were checked, and the outliers of the dataset, missing values, and incomplete data were excluded. Finally, 7235 women and their last-born under-five children were considered for the present study (Fig. 1).

Sample selection procedure for the present study

There were two outcome variables for the study: (i) nutritional status of mothers, measured by their body mass index (BMI), where BMI = weight (kg) / ((height (m))2. Mothers were defined as having chronic energy deficiency if their BMI was < 18.5 kg/m2, normal weight (18.5 ≤ BMI < 25 kg/m2), and over-nutrition if their BMI was ≥ 25 kg/m2 [29]; (ii) the nutritional status of under-five children was measured by three indicators: (i) stunting (height-for-age, z-score below − 2), (ii) underweight (weight-for-age, z-score below − 2), and (iii) wasting (weight-for-height, z-score below − 2). Each indicator was classified into two classes according to the cut-off point suggested by WHO; stunting (stunted: code 1, not stunted: code 0); underweight (underweight: code 1, not underweight: code 0); and wasting (wasted: code 1, not wasted: code 0) [30].

First, the main independent variable age at the first marriage (year) (AAFM) was divided into five groups; (i) AAFM ≤ 15, (ii) AAFM ≤ 16─<18, (iii) AAFM ≤ 18─≤20, (iv) AAFM ≤ 21─≤24, and (v) AAFM ≥ 25. Finally, it was divided into two groups according to the rule of the Bangladesh government: (i) child marriage (AAFM < 18 years) and (ii) not child marriage (AAFM ≥ 18 years) [28]. Child marriage was categorized as yes (code, 1), and not child marriage as no (code, 0). Some socioeconomic, demographic, and household factors were also considered independent variables in this study, as mentioned in Table 1. We followed some previous studies for selecting the variables and preparing their categories [28, 31, 32].

The background characteristics of the samples were summarized using a frequency distribution. We determined the proportion of each category of outcome variables and women’s child marriage using frequency distribution, and descriptive statistics was used to calculate the mean ± SD and median of AAFM of women and the mean ± SD of women and their under-five children age. The chi-square (χ2) test was utilized to examine the significance of the association between women’s child marriage and the nutritional status of women and their under-five children. The analysis of variance (ANOVA) was also used to find the variation in mean BMI among AAFM groups. As we mentioned, BDHS 2017–2018 collected data from overall Bangladesh using two-stage stratified cluster sampling; the data came from different levels of hierarchy. There was a cluster effect on the data set; a single-level statistical model would not be appropriate for analyzing this type of data set [33]. In this study, two-level logistic regression analysis was used for accounting cluster level variation and examining the association between mothers’ child marriage and the nutritional status of mothers and their under-five children. In two-level logistic regression, children/women was at the unit level and cluster at the second level. The cluster level variation was calculated using the formula mentioned in Nakagawa et al. study [34]. One of our outcome variables (women’s nutritional status) was ordinal, where ordinal logistic regression was more appropriate to analyze the data. However, the test of Parallel Lines showed that the assumption of the model was not satisfied. Alternatively, the two-level multinomial logistic regression model was used to find the association of child marriage with mothers’ nutritional status, uncontrolling/controlling the effect of selected socioeconomic, demographic, and other factors. We also used a two-level binary logistic regression model to examine the association between mothers’ child marriage and the nutritional status of their under-five children, uncontrolling/controlling the effect of the selected socioeconomic, demographic, and other factors. The variance inflation factor (VIF) was used to examine the multicollinearity problem among independent variables in logistic models; if 0 < VIF < 5, it was judged that there was no evidence of multicollinearity [35]. We did not find any multicollinearity problems among the independent variables for both multinomial and binary logistic models. The chi-square goodness-of-fit test and the Hosmer and Lemeshow test were used to examine the goodness-of-fit of multinomial and binary logistic regression models, respectively. Moreover, the accuracy of the models was checked by the receiver operating characteristic (ROC) curve. The adjusted odds ratio (AOR, OR = eβ, where β is the regression coefficient) with a 95% confidence interval of AOR and p-value were utilized to interpret the results coming from logistic regression models. Statistical significance was accepted at p < 0.05. Statistical analyses were carried out using SPSS software (version IBM 20).