Data analysis

Basic descriptive statistics were used to describe the data and the respondents’ characteristics (Table 1) and the prevalence of mothers’ and daughters’ FGM (Table 2). The bivariate binary logistic regression was used to compare the odds of daughters’ genital mutilation across the countries (Table 3). We also computed the prevalence of daughters’ FGM among mothers that had FGM and those without FGM across the countries (Table 3). Finally, the binary multivariable multilevel logistic regression models were used to identify the contribution of the individual, community and country-level factors associated with the FGM among female children (Table 4). We used the binary multivariable model embedded in the MLWin 3.03 module in Stata statistical package version 16. The first order marginal quasi-likelihood linearization (MQL1) was adopted for the estimation algorithms using the iterative generalized least squares (IGLS). We applied sampling weights to the data, and statistical significance was determined at 5%.

All numbers were weighted.

AOR Adjusted odds ratio, CI confidence interval, MOR median odds ratio, VPC variance partition coefficient. Tanzania has no records or religion and was excluded from the analysis. **Had all independent variables in the model.

However, Tanzania and Niger were not included in the multivariable analysis because the countries did not capture data on respondents’ religion and if FGM was a religions’ requirement or not. Respondents who did not give birth to any female child, or who has no living daughter or who were selected for male questionnaire according to the survey protocol were not asked questions on daughter’s FGM. Of the 95,507 women that were asked questions on daughter’s FGM, 93,063 (97,8%) provided valid responses (S1 Table). The remaining 2.2% with invalid responses on FGM among daughters were excluded from further analysis.

In an attempt to arrive at the best model that best predict mothers’ practice of FGM on their children, binary multilevel binary logistic regression models were fitted. We identified and fitted five models, and the results of these models were reported in Table 3. The first model was the null model (Model I) to assess the variation due to the neighbourhood- and the country-specific random effects without any explanatory variable. Model II had only the individual-level variables conditional on the neighbourhood and country-specific random effects while Model III had only the neighbourhood level variables conditional on the neighbourhood and country-specific random effects. In contrast, Model IV examined the country-level variables conditional on the neighbourhood and country-specific random effects. Finally, we developed Model V to estimate the odds of individual, neighbourhood and country-level factors conditional on the neighbourhood and country-specific random effects.

Typically, the multilevel analysis provides estimates for both fixed and random effects. The fixed effects are the lowest hierarchy in the multilevel model (individual mothers’ characteristics in this study) while the random effects are based on the contribution of the upper levels (higher hierarchies), these are the neighbourhood and country levels in this study. A three-level variance model for π_ijk= P(y_ijk = 1) as shown in Eqs (1) and (2) was determined.

In the Binomial function, there are two possible possibilities: had FGM for at least one daughter or not; y_ijk is the daughter i of neighbourhood j from country k, while the probability that the mother of daughter i of neighbourhood j from country k had FGM for at least a daughter is denoted by π_ijk.

The fitted models were based on the hierarchical logistic regression model with mixed outcomes consisting of the fixed and random parts, as shown in Eq (2).

Where πijk1−πijk is the odds that y_ijk = 1, β₀ is the fixed intercept, β_p are the regression coefficients of the covariates X_p, U_ojk is the random effect of daughters in the neighbourhood j in the country k, and U_ojk is the random effect of country k, e_ijk is the noise such that

Other statistical details of the model have been published earlier [29–31]. Note that p(β0)∝1,p(σU2)∼Γ−1(ε,ε) and p(σV2)∼Γ−1(ε,ε). We assumed diffuse priors and used an improper uniform prior for β₀ and a commonly used conjugate inverse Gamma prior for σU2 and σV2. We implemented the regression Eq (3) in the MLWin 3.03 module embedded in Stata statistical package version 16 with the individuals (i′s) as level 1, the neighbourhoods (j′s) as level 2 and countries (k′s) as level 3.

Fixed effects. The primary outcomes of all the five models were the measures of the association expressed as odds ratios (ORs) with their 95% confidence intervals (CI). It is an expression of the likelihood of the outcome variable in the categories of a variable compared to the reference category of that variable.

Random effects. This is a measure of variation that is explained by the higher levels of the hierarchies in the data. The measures of variations were explored using the variance partition coefficient (VPC) and the median odds ratio (MOR). We reported the random effects in terms of the odds using the methods proposed by Larsen et al. on neighbourhood effects [32]. The VPC is a summary of the degree of clustering in the data and it is a reasonable interpretation of the Intra-class correlation (ICC) which measures the extent to which the y_ijk′s in the same neighbourhoods and countries resemble each other as compared to those from other clusters [33]. Therefore, the VPC was used to measure the proportion of the total variance which are accounted for at the neighbourhood and the country levels and computed as σU2/(σU2+σV2+σe2) and σV2/(σU2+σV2+σe2) respectively.

The MORs are the measures of the variance of the odds ratio in higher levels (neighbourhood or country), and it estimates the probability of daughters’ FGM that can be attributed to any of the neighbourhood and country factors. A MOR of 1 is an indication that there is no neighbourhood or country variability. A MOR > 1 suggests that the contextual effects -neighbourhood/country variability- is significant. A higher MOR indicates that the contextual effects for understanding the probability of a mother to have FGM performed for her daughter is higher. This statistical analytic approach has been used and reported in the literature [29–31,34,35].

Multicollinearity within explanatory variables was assessed by examining the variance inflation factor (VIF) [36]. There was no evidence of multicollinearity. We used the model deviance computed from the -2loglikelihood to assess and identify the best model that fitted the data. Lower deviance indicates a better model fit.