Modeling timing of sexual debut among women in Zimbabwe using a Geoadditive Discrete-Time survival approach

Abstract

Background

Early sexual debut has undesirable health consequences such as an increased risk of contracting sexually transmitted infections (STIs) including HIV, mental health problems, pregnancy-related complications and death including abortion-related deaths. Despite a global decline in adolescent birth rates, Zimbabwe continues to face a high prevalence of underage pregnancies, highlighting significant early sexual debut among Zimbabwean adolescents. This study examined the spatial variation and the demographic and socio-economic determinants of the timing of early sexual debut among Zimbabwean women.

Methods

Data for 9,882 Zimbabwean women of reproductive age were drawn from the 2015 Zimbabwe Demographic and Health Survey (ZDHS). We defined early sexual debut as having first sexual intercourse before the 18 years of age. A fully Bayesian geoadditive discrete-time survival model was used. Adjustments for unequal sampling probabilities were done using the provided survey weights.

Results

Our findings show that women with primary education (aOR = 0.62,95% Crl:0.47–0.81), secondary education (aOR = 0.25,95% Crl:0.19–0.33) and higher education (aOR = 0.06,95% Crl:0.04–0.09) had lower odds of early sexual debut than those with no education. In comparison to those with middle household wealth index, women with higher household wealth index (aOR = 0.83, 95% CrI: 0.71–0.98) had lower odds of early sexual initiation. On the other hand, women with lower household wealth index had higher likelihood of early sexual debut (aOR = 1.13, 95% CrI: 1.03–1.26) than those with middle household wealth index. The type of place of residence and birth year cohort did not have a significant association with the odds of early sexual debut. The hotspots of early sexual debut were in Matabeleland North and Matabeleland South provinces.

Conclusion

To mitigate early sexual debut in Zimbabwe, targeted interventions are essential in Matabeleland North and Matabeleland South provinces as well as in the identified high-risk groups.

Introduction

Early sexual debut remains a concern in public health because it is related to a wide range of negative health outcomes in adolescence, which can last into young adulthood [1] or even the whole lifetime. Research has shown that early sexual debut is often illegal due to adolescent women not having the legal capacity to consent and may be associated with unequal power dynamics as partners tend to be older by at least five years [2, 3]. Unsafe early sexual intercourse has been linked with increased risk of contracting sexually transmitted infections (STIs) including HIV/AIDS, adolescent pregnancies which lead to unsafe abortion or abortion-related death [2, 4, 5], pregnancy-related maternal deaths, and other pregnancy complications such as preeclampsia, anemia, and obstetric fistula [37, 38]. Adolescent pregnancies are also associated with poor neonatal outcomes, including low birth weight, lower Apgar scores, and a higher prevalence of preterm deliveries [37, 38]. These pregnancies often lead to school dropouts resulting in lower levels of literacy and poor school performance [6, 7], substance use, depression, suicidal ideation, reduced levels of life satisfaction, and suicide attempts [4, 6, 8]. In the social sphere, early sexual debut tends to be associated with increased vulnerability to intimate partner violence, risky sexual behaviour and a higher number of lifetime sexual partners [2, 4, 5].

Each year, about 21 million girls aged 15–19 in developing countries become pregnant [9]. While the global adolescent birth rate (ABR) has dropped from 64.5 births per 1,000 women in 2000 to 42.5 in 2021 [9], it remains high in Sub-Saharan Africa at 101 births per 1,000 women. For Zimbabwe, the 2022 ABR was reported to be 93 births per 1,000 women aged 15–19 years [13]. This rate is lower than the rates in Zambia (116 births per 1000 women), Malawi (117 births per 1000 women) and Mozambique (165 births per 1000 women ). Although this rate has been declining since the 1960s and is lower than the rates in the aforementioned neighbouring countries, it is still higher than in South Africa (61 births per 1000 women as of 2022), Botswana (48 births per 1000 women), and Namibia (63 births per 1000 women), and significantly exceeds the global average (see Fig. 1). The high prevalence of adolescent pregnancy in Zimbabwe serves as an indicator of early sexual initiation among young women, underscoring the urgent need for targeted interventions to delay the age of sexual debut and reduce early pregnancies.

