Bayesian Multilevel Analysis of Utilization of Antenatal Care Services in Ethiopia

In sub-Saharan Africa, 72% of pregnant women received an antenatal care visit at least once in their pregnancy period. Ethiopia has one of the highest rates of maternal mortality in sub-Saharan African countries. So, this high maternal mortality levels remain a major public health problem. According to EDHS, 2016, the antenatal care (ANC), delivery care (DC), and postnatal care (PNC) were 62%, 73%, and 13%, respectively, indicating that ANC is in a low level. The main objective of this study was to examine the factors that affect the utilization of antenatal care services in Ethiopia using Bayesian multilevel logistic regression models. The data used for this study comes from the 2016 Ethiopian Demographic and Health Survey which was conducted by the Central Statistical Agency (CSA). The statistical method of data analysis used for this study is the Bayesian multilevel binary logistic regression model in general and the Bayesian multilevel logistic regression for the random coefficient model in particular. The convergences of parameters are estimated by using Markov chain Monte-Carlo (MCMC) using SPSS and MLwiN software. The descriptive result revealed that out of the 7171 women who are supposed to use ANC services, 2479 (34.6%) women were not receiving ANC services, while 4692 (65.4%) women were receiving ANC services. Moreover, women in the Somali and Afar regions are the least users of ANC. Using the Bayesian multilevel binary logistic regression of random coefficient model factors, place of residence, religion, educational attainment of women, husband educational level, employment status of husband, beat, household wealth index, and birth order were found to be the significant factors for usage of ANC. Regional variation in the usage of ANC was significant.