Fig. 1

Adolescent Fertility Rates (Number of Births per 1,000 Women, aged 15–19 years) from 1960 to 2022 for Zimbabwe, Namibia, Zambia, Malawi, Botswana, Mozambique, South Africa, Sub-Saharan Africa, and the World.

(Source: United Nations Population Division, World Population Prospects [13]).

Age of sexual consent, that is the minimum age at which a person is considered to have the legal capacity to consent to sexual intercourse, varies across countries globally, ranging from 12 years in Angola to 21 years in Bahrain [10, 11]. In sub-Saharan Africa, the legal age for marriage is set at 18 years in many countries but there is variation in the age of sexual consent ranging from 12 years to 18 years [11]. In Zimbabwe, the age of sexual consent had historically been set at 16 years, but it was subsequently raised to 18 years in 2022 [12]. Currently, Zimbabwean law regards sex with a person below the age of 18 years as statutory rape [34].

To understand and support the development of early sexual debut preventive strategies, it is crucial to understand the factors influencing the timing of sexual debut for adolescents in Zimbabwe. Conducting such a study may lead to identifying key determinants to facilitate effective interventions and evidence-based policies aimed at promoting healthy, consented sexual behaviour. There are a few studies that have focussed on the determinants of early sex in Zimbabwe [14, 15]. Some of the identified risk factors for early sexual debut from previous studies include living in rural areas with lower exposure to media, larger family size and lower socioeconomic status, which are often associated with poverty and family pressure [1, 2, 6, 8]. None of these studies has investigated the possibility of spatial variation i.e. differences across different geographic locations in the timing of first sexual intercourse in Zimbabwe.

The age of sexual debut in Zimbabwe likely differs by location due to differences in cultural norms and other contextual factors. Since policy and governance are regionally implemented, spatial analysis is important for targeted interventions to geographically align with the pre-existing infrastructure. This study, therefore, employs a geoadditive discrete time-survival model to identify districts with high rates of early sexual debut and factors associated with early sexual debut among Zimbabwean women.

A discrete survival approach is chosen because, in the DHS surveys, age at first sex is reported in complete years. In comparison to continuous survival models, the discrete survival model offers several advantages. Hazards for discrete time-to-event models can be formulated as conditional probabilities and are therefore more easily interpretable. Models for discrete time-to-event data do not cause problems with ties and they can be embedded into the family of generalized linear models, making estimation easier [21, 27].

The geoadditive discrete-time survival models used in this study have previously been used in other contexts. For instance, they were used in the analysis of child mortality in Malawi [16] and Nigeria [17] as well as to model age of sexual initiation in Nigeria [18]. Nguyen and Eaton (2022) [19] carried out a study focusing on trends and country-level variation in age at first sex in 42 African countries including Zimbabwe. The spatial unit of analysis in this study was the country therefore the observed spatial variation is likely to mask a lot of variation among smaller spatial units such as districts. This study leverages district-level spatial resolution to shed light on localized variations in early sexual debut, thereby informing targeted interventions aimed at addressing this issue at the local level.

Materials and methods

Data

The data used in this study were obtained from the Zimbabwe Demographic and Health Survey (ZDHS) for 2015 [20]. The survey was implemented by the Zimbabwe National Statistics Agency (ZIMSTAT) from July to December 2015.

Administratively, Zimbabwe is divided into 60 districts spread in ten provinces. Each province in Zimbabwe is divided into districts, and each district is divided into wards. In the 2012 Zimbabwe Population Census, every ward was partitioned into smaller units referred to as census enumeration areas (EAs). For the 2015 ZDHS sample selection, a stratified, two-stage cluster design was employed, using EAs as the sampling units for the initial stage and households as the sampling units for the subsequent stage.

Four hundred EAs, of which 166 were in urban areas and 234 in rural areas were selected and from these, 11,196 households were selected. This study drew a sample of 9,882 women aged 15–49 years from the main data. For each woman, the outcome of interest was the reported age at first sexual intercourse. All women who reported to have had a sexual debut before the age of 18 were regarded as having experienced the event of interest at the reported age at first sex. Women aged 15–17 years at the time of interview who reported to have never had sex were censored at their current age. Those who were aged $\:18-49$ years at the time of the interview and reported to have never had sex were censored at 18 years while those who were in the same age group at the time of the interview and reported having had sex at 18 years or later were also censored at 18 years.