Background of the Study
Maternal health refers to the well-being of women through pregnancy, childbirth, and the postpartum period. In rich countries, where women have access to basic health care, giving birth is a positive and fulfilling experience. On the other hand, for many women that live in poor countries, it is associated with suffering, ill-health, and even death [1]. Maternal health is advantageous to the communities, families, and nation due to its profound effects on the health of women, immediate contact of survival of the newborn, and long-term well-being of children, particularly girls and the well-being of families.
The recent World Health Organization (WHO) estimates on maternal mortality showed that developed countries had a consistent low maternal mortality ratio that averaged less than 10 deaths per 100,000 live births for over a decade [2]. Studies showed that Poland, Sweden, Japan, United Kingdom, and New Zealand have had the best performance in reducing and maintaining low maternal mortality levels over the past decades [3].
In sub-Saharan Africa, there are still more maternal health care problems. It is approximated that about 800 girls and women died every day as a result of pregnancy and childbirth-related complications [4]. Maternal deaths are now increasingly concentrated in sub-Saharan Africa, where high fertility rates combine with inadequate access to quality antenatal care and skilled attendance at birth to substantially elevate the risk of death in this region. The WHO estimates that over 500,000 women and girls die from complications of pregnancy and childbirth each year, worldwide, with approximately 99% of these deaths occurring in developing countries [5].
In sub-Saharan Africa, countries like Burundi, Cameroon, and Senegal registered about 740, 590, and 320 deaths per 100,000 live births, respectively. Variations in maternal health outcomes and any issues in sub-Saharan Africa have attributed to pregnancy and childbirth complications with unbalanced socioeconomic status [6].
The use of antenatal care services among women in sub-Saharan Africa has shown that 72% of pregnant women received an antenatal care visit once or more times. The very low maternal and infant morbidity and mortality rates reported for developed countries compared with the high figures in developing countries have been attributed to the higher utilization of modern maternal health services by the former [7].
Ethiopia has one of the highest rates of maternal mortality in sub-Saharan African countries. So, these high maternal mortality levels remain a major public health problem. The proportion of women reproductive age from 15 to 49 reports having problems in accessing health care services, i.e, the services are decreased from 96% in 2005, 94% in 2011, and 70% in 2016, despite the country's quest to reduce maternal mortality of ratio to 199 per 100,000 live births by 2020 based on the suggestion made with the goal of the reproductive health program. According to the Ethiopian Demographic and Health Survey, the antenatal care (ANC), skilled delivery care, and postnatal care (PNC) were 62%, 73%, and 13%, respectively, indicating that ANC is low [8].
Considering the goal of reducing maternal mortality and improving maternal health care services, it is essential to understand the factors that affect antenatal care service utilization in order to find the possible solutions for improving maternal health care service in Ethiopia.
Studies conducted in Nigeria and Ethiopia tried their effort to establish correlates having an effect on antenatal care service utilization using binary logistic regression [9]. A study conducted by [9] which used logistic regression to study the utilization of maternal health care services in Kenya found out a number of socioeconomic and demographic characteristics of women greatly influence the usage of maternal health care services such as age of women, education, work status, place of residence, birth order, wealth index, and region. However, their study did not investigate the relative contribution at the individual and regional level of factors influencing women's use of health care services using the Bayesian multilevel.
Despite the country's mission to reduce maternal mortality ratio to 199 maternal deaths per 100,000 live births by 2020 on the basis of the recommendation made with the goal of the reproductive health program, the ratio is still very high [8]. This gap needs to be investigated to ensure an increase in the utilization of maternal health care among mothers through identifying independent variables influencing the utilization of antenatal care services, hence a decrease in child mortality and improvement in maternal care. To fill this gap, the Bayesian multilevel modeling has been employed for identifying factors affecting antenatal care service utilization in Ethiopia at individual and regional levels.
There are a number of reasons for using multilevel models. The classical logistic regression analysis treats the units of analysis as independent observations. One consequence of not considering hierarchical structures is that standard errors of regression coefficients will be underestimated, leading to an overstatement of statistical significance and wrong conclusion. Standard errors for the coefficients of higher-level predictor variables will be the most affected by ignoring grouping.
In a situation where research questions concern the extent of grouping in individual outcomes, and the identification of "outlying" groups, multilevel analysis is highly advised. For evaluating health care intervention, obtaining the effect of antenatal care on maternal mortality can be examined. Additionally, for estimating group effects simultaneously with the effects of group-level predictors, the use of multilevel analysis is mandatory. Moreover, making an inference to a population of groups can be attained using the multilevel analysis [10]. The multilevel analysis enables the proper investigation of the effects of all explanatory variables measured at different levels on the response variable. They employed the model on HIV testing and got a more accurate result [11].
The 2016 Ethiopian Demographic and Health Survey (EDHS) sample was a stratified cluster and selected in two stages [8]. This stratified cluster sampling scheme often introduces multilevel dependency or correlation among the observations that can have implications for model parameter estimates. For multistage clustered samples, the dependence among observations often comes from several levels of the hierarchy. In this case, the use of single-level statistical models is no longer valid and reasonable. Hence, in order to draw appropriate inferences and conclusions from multistage stratified clustered survey data, we may require techniques like multilevel modeling [12].
The purpose of this study is to understand the current status of utilization of antenatal care services in Ethiopia and to elucidate the various explanatory variables having an effect on the use of ANC services in Ethiopia. Moreover, the study is aimed at comparing the predictive performance of the Bayesian multilevel and the classical binary logistic regression models. The results of the study would help the policymakers in understanding the determinants of maternal and child mortality in the country and serve as an important input for any possible intervention aimed at improving the low utilization of antenatal care services in the country. This study covered women in the reproductive age from 15 to 49 years who gave the outcome of antenatal care service in Ethiopia.

Study Area and Variable
Description. Ethiopia has 9 regional states and two city administrations with 611 districts and 15,000 lower administrations. According to the 2016 Ethiopia Demographic and Health Survey (EDHS), the total population of Ethiopia was around 102 million with 4.5 (1) Dependent Variable. The outcome of interest in this study was the utilization of ANC service. A woman is considered to have used ANC if she was checked by a health professional at least once during her pregnancy period. Therefore, the outcome is a categorical variable which is measured as use or not use of service care.
With P ij = P ðY ij = 1/X ij Þ is the probability of using service care for woman i in region j, X ij is the observed characteristic of the i th women in the j th region and 1 − P ij = P ðY ij = 0/X ij Þ is the probability of not using service care for woman i in region j.
(2) Independent Variable. The utilization of antenatal care services is influenced by socioeconomic and demographic factors. These include age of respondent employment status of women, employment status of the husband, birth order, husband's education status, place of residence, religion, attitudes towards wife-beating, women's decision-making power, marital status, mothers' educational status, women exposure to media, and wealth index.

Method of Statistical Data Analysis
2.2.1. Multilevel Linear Models. When the data was collected in hierarchical or clustered structures, the suitable model is multilevel models. Multilevel models are used to account for the correlation of observations within a given group by incorporating group-specific random effects. These random effects can be nested (individuals nested in regions, with random effects at the women and region levels) [13]. The dependent variable must be examined at the lowest and highest level of analysis.

Bayesian Multilevel Logistic Regression Model
(1) Two-Level Model. A multilevel logistic regression model can account for lack of independence across levels of nested data (i.e., individuals nested within regions). It assumes that all experimental units are independent meaning that any variable which affects the utilization of antenatal care services has the same effect in all regions, but multilevel models are used to assess whether the effect of predictors varies from region to region. The probability of "success" or "failure" is the same for all individuals in the group.
The standard assumption is that Y ij has a Bernoulli distribution. Like the logistic regression, the p ij is modeled using the link function logit. The two-level models are given by where X hij = ðX 1ij , X 2ij , ⋯ ⋯ :,X kij Þ represents the first and the second level of covariates for variable k, β = ðβ 0 , β 1 , ⋯ ⋯ ⋯ ⋯ :,β k Þ are regression parameter coefficients, and U 0j , U 1j , ⋯ ⋯ ⋯ ⋯ U kj are the random cluster effects of a model parameter at higher levels with the assumption that follows a normal distribution with mean zero and variance σ 2 u . Therefore, conditional on U 0j , U 1j , ⋯ ⋯ ⋯ U kj , the Y ij can be assumed to be independently distributed as Bernoulli random variables.
(i) Bayesian Multilevel Analysis of Empty Models. The empty two-level model also called a null two-level model for a dichotomous outcome variable refers to a population of groups (level-2units) and specifies the probability distribution for group dependent probabilities p j in Y ij =p j +ε ij without taking further explanatory variables into account. We focus on the model that specifies the transformed probabilities f ðp j Þ and the logit function; then, f ðp j Þ is just the log-odds for group j. Thus, the log-odds have a normal distribution in the population of groups, and it is given by where p j = e β 0 +u 0j /ð1 + e β 0 +u 0j Þ, β 0 is the population average of the transformed probabilities, and u 0j are regional level random effects that are i.i.d. normally distributed with zero means and constant variances σ 2 u . This model does not include a separate parameter for the level-one variance, because the level-one residual variance of the Y ij follows the Bernoulli distribution directly from the success probability, as indicated by Varðϵ ij Þ = p j ð1 − p j Þ.
The likelihood function of the empty model is 3 Computational and Mathematical Methods in Medicine (a) Prior Distribution. For the parameters β 0 and σ 2 u0 , the prior distributions are default noninformative uniform distribution for the intercept and gamma for the random part of the model. It is given by P ðβ 0 Þ ∝ 1 and P ðσ −2 u0 Þ ∝ gamma ðα, βÞ where α and β are scale and shape parameters which are fixed constants.
(b) Posterior Distribution. By using MCMC methods, we can estimate random parameters from the posterior distribution P ðσ 2 u0 , β 0 /y ij Þ.
For parameter β 0 , the full conditional distribution is For the parameter σ 2 u0 , the full conditional distribution is denoted as where N = ∑nj is the total number of individual respondents interviewed in all regions of the country and α and β are scale and shape parameters, respectively, which are both fixed constants.
(2) Variance Component Model. The variance component model has the form where Y ij is the utilization of health care services of the i th women in j th region; u 0j is regional level random effects that are i.i.d. normally distributed with mean zero and constant variances σ 2 u . ϵ ij is errors that are i.i.d. normally distributed with zero means and constant variances σ 2 e and β 0 is the overall mean.
The model decomposes the total variance into two the region and women levels, representing between and within region variabilities in the utilization of antenatal care services [14]. The interclass correlation (ICC) measures the correlation between observations within cluster (region) as where σ 2 u is the "between-region" variance and σ 2 ε is the "between-mothers" variance.
(i) Bayesian Multilevel Analysis of Random Intercept Models. In the random intercept logistic regression model, the intercept is the only random effect meaning that the groups differ with respect to the average value of the response variable, but the relation between explanatory and response variables cannot differ between groups. It represents the heterogeneity between groups in the overall response.
We now assume that there are variables which are a potential explanation for observed success or failure. These variables are denoted by The logistic random intercept model expresses the logodds, i.e., the logit of p ij , as a sum of a linear function of the explanatory variable and a random group-dependent deviation u 0j .
where β 0 +∑ k h=1 β h x hij is the fixed part of the model, because the coefficients are fixed, the remaining part u 0j is known as the random part of the model [15] and the success probability is For the prior and posterior distribution, it is also the same as in the above equation.
(ii) Bayesian Multilevel Analysis of Random Coefficient Models. Random coefficient models are ones where the coefficient of lower-level predictors is modeled as well. The random coefficient model explains unobserved heterogeneity in the effects of explanatory variables on the response variable.
Consider explanatory variables which are potential explanations for the observed outcomes. Denoting these variables by X h =fX hij, i = 1, 2, ⋯:n j , h = 1, 2, ⋯ ⋯ , k, and j = 1, 2, ⋯ ::, N,g. Since some or all of these variables could be levelone variable, the success probability is not necessarily the same for all individuals in a given group. Therefore, the success probability depends on the individual as well as the group and is denoted by p ij .
The first part of the above model, β 0 + ∑ k h=1 β hj x hij , is the fixed part, ðu 0j, u 1j Þ follows a bivariate normal distribution, Computational and Mathematical Methods in Medicine (a) Prior Distribution. Let us denote the parameters β 0 , β 1 , ⋯ ⋯ ⋯ ⋯ :β k and Ω u as prior distribution as follows.
Pðβ 0 Þ ∝ 1, Pðβ 1 Þ, ⋯ ⋯ ⋯ ⋯ Pðβ k Þ ∝ 1 and pðΩ u Þĩ nverse-Wishartk ðs u , hÞ distribution, where Ω u is the variance-covariance matrices, s u is an estimate for the true value of Ω u , and k is the dimension of Ω u that denotes the degree of freedom, so this prior is essentially as uninformative prior. In statistics, the Wishart distribution is a generalization to multiple dimensions of the chi-squared distribution or, in the case of noninteger degrees of freedom, of the gamma distribution [16]. The Wishart distribution is the sampling distribution of the matrix of sums of squares and products of normal distributional assumption. Then, Ω u is positive definite with probability density function given by (b) Posterior Distribution. For the parameter β 0 , β 1 , β 2 , ⋯ ⋯ ⋯ ⋯ :,β k , the posterior distribution is given by whereĥ = 0, 1, 2, ⋯::, k. For the parameter Ω u ,the full conditional distribution is given by (c) Estimation Techniques for MCMC. Bayesian estimation was performed using the software MLwiN, a specialized program for performing multilevel analysis. The Markov chain Monte-Carlo (MCMC) method is a general simulation method for sampling from posterior distributions and computing posterior quantities of interest.
For Bayesian model selection, the deviance information criterion (DIC) that is a hierarchical modeling generalization of the AIC is used. The definition of deviance is a -2loglikelihood ratio of a reduced model compared to the full model. In Bayesian, the lowest expected deviance has the highest posterior probability. Assessing goodness of fit involves investigating how close the values are predicted by the model with that of observed values [17]. The comparison of observed to predicted values using the likelihood function is based on the statistic called deviance.
where p i is the predicted value for observation i. The larger deviance, the poorer fitness of the model.  Table 1 illustrate the relationship between socioeconomic and demographic variables with the use or not use of antenatal care services. Based on the P value in the last column, all the variables have an association with the use or not use of antenatal care.