Background characteristics such as place of residence (urban or rural), age, level of education, wealth status, and birth cohort were considered covariates in this study. Place of residence was determined based on the geographic classification provided in the DHS dataset, which categorizes respondents’ residences as either urban or rural. Age was recorded in completed years at the time of the survey. The respondent’s level of education was measured based on the highest level of education attained, typically categorized in the DHS dataset (e.g., no formal education, primary, secondary, higher education). The wealth level was assessed using a wealth index constructed from household asset data, with respondents classified into wealth quintiles based on a composite score from ownership of various assets. Lastly, birth cohorts were categorized according to the year of birth of the respondents, allowing for the analysis of trends across different generational groups.

Data Availability

The data analysed in this study are publicly available at https://dhsprogram.com/.

Statistical analysis

Frequencies and proportions were obtained using STATA software version 17, while adjusting analysis using sampling weights to account for stratification and clustering in the survey design using the “svyset” command. To estimate the model parameters, we use a fully Bayesian framework in BayesX (RBayesX) [22,23,24].

The geoadditive bayesian discrete-time survival model

The geoadditive Bayesian discrete-time survival model specification is as follows:

Let $\:{[a}_{0},{a}_{1}),\:\left[{a}_{1},{a}_{2}\right),\dots\:,,\:{[a}_{q-1}$,$\:{a}_{\infty\:})$ be yearly intervals with $\:{a}_{0}=0$ being the time of birth. Let $\:T\in\:\{1,2,\dots\:,q=18\}$ represent the follow up time in years with $\:T=t$ representing an observation made in the year interval [$\:{a}_{t-1}$,$\:{a}_{t}$). For a woman with covariate vector $\:{\varvec{X}}_{i}=({X}_{i1},{X}_{i2},\dots\:,{X}_{ip,\:})$ the hazard of early sexual debut at $\:T=t$ is given by

$$\:h\left(t|{\varvec{X}}_{\varvec{i}}\right)=P\left({T}_{i}=t|{T}_{i}\ge\:t,{X}_{i}\right)$$

(1)

which is the conditional probability of early sexual debut in the year interval [$\:{a}_{t-1}$,$\:{a}_{t}$) given that the woman has reached the interval without experiencing the event. The discrete survivor function for surviving beyond the year interval [$\:{a}_{t-1}$,$\:{a}_{t}$) is given by:

$$\:S\left(t|{\varvec{X}}_{\varvec{i}}\right)=P\left({T}_{i}>t|{\varvec{X}}_{i}\right)=\prod\:_{k=1}^{t}\left(1-h\left(k|{\varvec{X}}_{\varvec{i}}\right)\right)$$

(2)

The survival information for each woman is recorded as ($\:{t}_{i},{\delta}_{i}$), where $\:i=\text{1,2},\dots\:,9882$. Here, $\:{t}_{i}$ represents the follow-up period. For women who had their first sexual encounter before the age of 18, $\:{t}_{i}$​ is the reported age at sexual debut. For the rest of the women, $\:{t}_{i\:}$is the minimum of the age at the time of the survey and 18 years. A value of $\:{\delta}_{i}=1\:$denotes a woman who had early sexual debut at $\:{t}_{i}$ while $\:{\delta\:}_{i}=0$ denotes a woman who was censored at time $\:{t}_{i}$.

From the pair $\:({t}_{i},{\delta}_{i}$), we create a binary event indicator $\:{y}_{it}$$\:\left\{i=\text{1,2},\dots\:,9882,t=\text{1,2},\dots\:,{t}_{i}\right\}$ as follows:

$$\:{y}_{it}=\left\{\begin{array}{cc}1,&\:if\hspace{0.25em}t={t}_{i}\hspace{0.25em}and\hspace{0.25em}{\delta}_{i}=1\\\:0,&\:\text{elsewhere.}\end{array}\right.$$

In this way, the experience of an individual at different time points is represented as (0,0,….,0,1) representing sexual debut at time t ($\:{y}_{it}=1)\:$or (0,0,…,0,0) denoting remaining without experiencing the event beyond time t $\:\left({y}_{it}=0\right).$ The hazard of the $\:{i}^{th}$ woman having a sexual debut at time $\:t$ can then be expressed as a function of the covariate vector $\:{X}_{i}$, through an equation of the form

$$\:h\left(t|{X}_{i}\right)=P\left({y}_{it}=1|{X}_{i}\right)=F\left({g}_{0}\left(t\right)+{X}_{i}^{T}\varvec{\beta\:}\right)\:\:\:\:\:\:\:\:\:\:\:\:\left(4\right)$$

where $\:F\:\:$ is the (inverse) link function, $\:{g}_{0}\left(t\right)$ is the baseline effect which can be modelled using a smooth function and $\:\varvec{\upbeta\:}$ is a vector of fixed effects. Commonly used link functions are the probit, logit and cloglog link functions. In this study, the logit link,

$$\:F\left({g}_{0}\left(t\right)+{X}_{i}^{T}\varvec{\beta\:}\right)=\frac{{e}^{{g}_{0}\left(t\right)+{X}_{i}^{T}\varvec{\beta\:}}}{1+{e}^{{g}_{0}\left(t\right)+{X}_{i}^{T}\varvec{\beta\:}}}$$

was chosen because it offers convenient interpretation of covariate effects in terms of odds ratios. The covariates considered in the model were the type of place of residence (rural or urban), wealth level of the family, religion, education, and birth year cohort of the woman.

The systematic component, $\:{g}_{0}\left(t\right)+{X}_{i}^{T}\varvec{\beta\:}$, may be generalized to a geoadditive predictor,

$$\:{g}_{0}\left(t\right)+{\varvec{X}}_{\varvec{i}}\varvec{\upbeta\:}+{g}_{spat}\left({s}_{i}\right)$$

where the effects $\:{g}_{0}\left(t\right)$ and $\:{X}_{i}$ are as defined before and $\:{g}_{spat}\left({s}_{i}\right)$ is the non-linear effect of district $\:{s}_{i}\in\:\{1,2,\dots\:,S\}$, where woman $\:i$ lives. The spatial non-linear effects $\:\left({g}_{spat}\left({s}_{i}\right)\:\right)$ are further divided into two components: the spatially correlated or structured component ($\:{g}_{str}\left({s}_{i}\right))$ and the uncorrelated or unstructured effects ($\:{g}_{unstr}\left({s}_{i}\right)).$ There were 60 district spatial units considered in this study.

For fixed effect parameters, $\:\varvec{\upbeta\:}$ diffuse priors are assumed, and for the baseline effect $\:{g}_{0}\left(t\right)$, we assumed Bayesian P-Spline and second order random walk penalty smoothing priors. For the structured spatial effects $\:{g}_{spat}\left(s\right)$, we choose a Gaussian Markov random field prior

$$\eqalign{& \>{g_{spat}}\left( s \right)|\{ {g_{spat}}\left( l \right);s \ne \>l\}, \cr & \sigma {\>^2} \sim \>N\left( {\sum {{\>_{{\rm{l}} \in \>{\rm{Ds}}}}} {g_{str}}\left( l \right)/{N_s},\sigma {\>^2}/{N_s}} \right) \cr} $$

(6)

where $\:{N}_{s}$ is the number of adjacent districts and $\:{D}_{s}$ is the set of districts neighbouring s. Thus, the conditional mean of $\:{g}_{str}\left(s\right)$ is an average of function evaluations $\:{g}_{str}\left(l\right)$ of neighbouring districts $\:l$. The variance $\:{\sigma\:}^{2}$ controls the amount of spatial smoothness. The unstructured spatial effects are independent and identically distributed random effects

$$\:{g}_{unstr}\left(s\right)\sim\:N\left(0,{\varphi\:}^{2}\right).$$

(7)

For all variance components, an inverse gamma distribution with hyperparameters $\:a$ and $\:b$ is assumed, that is, $\:{\varphi\:}^{2}\sim\:IG\left(a,b\right)$.

Fully Bayesian inference is based on the posterior distribution of the model parameters, hence MCMC sampling from full conditionals for non-linear effects, spatial effects, fixed effects, and smoothing parameters is used for posterior analysis. Deviance Information Criteria were used to compare four different specifications of Eq. 4: the no spatial effects model, unstructured spatial effects model, structured spatial effects model and the full spatial model. The models were fitted based on the following specification: 10,000 iterations, a thinning effect of 10, burn-in of 2000 and various priors were tested to assess the sensitivity of the model estimates to the choice of prior. The model diagnostics were done using the trace plot traversed rapidly and the autocorrelation plots, which had a steep fall towards zero, indicated good performance of the model.

Lastly, the prevalence of early sexual debut was mapped for each enumeration area for the 400 clusters in the DHS survey. To allow estimation of the aerial polygon, kriging was used to smooth out the prevalence estimates. The smoothed map was done using the geoprocessing tool in ArcMap. We performed kriging using cluster prevalences and overlaid province and district boundaries using ArcGIS.

Results

Descriptive statistics

Table 1 shows the descriptive information for each of the covariates. Each of the 10 provinces of Zimbabwe contributed roughly 10% of the sample used in this study. The prevalence of early sexual debut appears to be higher on the Western and Northern parts of the country. Matabeleland North province appears to have the highest prevalence. The most prevalent districts are Hwange, Binga, Lupane and Tsholotsho in Matabeleland North province, Bulilima in Matabeleland South province and also Guruve, Centenary and Rushinga in Mashonaland Central province. The prevalence of early sexual debut is lower in the South Eastern parts of the country (Fig. 2).

Fig. 2

Smoothed early sexual debut prevalence map for adolescent women in Zimbabwe.

(Source: Author’s own calculation from 2015 ZDHS data).

The percentage of women who had sex by the age of 18 years was higher in rural as compared to urban areas and decreased with education, wealth index and age. The study also indicated possible variations based on religious affiliations. Participants identifying as Traditional/Muslim/None religion showed a higher percentage of ESI (Early Sexual Debut) at 60%, compared to those from Apostolic Faith religion at 45% and those from other Christian churches at 26% (Table 1).

Table 1 Descriptive statistics of covariates used in the study

Figure 3 shows the empirical baseline hazard function for sexual debut. The hazard of first sex rises steeply between the ages of 14 and 18 years then peaks at around the age of 22 years and begins to drop from 25 years onwards. The smallest reported age at sexual debut was 8 years. The median age at sexual debut was 17.74 years. This was obtained by interpolating the lifetable survival probabilities to determine the time at which the cumulative survival probability is closest to 0.5. In what follows, all women who experienced sexual debut after the age of 18 years were censored at 18 years.

Fig. 3

Empirical baseline hazard function

Geo-spatial survival model results

Four models (no spatial model, unstructured spatial effects, structured spatial effects and the full spatial model which comprises of both the structured and unstructured components) were fitted to identify the factors associated with the timing of the early sexual debut (Table 2). The model fit was assessed using the effective numbers of parameters (PD) and the Deviance Information Criterion (DIC). The model with the least DIC is the full spatial model (DIC = 21748.28). The discussion that follows is therefore based on this model.

Table 2 Geo-additive survival model results

Fixed effects

Table 2 shows the posterior means of the adjusted odds ratios together with the 95% credible intervals (CrI). Increased levels of education were associated with lower odds of early sexual debut. For a woman with primary education, the odds of reporting early sexual debut were about 62% (95% CrI: 47–81%) of the odds for a woman with no education. On the other hand, for a woman who attained higher education, the odds of reporting early sexual debut are only about 6% (95% CrI: 4–9%) of the odds for a woman with no education.

Higher family wealth index was also associated with lower odds of early sexual debut. For a woman from a rich family, the odds of reporting an early sexual debut were 0.83 (95% CrI: 0.71–0.98) times those for a woman with medium family wealth status whereas for a woman from a poor family, the odds of early sexual debut were 1.1 3 (95% CrI: 1.03–1.26) times those for a woman with medium family wealth status. The odds of early sexual debut do not differ significantly across year cohorts. For women affiliated with the Muslim/Traditional/No religion, the odds of early sexual debut were about one and half times those for women affiliated to the Apostolic Faith churches while for women affiliated to the other Christian churches, the odds of early sexual debut were 0.74 times those for women affiliated to the Apostolic Faith churches. Women from rural and urban areas had the same odds of early sexual debut.

Spatial effects

The spatial effects are displayed in Fig. 4 for the unstructured (left) structured (centre) and total (right) spatial effects on the scale of the log odds of sexual debut. The random effects variance components in our model quantify the district-level variability in the risk of early sexual debut that is not explained by the fixed effects. Specifically, the structured random effects account for spatially correlated variability between districts (i.e., neighbouring districts that share similar risk profiles), while the unstructured random effects capture spatially uncorrelated variability, reflecting any residual, district-specific risk that cannot be explained by either the fixed effects or the structured spatial patterns.

The structured effects were more pronounced, contributing 79% to all spatial variability, indicating that early sexual debut may be driven by district specific characteristics. Figure 4 below also highlights the dominance of structured spatial effects on the risk of early sexual debut. This dominance suggests that neighbouring districts exhibit similar risks of early sexual debut, potentially due to shared environmental, cultural, or socio-economic factors that affect sexual behaviour patterns across these districts.

The structured random spatial map shows a West-East divide with districts towards the West (especially the South West such as Bulilima and Mangwe) having a higher risk of sexual debut in comparison to districts located towards the eastern part of the country. The unstructured spatial map shows that districts to the South West of the Country and some districts in the Northern Central part of the country such as UMP, Kadoma and Chikomba districts have a relatively higher risk of sexual debut in comparison to the rest of the country. This is the same pattern observed in the total spatial effect. In general, early sexual debut was prominent in Matabeleland North province which borders Botswana and Zambia.

Fig. 4

Unstructured spatial effects(left), Structured spatial effects (center) and Total spatial effects(right) highlight distribution of spatial risk for early sexual debut

Discussion of the results and implications of the findings

This study used a Bayesian geoadditive discrete survival model to estimate geographical variations of the time to early sexual debut among women in Zimbabwe using the 2015 DHS data. We explored various factors including birth year cohort, household wealth index, education, religion, type of place of residence, and their associated correlation with early sexual debut for this population.

After adjusting for spatial effects, household wealth index, level of education and religion were found to be significant. Although not statistically significant, our results indicate that the odds of early sexual debut among rural women were 15% higher than those for women in urban areas. This finding aligns with previous studies that reported living in rural areas as a risk factor for early sexual debut [1, 2, 6, 8, 25].

Our results show that the odds of early sexual debut do not significantly vary by year cohort. On the other hand, Nguyen and Eaton (2022) [19] reported that in Zimbabwe the trend was towards older age at sexual debut in recent generations. However, their results were obtained after correcting for reporting bias whereby younger women tend to report a later age at sexual debut than when interviewed at later ages [19, 26]. It is possible that the lack of bias correction in our work may explain the difference in findings on cohort effects. Our study is based on cross sectional data; hence it was not possible to correct for reporting bias as Nguyen and Eaton did. It is also possible that the differences in results from the two studies could be emanating from the fact that the models were adjusted for different covariates.

We also found that the odds of early sexual debut decrease as education levels rise. This finding on the positive effect of education could be interpreted alongside other studies conducted in Malawi [27] and the USA [28] which reported lower odds of engaging in early sexual debut among more educated women. The life course theory [35] provides a useful framework for interpreting this result. This theory emphasizes that the timing and sequencing of life events (such as the onset of sexual debut) are influenced by various factors, including age, social norms, historical period, and cohort effects. For instance, the association between higher levels of education and delayed sexual debut can be explained through the Life Course Theory. Women with more education may have better opportunities, higher aspirations, and more resources to delay engagement in risky behaviours. However, interpreting the effect of education requires caution due to the bidirectional relationship between early sexual debut and educational attainment. Early sexual debut can lead to teenage pregnancy and subsequent school dropout, indicating that our findings should be understood as associations rather than direct causal effects. Further research using longitudinal data is necessary to clarify these relationships. Longitudinal data would allow for the incorporation of educational attainment as a time-varying covariate, providing a clearer understanding of the links between educational attainment and early sexual debut. However, due to the cross-sectional nature of our dataset, educational attainment was treated as a fixed covariate, which limits the ability to make causal interpretations.

Our findings on the association between wealth and early sexual debut are also consistent with the Life Course Theory. Women from wealthier families have lower odds of early sexual debut, which may reflect the greater access to education, healthcare, and parental supervision that often accompanies higher socio-economic status. Similar to education, interpreting the effect of wealth requires caution due to the bidirectional relationship between socio-economic status and early sexual debut; early sexual debut can lead to outcomes such as teenage pregnancy and lower economic stability. Additionally, it is possible for a woman from a poor background to marry into an affluent family. The true nature of the relationship between these variables can only be clarified through longitudinal data. Therefore, future research should prioritize longitudinal studies to better understand the dynamic interplay between wealth, education, and early sexual debut.

A common belief in Zimbabwe associates adolescent marriages and pregnancies with membership of the Apostolic Faith churches [29, 30]. This study’s results indicate that, despite higher risks compared to women from other Christian denominations, Apostolic Faith church members have significantly lower odds of early sexual debut compared to women from Muslim or non-religious backgrounds. These findings are consistent with previous studies suggesting that religion is a significant predictor of sexual abstinence [31] and delayed initiation of sexual intercourse could be explained by religious affiliation [32]. It must however be noted that combining the Muslim and non-religious groups was done, not because they had similar outcomes, but because each group had very few individuals, which could have made the parameter estimates less stable. It is therefore possible that if the two groups were to be separated, the picture in terms of how each group compares to the other religions may be different from what is portrayed in the current results.

The most significant finding of this study is the spatial variation of the risk of early sexual debut in the country. Women from districts in Matabeleland North and Matabeleland South have the highest risk of early sexual debut. It is possible that the disparity between these two provinces and the rest of the country is due to cultural differences. The Ecological Systems Theory [36] offers a useful lens through which this result could be interpreted. This theory posits that individual behaviour is influenced by multiple layers of environmental factors, including personal, family, community, and societal contexts. The West-East divide observed in our spatial results suggests that district-level characteristics, including socio-economic and cultural differences, significantly influence sexual behaviours. For example, districts in Matabeleland are characterized by specific socio-cultural practices that may contribute to earlier sexual debut. The structured spatial effects, which contributed 79% of the spatial variability, indicate that neighbouring districts share similar risk profiles, possibly due to shared cultural norms or environmental factors.

It is also interesting to note the two provinces with the highest odds of sexual debut were also the provinces with the highest prevalence of HIV [33]. It is worth noting that the province of Bulawayo shows a lower risk of early sexual debut, which is surprising considering its location within the Matabeleland region. Being a metropolitan province Bulawayo is inhabited by people of different cultures which could explain the disparity with the rest of Matabeleland provinces. Women from the Midlands province and some parts of Mashonaland central, East, and West provinces also showed a fairly high odds of early sexual debut. These results are consistent with the findings of the Zimbabwe Multiple Indicator Cluster Survey 2019 [25] which reported that Mashonaland Central and Matabeleland South provinces were among the provinces with the highest rates of early sexual debuts ranging between 6.2 to 6.8%.The strength of this study is that it leverages district-level spatially correlated survival data to gain deeper insights into localized variations in and determinants of early sexual debut. However, its limitations are that the retrospective data used, by nature, can be susceptible to inaccuracies due recall bias. Additionally, as reported by Nguyen and Eaton (2022) [19] and Zaba et al. (2004) [26], the sensitivity surrounding this topic may lead participants to provide responses that they deem socially acceptable.

Conclusion

This study used a geoadditive discrete-time survival model estimated in a Bayesian framework to investigate the geographical variations in the timing of early sexual debut for women in Zimbabwe and to evaluate the potential influence of socio-demographic factors on the timing of these early first sexual encounters. The findings indicate that higher levels of education and wealth status correlate with decreased odds of early sexual debut, and women from the Christian religions face lower odds compared to other religious groups. Notably, spatial analysis revealed that the Western part of the country and some central and northeastern districts exhibit higher odds of early sexual debut. We recommend that any interventions may be targeted towards these high-risk areas.