Bayesian Multilevel Logistic Regression Analysis Output.
In the multilevel analysis, a two-level structure is used with regions as the second-level units and women as the first-level units. The expectation is that there would be differences in the usage of ANC services among regions. The nesting structure is women within regions with a total of 7171 women. Table 2, the deviance-based chi-square test of the random coefficient of the Bayesian multilevel model is statistically significant (P value = 0.000 * ) and the deviance-based chi-square test is 290.47. This implies that as compared to the model with a random intercept with the fixed explanatory variables and null Bayesian multilevel model, the Bayesian multilevel model with random coefficient is better fit. Therefore, the study has been focused on the Bayesian multilevel logistic regression for the random coefficient model.

Comparisons of Model. As it is depicted in
From Table 3, the measure D indicates the average deviance from the complete set of iterations. DðθÞ indicates the deviance at the expected value of the unknown parameters, and pD is the estimated degrees of freedom consumed in the fit or it is the difference between the average deviance from the complete set of iterations and the deviance at the expected value of the unknown parameters. Also, the deviance information criterion is the sum of D and pD. For comparing models, when we use the DIC (deviance information criteria). The random intercept model reduced by 942.39 than the null model. Thus, the random intercept model with the fixed explanatory variables shows a highly significant improvement suggesting that it is better than the null model and the DIC diagnostics of the random coefficient model is reduced by 950.9 and 8.51 from the random intercept model with the fixed explanatory variables and the null model, respectively. Therefore, this smallest value of DIC indicates that Bayesian multilevel logistic regression for the random coefficient model is the most significant model as compared to the null and random intercept model with the fixed explanatory variables. Therefore, this study employs the Bayesian multilevel logistic regression for the random coefficient model. The null model contains no explanatory variables, and it can be focused on assessing the heterogeneity of utilization of antenatal care services among regions. As it is displayed in Table 4, the multilevel empty model fits better than the single-level empty model for antenatal care services. The result from Table 4 indicates that the variance of the random factor is significant which indicates that there is regional variation in receiving antenatal care services.
From the result presented in Table 4 above, it is observed that the overall mean of receiving antenatal care services   Computational and Mathematical Methods in Medicine from a health professional is estimated as β0 = 1:215 and the regional variance σ 2 u0 is estimated as 2.830. From the result presented in Table 4 above, it is observed that the overall mean of receiving antenatal care services from a health professional is estimated as β0 = 1:215 and the regional variance σ 2 u0 is estimated to be 2.830. The variances σ 2 ε and σ 2 u in Table 4 estimate the variation among individual mother and among regions of the country. The individual (level-1) variance was fixed to π 2 /3 ð3:29Þ for the logit model. Also, to examine how much variation in the utilization of ANC services of women was attributable to the region level factors, it is useful to see the intraregion correlation coefficient ðICCÞ = σ 2 u /ðσ 2 u + σ 2 ε Þ = 2:830/ð2:830 + 3:29Þ = 0:4624, which measures the proportion of variance of receiving antenatal services that is between regions, not within regions. The intraregion correlation coefficient (ICC) in the intercept-only model is 0.4624 which is significant at a 5% level of significance. This means that 46.24% of the variation in utilization of ANC service exists between regions, whereas the remaining is 53.76% attributable to an individual level, i.e., within region differences.

Bayesian Multilevel Logistic Regression of the Random
Coefficient Model. This model contains a random slope for the educational attainment of women, which means that it allows the effect of the coefficient of this variable to vary from region to region. This model is more appropriate than the random intercept model for the variables being used since it is intuitive to assume that the educational attainment of women varies from region to region.
As shown in Table 5, the Bayesian multilevel random coefficient model, some of the independent variables have a significant effect on the usage of ANC services in Ethiopia. These are the place of residence, beat, birth order, religious, husband occupation, wealth index, and husband education.
For the categorical variable place of residence in Table 5, the odds of using antenatal care service for women living in rural areas is 0.206 times than women living in urban areas. In other words, the use of antenatal care services for women living in rural areas was 79.4% less than that of women in urban areas. This indicates that women who were living in rural areas were less likely to receive ANC services than living in urban areas. The educational level of the husband was also found to be a significant contributing factor for receiving ANC services at a 5% level of significance. The odds of the educational level of the husband who had primary education were 1.547 indicating that the usage of antenatal care services for those women in this category is 54.75% more than the women having a husband with no education. Similarly, women having a husband with secondary education and higher-level education have better odds of using antenatal care services than women having a husband with no education. Mothers whose husband's educational levels were primary, secondary, and higher level were more likely to receive ANC services than women whose husbands were not educated.
The birth order in Table 5 is the other significant factor for receiving ANC services at a 5% level of significance. The odds of using antenatal care service for birth order category of women having children between 5 and 8 were 0.795 indicating a 20.5% decrease than women having children from 1 to 4. Similarly, women having children greater than eight has an odd ratio of 0.719 indicating a 28.1% decrease in the use of antenatal care services compared to women having 1-4 children. In other words, women having children between 5 and 8 and greater than eight were less likely to receive ANC services than women having children from one up to four.
The effect of attitudes towards wives beating on the use of antenatal care services was the other concern of this study illustrated in Table 5. From the analysis, it was found that the variable was found to be a significant determinant variable of receiving ANC services. The odds of using antenatal care services for beating wives that argues with husband were 0.864 indicating a 13.6% decrease compared to not beating wives that argues with husband.
The household wealth index displayed in Table 5 was also the other significant factor that affects the use of antenatal care services at a 5% level of significance. The odds of using ANC for women in the medium wealth index category were 1.571 which indicates that women in the medium wealth index category use antenatal care services 57.1% higher than women who are in the poor index category. Similarly, those women in the rich wealth index category use the antenatal care services 90% higher than the women in the poor wealth index category. This indicates that women who are in the medium and rich wealth index category were more likely to receive ANC service care as compared to women in the poor wealth index category.
The husband's occupation in Table 5 was found to be a significant determinant variable of receiving ANC services.

Computational and Mathematical Methods in Medicine
The odds of using antenatal care services for women having a working husband were 1.421 showing a 42.1% increase in usage than women whose husband was not working. In other words, women whose husband was working were more likely to receive ANC services than women whose husband was not working.
The other concern of this study was to examine the effect of religion on the use of antenatal care services. As it is depicted in Table 5, women whose religion was Catholic, Protestant, Muslim, and other followers have odds ratios of 0.335, 0.493, 0.623, and 0.252, respectively. This indicates that women whose religion are Catholic, Protestant, Muslim, and others use antenatal care services 66.5%, 50.7%, 37.7%, and 74.8%, respectively, less than women with Orthodox religion. Thus, women that follow Catholic, Protestant, Muslim, and other followers were less likely to receive ANC services than women whose religion was Orthodox.
From Table 6, we can clearly observe that the educational attainment of women has significant contributions on receiving ANC services at a 5% level of significance. The odds of educational attainment of women who had primary, secondary, and higher education were 2.052, 3.408, and 10.612, respectively. This points out that the use of antenatal care services for women with primary, secondary, and higher education was 2.052, 3.408, and 10.612 times than women who had no education, respectively. This illustrates that women whose educational levels were primary secondary and higher level were more likely to attend ANC services than women who were not educated.
The estimated variance of intercept and slope of educational attainment varies significantly. This showed that the factor women educational level has brought a variation in the usage of antenatal care services across regions of the country. Some of the variances of the interaction terms between intercepts and slopes of explanatory variables are also found significant. Interpretation of significant covariance terms can be easily made in terms of the correlation coefficients displayed in Table 6, and the correlation matrix contains the estimated correlation between random intercepts and slopes which are negatively correlated. The negative sign for the correlation between intercepts and slopes implies that regions with higher intercepts tend to have on average lower slopes on the corresponding predictors.
3.6. Binary Logistic Regression Analysis Output. Binary logistic outcomes (dependent variables) are binary (dichotomous) and can be coded 0 (no use of ANC) and 1 (use of ANC).
The predictors can be continuous or dichotomous, as in regression analysis. The output of the binary logistic regression is depicted in Table 7.
Among the covariates in the model, residence, birth order, religion, husband occupation, wealth index, and husband education have a significant effect on the odds of using antenatal care services. But the remaining covariate respondent's opinion that a husband is justified in hitting or beating his wife and age of the respondents do not have a significant effect on the odds of using antenatal care services at a 5% level of significance.

Comparison of Bayesian Multilevel Random Coefficient
and Binary Logistic Models. The goodness of fit of various statistical models can be examined using various statistical measures. These measures are used to select a suitable model among the possible candidate models. In Table 8, some model comparison measures have been listed. Among these, the most commonly used statistical measures used to select a model having good predictive capacity are AIC and BIC. The model having the smallest AIC (DIC) and BIC is considered to have a good fit and predictive performance that can minimize the difference between the actual and predicted observations.

Computational and Mathematical Methods in Medicine
On the basis of the AIC and DIC values in Table 8, the multilevel random coefficient model has a good fit or better predictive performance than the binary logistic regression model as both the DIC for it is much smaller compared to the AIC for the binary logistic model. The other comparison measures also confirmed that the multilevel random coefficient model has a good fit.

Discussion
This study employed the Bayesian multilevel logistic regression analysis. Women are nested within the various regions in Ethiopia in order to explain regional variations in the usage of antenatal services. We employed three multilevel logistic regression models for the response variable on the use or not use of antenatal care services. The Bayesian multilevel logistic regression empty model, Bayesian multilevel logistic regression random intercept, and Bayesian multilevel logistic regression for the random coefficient model were applied to explain regional differences in the usage of antenatal care services among women.
The place of residence was significantly associated with the usage of antenatal care services. This shows that women who were living in rural areas were less likely to receive ANC services than living in urban areas. This difference might be due to the fact that urban women are more accessible to health services and have information and education about ANC service care than women in the rural area. The findings from this study and previous studies in Ethiopia showed that there exists an association between ANC utilization and residence of mothers [18][19][20].
In this study, it was found that women in primary, secondary, and higher education were more likely to receive antenatal care services than women who are not educated. The study that was conducted in Ghana also showed a similar finding with this study [21]. Moreover, another study conducted in Assosa, Ethiopia, concluded that women's education is the determinant factor for the use of antenatal care services [22]. Thus, this study also confirmed that women who attended formal education are better in using antenatal care services than those who are not educated.
This study also concluded that the husband's educational level were a significant contributing factor for receiving ANC services. This indicates that mothers whose husband's educational levels were primary, secondary, and higher levels were more likely to receive ANC services than those not educated. The other studies conducted in Ethiopia also verified that women with partners having primary or higher education used ANC services more than those women who had partners/husbands with no education in Ethiopia [23].
The other factor of interest for this study was the household wealth index. It was found that the household wealth index was an important determinant factor on the usage of ANC service. This indicates that women who are from a household with medium and higher wealth quintile are more likely to utilize ANC service care than those who are from poor wealth households. This positive relationship between wealth index and antenatal care usage was also supported by a study conducted by various researchers [18,24].
The birth order was the other determinant factor that influences the usage of antenatal care services. The study conducted by Kamau and Habtom [25,26] also confirmed this truth. Moreover, research conducted by Habtom [26] showed that the mothers used more antenatal care during the second-order birth compared to the first-order birth but the use of antenatal care decreased with birth order more than two. This research also depicted that the higher the birth order, the lesser the usage of antenatal care.
The husband's occupation was a significant determinant variable on the usage of ANC services. This indicates that women whose husbands are working were more likely to receive ANC services than women whose husbands are not working. This finding is in agreement with the results of the previous study by Terfasa et al. [27] stating that women with husbands in working status were more likely to receive antenatal service care.
The study conducted by Abosse et al. concluded that the attitude towards beating wives that argues with husband was a significant determinant factor on receiving ANC services. This study also approved that women who believed that the 9 Computational and Mathematical Methods in Medicine husband is justified in hitting or beating his wife are less likely to receive antenatal care services. Thus, women who have positive attitudes towards beating wives that argues with husband use less antenatal care services that is also confirmed by Abosse et al. [28].

Conclusion
The study identified some socioeconomic and demographic factors that affect the utilization of antenatal service cares. Based on the findings, the upgrading of women and their husband's education is a very important factor on the utilization of antenatal care services. In this regard, all stakeholders particularly the government should act to raise women's and their husband's education level. Moreover, the government and other stakeholders working on health and related issues should strengthen the effort to improve the accessibility of health facilities, constructing roads and providing transportation services for pregnant women in rural areas. Finally, health extension workers should be given better training to identify and transfer high-risk women for better antenatal services in a health institution.

AIC:
Akakie information criterion ANC: Antenatal care BIC: Bayesian information criterion CSA: Central Statistical Agency DC: Delivery care DIC: Deviance information criterion EDHS: Ethiopian Demographic and Health Survey ICC: Interclass correlation MCMC: Markov chain Monte-Carlo PNC: Postnatal care SPSS: Statistical Package for Social Sciences WHO: World Health Organization.

Data Availability
In this study, we used the information from the fourth Demographic and Health Survey conducted of Ethiopia from January 18, 2016, to June 27, 2016 by Central Statistical Agency (CSA) focusing on all women who are in the reproductive age group (aged from 15-49 years).

Conflicts of Interest
The author declares no conflict of interest