Spatial Analysis of the Determinants of HIV/AIDS Among Females Aged 15–34 in KwaZulu Natal, South Africa Using a Bayesian Spatial Logistic Regression Model

1. Background

Human immunodeficiency virus (HIV) remains a significant public health concern in South Africa, home to the largest population of individuals living with HIV globally. KwaZulu-Natal (KZN) is the most affected province, with the highest prevalence rates and profound socio-economic impacts [1]. Women aged 15–34 in KZN are a key demographic in the fight against HIV/AIDS, accounting for a substantial proportion of new infections due to biological vulnerability, societal pressures, and economic challenges [2,3]. Literature documents that women in this age group face heightened risks of HIV transmission, compounded by factors such as inconsistent condom use, early sexual debut, substance abuse, and concurrent sexually transmitted infections (STIs) [4,5]. Poverty, unemployment, transactional sex, and intergenerational relationships further contribute to their vulnerability [6]. In addition, intimate partner violence (IPV) exacerbates the risk by limiting women's ability to negotiate safer sexual practices [7]. Geospatial disparities in HIV prevalence in KZN, particularly in peri-urban and rural areas, highlight barriers such as limited healthcare access and high poverty rates [8]. Spatial epidemiology has identified clusters of high prevalence, underscoring the need for targeted interventions [9].

Bayesian spatial logistic regression offers a robust framework for examining spatial, demographic, and individual-level factors influencing HIV prevalence. This approach incorporates spatial dependencies and heterogeneity, essential for understanding geographic variability in HIV burden across KZN [10]. The application of Bayesian spatial logistic regression has been demonstrated in several contexts. For instance, [11] used Bayesian semi-parametric regression to analyse HIV prevalence among men in Kenya, integrating structured and unstructured spatial effects. Similarly, [12] applied Bayesian spatial modelling to study tuberculosis-HIV co-infection in Ethiopia, uncovering significant geographical heterogeneity. These studies underscore the utility of Bayesian approaches in informing public health interventions through spatially detailed risk analyses.

Structured additive models (SAMs) further enhance Bayesian approaches by accommodating non-linear effects of continuous variables alongside spatial random effects [13,14]. This flexibility allows for integrating individual and area-level risk factors, enabling more accurate spatial estimates and better identification of high-risk clusters [15]. The structured additive approach has proven valuable in addressing spatial autocorrelation while elucidating HIV drivers in regions like KZN, where complex interactions between socio-demographic and geographic factors exist. Spatial models like the Besag-York-Mollié (BYM) model have been widely used in HIV research to analyse structured and unstructured spatial effects, often implemented efficiently using integrated nested Laplace approximation (INLA) [16,17]. Such models have revealed significant geographic disparities in HIV prevalence, aiding in the identification of high-risk clusters and guiding public health interventions [15].

The study aimed to investigate the spatial distribution and key demographic, behavioural, and socio-economic factors associated with HIV prevalence among female youth in KwaZulu-Natal, using a Bayesian spatial logistic regression framework with a structured additive model. To the best of our knowledge, there are few studies integrating advanced Bayesian spatial logistic regression framework to exclusively capture both micro-level (individual risk factors) and macro-level (spatial dependencies) determinants of HIV prevalence targeting female youth in KwaZulu Natal. As a result, this study was the first of its kind.

Youth is often defined as individuals aged 15–24 (United Nations). However, this study adopts a broader definition, encompassing individuals aged 15–34, in alignment with regional demographic trends and epidemiological significance. This age range captures critical life transitions that influence HIV risk, including adolescence, early adulthood, and early middle age, as highlighted by [18]. The expanded age range also reflects the South African context, where young adults up to the age of 34 face significant socioeconomic and health vulnerabilities.

2. Methodology

2.1. Study Area Location

Figure 1A and 1B below depict the location of the study area within uMgungundlovu District and the location of the two sub-districts of KwaZulu-Natal Province, namely: Vulindlela (Western part) and the Greater Edendale (Eastern part), respectively.

2.2. Sources of Data and Study Population

We conducted a secondary analysis of data collected from participants in the HIV Incidence Provincial Surveillance System (HIPSS). Briefly, HIPSS was undertaken from 11 June 2014 to 18 July 2015 in rural Vulindela and the peri-urban Greater Edendale areas in uMgungundlovu District (Figure 1A and B) of KZN, South Africa.

From a total of 600 enumeration areas, 591 enumeration areas with more than 50 households were included in the sample. Of these, 221 enumeration areas were drawn randomly. Within an enumeration area, the households were drawn systematically with a random start. The study staff identified households using the Global Positioning Systems receiver to record the geographic coordinates of each randomly selected household. Only one age-eligible individual per household was randomly selected and enrolled following written informed consent. Questionnaires were administered to obtain household- and individual-level data on demographics, socio-economic status, and health-related information. All enrolled participants provided peripheral blood samples for laboratory measurements for HIV.

In the South African context, using the age group 15–34 as a definition of "youth" aligns with official policy frameworks [18]. This range is designed to account for the extended transitional period many young people experience in South Africa due to socioeconomic factors such as unemployment, prolonged education, and delayed family formation​. This age group was chosen to capture both adolescent and young adult populations, reflecting the epidemiological realities of HIV risk in KwaZulu-Natal and aligning with local public health policies. A total of 3324 female participants aged between 15–34 years were involved in our study.

2.3. Study Variables

The dependent variable in this study was “HIV prevalence” which was defined as the ratio of the number of HIV positive participants in an enumeration area to the total number of participants in the same enumeration area. In our analysis, we used unweighted HIV prevalence since we were focused on detecting geographic locations where spatial clustering of HIV prevalence occurs.

HIV status among participants in the study population was categorised as a binary outcome:

$$y_{ij}=\left\{\begin{array}{c}1 for a participant who tested positive\\ 0 for a participant who tested negative\end{array}\right.$$

(1)

The covariates included in the study comprised sociodemographic, behavioural, and biological variables. These included age, level of education, marital status, main income, ever consumed alcohol, ever diagnosed of TB, ever diagnosed of STI, number of sexual partners, condom use, forced first sex, ever pregnant, run out of money, meal cuts, being away from home, duration in the community and access to health facilities. We applied the variance inflation factor (VIF) to check for multicollinearity before fitting the models. All the variance inflation factor values were very small, all less than 1.5, indicating that multicollinearity was not a major issue in the fitted model.

2.4. Spatial Autocorrelation

Spatial analysis was used to identify HIV prevalence clustering and examine the influence of related variables. Enumeration Areas (EAs), with geo-referenced boundaries, were the spatial units linked to HIV prevalence data. Observations that are close in space tend to have similar values and exhibit spatial autocorrelation. Spatial models account for this to distinguish between general trends driven by covariates and spatial random variation [19]. The global Moran’s index and Geary’s C statistic were employed to assess spatial autocorrelation, determining whether HIV prevalence is dispersed, random, or clustered.

2.4.1. Global Moran’s Index Statistic

The Global Moran’s Index measures overall spatial autocorrelation across a study area, indicating the presence, strength, and direction of spatial patterns. Positive autocorrelation occurs when neighboring enumeration areas have similar values, while negative autocorrelation suggests contrasting values. When spatial patterns are random, the index approaches zero [20,21,22,23].

The Moran’s Index is calculated as:

$$I={\displaystyle \frac{n}{\mathit{\Theta }}}\times {\displaystyle \frac{\sum _{i=1}^n{\sum }_{j=1}^n{\mathsf{\Theta }}_{ij}\left(y_i-\overline{y}\right)\left(y_j-\overline{y}\right)}{\sum _{i=1}^n{\left(y_i-\overline{y}\right)}^2}}$$

(2)

where $n$ is the total number of enumeration areas, $y_i$ is the value of the variable at location $i$, $\overline{y}$ is the mean of the variable $y$ across all enumeration areas, ${\mathsf{\Theta }}_{ij }$is the spatial weight between enumeration area $i$ and enumeration area $j,$ and $\mathsf{\Theta }$ indicates the sum of all spatial weights.

2.4.2. Geary’s C Statistic

Geary’s C evaluates spatial autocorrelation by assessing similarity or dissimilarity between values at neighboring locations. Unlike Moran’s Index, it is sensitive to local variations [24,25]. The formula for Geary’s C is:

$$C={\displaystyle \frac{(N-1){\sum }_{i=1}^N\sum _{j=1}^N{\mathit{\Theta }}_{ij}{\left(y_i-y_j\right)}^2}{2\mathit{\Theta }\sum _{\mathrm{i}=1}^{\mathrm{N}}{\left({\mathrm{y}}_{\mathrm{i}}-\overline{\mathrm{y}}\right)}^2}}$$

(3)

where N is the total number of enumeration areas (locations), $y_i$ and $y_j$ are the values of the variable of interest at locations i and j, $\overline{\mathrm{y}}$ is the mean of the variable across all locations, ${\mathit{\Theta }}_{ij}$ is the spatial weight between location i and location j, and $\mathit{\Theta }$ is the sum of all ${\mathit{\Theta }}_{ij}$. Values of C < 1 indicate positive spatial autocorrelation, C > 1 indicate negative spatial autocorrelation, and C = 1 implies no autocorrelation [26].

While Moran’s Index and Geary’s C identify spatial autocorrelation, they cannot differentiate between hot-spots and cold-spots within clusters. For this, Kulldorff’s spatial scan statistic (SaTScan) can be used. It detects significant spatial clusters of risk factors, identifying areas with higher or lower HIV prevalence in circular windows, enabling hotspot and cold-spot analysis [27].

2.5. Bayesian Logistic Regression Models

Bayesian logistic regression is a powerful method for modelling binary outcomes, such as disease presence, by estimating posterior distributions of regression parameters. This approach integrates prior beliefs with observed data, producing posterior distributions that reflect both sources of information [16,28].

The binary outcome $Y_i\in \{0;1\}$ follows a Bernoulli distribution:

$Y_{i }~$ Bernoulli ($p_i)$, where $p_i$ is the probability that $Y_i=1$, linked to the linear predictor ${\phi }_i$ by the logistic function:

$$p_i={\displaystyle \frac{\mathrm{e}\mathrm{x}\mathrm{p}\left({\phi }_i\right)}{1+ exp\left({\phi }_i\right)}}$$

(4)

and

$${\phi }_i={\beta }_0+X_i^T\mathsf{\beta }$$

(5)

where ${\beta }_0$ is the intercept, ${\mathbf{X}}_i$ is the vector of covariates for observation $i$ and $\mathsf{\beta }$ is the vector of regression coefficients. The likelihood function of N observations is expressed as:

$$L\left(\beta \right)=\prod _{i=1}^N{p}_i^{Y_i}{\left(1-p_i\right)}^{1-Y_i}$$

(6)

Bayesian spatial logistic regression extends this framework by incorporating spatial dependencies, enabling the analysis of structured and unstructured spatial variability in binary data such as disease prevalence [29].

The model is given by: $Y_{i }~$ Bernoulli ($p_i)$,

with the probability$p_i$ linked to the linear predictor ${\phi }_i$:

$$p_i={\displaystyle \frac{\mathrm{e}\mathrm{x}\mathrm{p}\left({\phi }_i\right)}{1+exp\left({\phi }_i\right)}}$$

(7)

The linear predictor includes spatial random effects:

$${\phi }_i={\beta }_0+X_i^T\mathit{\beta }+{\mathit{\theta }}_i$$

(8)

where ${\mathit{\theta }}_i$ represents the spatial random effects.

Spatial dependencies are captured using priors like the conditional autoregressive (CAR) model, intrinsic CAR (ICAR) model, or Gaussian process (GP) model, allowing robust modelling of spatially correlated binary outcomes [30,31].

2.6. Prior Distributions

In Bayesian analysis, prior distributions represent beliefs about parameters before observing data and are combined with likelihood functions to obtain posterior distributions [16,31,32]. Priors are essential in hierarchical spatial models, especially with small sample sizes or variable data, and help regularize the model [33].

Choosing priors involves balancing prior knowledge and non-informativeness. Informative priors guide inference when prior knowledge is available, while weakly informative or non-informative priors are used when prior knowledge is absent. In this study, non-informative priors were used for regression coefficients and random effects variances due to lack of prior knowledge.

Penalized complexity (PC) priors were applied to the precision parameter of the random effects. These priors balance model simplicity and complexity, avoiding issues like overfitting and computational problems associated with flat priors [34,35]. The PC prior for precision $\rho$ is expressed as:

$$\pi \left(\rho \right)=\nu e^{-\nu \sqrt{\rho }}$$

(9)

with

$$\nu ={\displaystyle \frac{log\left(\alpha \right)}{U}}$$

(10)

and

$$\rho ={\displaystyle \frac{1}{{\sigma }^2}}$$

(11)

where $\rho$ is the precision, $U$ is the upper bound for the standard deviation $\sigma$ of the random effect, and $\alpha$ is the probability that $\sigma >U$.

2.7. Posterior Distributions and Point Estimates

The posterior distribution contains complete information about parameter estimates, summarized using point estimates and credible intervals. Point estimates include the posterior mean, posterior mode, and posterior median, which are used for inference and prediction.

The posterior mean is the expected value of the parameter under the posterior distribution. It is a common estimate, especially when the posterior is symmetric. For a parameter $\beta$, it is given by:

$${\widehat{\beta }}_{mean}=\mathbb{E}\left[\beta |y\right]=\int \beta P\left(\beta |y\right)d\beta$$

(12)

The posterior mode also known as the maximum a posteriori (MAP) estimate is the mode of the posterior distribution, i.e., the value of $\beta$ that maximises $P\left(\beta |y\right)$ and is expressed as:

$${\widehat{\beta }}_{MAP}={argmax}_{\beta }P\left(\beta |y\right)$$

(13)

where ${arg max}_{\beta }$ indicates finding the value of $\beta$ that maximises this posterior probability. The MAP estimate is often used when the posterior is skewed, but it can be sensitive to the choice of the prior.

The posterior median is a robust point estimate that divides the posterior distribution into two equal parts. It is less sensitive to outliers compared to the mean or mode.

Credible intervals provide the range where the parameter likely falls with a given probability. A 95% credible interval means there is a 95% probability that the true parameter lies within the interval:

$$P\left({\beta }_{lower}<\beta <{\beta }_{upper}|y\right)=0.95$$

(14)

Unlike frequentist confidence intervals, credible intervals offer direct probabilistic interpretation.

2.8. Bayesian Spatial Logistic Regression Models Applied

Bayesian logistic regression incorporates prior beliefs and spatial dependencies. Below are the applied models.

2.8.1. Unstructured Bayesian Spatial Logistic Regression Model:

This model accounts for heterogeneity by incorporating independent and identically distributed random effects, assuming no spatial dependency [16,36]. It is defined as: $Y_{i }~$ Bernoulli ($p_i)$,

with

$$logit\left(p_i\right)={\beta }_0+X_i^T\mathit{\beta }+u_i$$

(15)

where $u_i$ denotes the unstructured random effects and $u_i~\mathsf{{\rm N}}\left(0,{{\sigma }^2}_u\right)$.

2.8.2. Structured Bayesian Spatial Logistic Regression Model:

This model incorporates spatial dependence using a structured random field, improving predictions by considering the influence of nearby locations. It is defined as:

$$Y_{i }~\mathrm{Bernoulli} \left(p_i\right), \mathrm{with} logit\left(p_i\right)={\beta }_0+X_i^T\mathit{\beta }+{\mathit{\theta }}_{\mathit{i}}$$

(16)

where $\mathit{\theta }~CAR\left(\mathit{W}, {\rm T}\right)$ and the conditional autoregressive (CAR) model for $\mathit{\theta }$ assumes:

$${\theta }_i|{\theta }_{-i},\mathsf{{\rm T}}~\mathsf{{\rm N}}\left({\displaystyle \frac{1}{{\eta }_i}}{\sum }_{j\in neigh\left(i\right)}{\theta }_j,{\displaystyle \frac{1}{{\mathsf{{\rm T}}\mathsf{\eta }}_i}}\right),$$

where ${\mathit{\theta }}_{\mathit{i}}$ indicates the spatially structured random effect at location $i$, ${\eta }_i$ represents the number of neighbours of location $i$, $\mathit{W}$ is a spatial adjacency matrix and $\mathsf{{\rm T}}$ is the precision parameter [13,30,37].

2.9. Model Selection Criteria

After constructing the two Bayesian spatial logistic regression models to capture different spatial structures, we compared their performance. Model selection was based on the following model selection criteria: deviance information criteria (DIC) [38], the effective number of parameters (pD), the mean deviance ($\tilde{D}$) and the Watanable-Akaike information criteria (WAIC) [39]. Lower DIC, $\tilde{D}$ and WAIC values and a higher pD value suggest a better model fit. Hence the best-fitting model was selected based on the smallest DIC, $\tilde{D}$ and WAIC and the highest pD.

2.10. Model Diagnostics

After selecting the best-fitting model, we assessed its adequacy using residual plots and normal Q-Q plots. A well-fitted model should have residuals symmetrically distributed around zero, with no clear pattern or trend and constant variance [40,41,42]. Deviations from normality in the Q-Q plot suggest that residuals do not follow a normal distribution. We also examined spatial autocorrelation in residuals using Moran’s I statistic, Geary’s C statistic, and the variogram plot to verify whether the spatial structure was adequately captured. High spatial autocorrelation in residuals indicates the model failed to fully account for spatial dependencies [43,44]. Significant Moran's I and Geary’s C values suggest poor model fit. Increasing semi variance with distance indicates spatial autocorrelation, suggesting an inadequate model. Flat variogram suggests spatially uncorrelated residuals, indicating a well-fitted model [43,45]. Additionally, posterior density plots were examined for model validity, reliability, and stability. A smooth, unimodal density plot indicates a well-fitting model, while a multimodal plot may suggest model ambiguity or data issues [16].

2.11. Software and Implementation

The Bayesian spatial logistic regression models were implemented using the Integrated Nested Laplace Approximation (INLA) method [14,46] in R (version 4.4.0). The following R packages were used: “INLA”, “sf”, “sp”, “spdep”, and “dplyr” packages. Spatial relationships between enumeration areas were established using a spatial weight matrix, with neighbours identified via Queen’s contiguity. Additionally, Kulldorff’s spatial scan statistics were applied using SaTScan (version 10.1.3).

3. Empirical Results

Summary statistics for the HIV prevalence rates for all the covariates included in the study are depicted in Table 1. While the summary statistics provide an initial indication of associations, the Bayesian model results are prioritized due to their robustness in adjusting for spatial correlations and confounding effects. This approach ensures that our conclusions are based on a more comprehensive analysis of the data.

There were 3324 females who were included in this research and 1576 individuals were HIV positive giving us an overall HIV prevalence of 47.4% (95% CI: 45.7–49.1), (p-value < 0.0001). We noticed that HIV prevalence increased as age increased, and it was 20.4% (95% CI: 16.8–24.5), 37% (95% CI: 34.2–40.0), 54% (95% CI: 50.8–57.1), and 67.5% (95% CI: 64.2–70.8) for age groups 15–19, 20–24, 25–29 and 30–34 respectively, (p-value < 0.0001). Considering education level, individuals with primary education had the highest HIV prevalence of 70.6% (95% CI: 59.7–80.0) followed by those with no schooling with 55.6% (95% CI: 44.1–66.6), (p-value < 0.0001). Participants who had no source of income had the highest HIV prevalence of 50.5% (95% CI: 43.4–57.6), (p-value = 0.169432).

Looking at the marital status covariate, participants who were divorced and those who were separated, but still legally married had the highest HIV prevalence of 100% (95% CI: 15.8–100.0) followed by participants who were single but had been living with someone as a husband/wife before with an HIV prevalence of 63.7% (95% CI: 55.5–71.8). The p-value for the marital status covariate is 0.000181. The HIV prevalence was higher among participants who were once diagnosed with TB, 63% (95% CI: 54.6–70.8) compared to those who were not diagnosed with TB, 46.7% (95% CI: 45.5–48.5), (p-value = 0.000365). Participants who indicated that they were not using condoms as a prevention method had a higher HIV prevalence, 58.1% (95% CI: 47.0–68.7), compared to those who were using condoms, 47.1% (95% CI: 45.4–48.9), (p-value = 0.056253). Classified by number of sexual partners, HIV prevalence increased as the number of partners increased, and it was 45.5% (95% CI: 43.6–47.3) for participants with 1 partner, 51.6% (95% CI: 45.9–57.3) for participants with 2 partners, and 67.5% (95% CI: 60.6–73.8) for participants with 3 partners, (p-value < 0.0001). HIV prevalence for participants who did not consume alcohol was slightly lower, 45.8% (95% CI: 43.9–47.6), compared to those who were consuming alcohol, 58.5% (95% CI: 53.7–63.3), (p-value < 0.0001). Participants who were diagnosed with STIs had a higher HIV prevalence of 63% (95% CI: 56.2–69.4), compared to 46.3% (95% CI: 44.5–48.1) for participants who were not diagnosed with STIs, (p-value < 0.0001). Based on the forced first sex covariate, participants who had forced first sex had the highest HIV prevalence of 54.7% (95% CI: 43.5–65.4), (p-value = 0.246837). Participants who were away from home had a higher HIV prevalence of 50.1% (95% CI: 44.8–55.5), compared to those who were not away from home, (p-value = 0.407053). For the length in community covariate, the highest HIV prevalence of 60% (95% CI: 14.7- 94.7) was recorded for participants who did not respond and the lowest HIV prevalence of 46.7% (95% CI: 44.8–48.7) was observed for those participants who were always in the community, (p-value = 0.447987). The HIV prevalence was higher among participants who accessed health care, 50.5% (95% CI: 47.7–53.4) compared to those who did not respond, 33.3% (95% CI: 4.3–77.7) and those who did not access health care, 45.6% (95% CI: 43.4–47.8), (p-value = 0.018296). Considering run out of money covariate, participants who ran out of money had the highest HIV prevalence of 49.1% (95% CI: 45.2–52.9), (p-value = 0.618173). The HIV prevalence was slightly higher among participants who had meal cuts, 47.7% (95% CI: 43.7–51.8) compared to those who had no meal cuts, 47.5% (95% CI: 45.6–49.4), (p-value = 0.515632). Lastly, participants who once became pregnant had a higher HIV prevalence of 50.4% (95% CI: 48.4–52.3) compared to those who had never became pregnant, 37.4% (95% CI: 33.9–40.9), (p-value < 0.0001).

The HIV prevalence also varied among enumeration areas (ranging between 0–100%). The geographical distribution of HIV prevalence by enumeration areas is shown in Figure 2. This map was created using ArcGIS software with the application of the “tidyverse”, “sf”, and “tmap” packages in R software.

The result for Moran’s index statistic of HIV prevalence was 0,8157 with a p-value < 0.001, indicating a very strong positive spatial autocorrelation in the wards of uMgungundlovu District (Table 2). The positive and statistically significant Moran’s index value supports that there are clusters of high and low HIV prevalence areas within the study region, suggesting a non-random spatial pattern. The positive Moran’s index also suggests that the HIV prevalence in any two spatial neighbouring wards tended to have similar HIV prevalence.

Furthermore, findings from Geary’s C test statistics support the results from Moran’s index statistic as they both reveal consistent evidence of spatial heterogeneity in HIV prevalence within uMgungundlovu District. The summary statistics results for Moran’s index statistic and Geary’s C statistic are displayed in Table 2 below.

As shown in Table 2, both Moran’s I and Geary’s C indicate significant and strong positive spatial autocorrelation in HIV prevalence. These results confirm spatial heterogeneity, suggesting that HIV prevalence is not randomly distributed but influenced by underlying spatial processes or risk factors in uMgungundlovu District. However, while Moran’s I and Geary’s C detect spatial autocorrelation, they do not differentiate hotspots from cold spots. To address this, Kulldorff’s spatial scan statistics were applied, identifying two clusters of HIV prevalence. The spatial distribution of these clusters is visualized in Figure 3.

Cluster 1, a hotspot with a 2.53 km radius, had an HIV prevalence of 48.4%, a relative risk (RR) of 1.22, and a p-value of 0.025, indicating a 22% higher risk inside the cluster compared to outside. This cluster was located around Greater Edendale. Cluster 2, another hotspot with a 2.28 km radius, had an HIV prevalence of 49.6% and an RR of 1.28, meaning the risk was 28% higher within the cluster. This cluster covered Nadi, KwaMbanjwa, Zayeka, KwaMtogotho, KwaNxamalala, and Henley. However, it was not statistically significant (p = 0.467).

To identify factors associated with HIV prevalence, Bayesian spatial logistic regression was applied, considering sociodemographic, behavioural, and biological factors. Most covariates were statistically significant at the 5% level across all three models. Model selection was based on DIC, pD, $\tilde{D}$, and WAIC, as shown in Table 3.

Based on WAIC, DIC, and pD, the structured model emerged as the best model. It has the lower DIC, pD and WAIC values as shown in Table 3 compared to the unstructured model. The structured model strikes the best balance between model fit, complexity, and predictive accuracy, making it the optimal choice. Hence the results of this research are based on the structured model.

Adjusted odds ratios (OR) together with their corresponding 95% credible intervals (CI) for the participants’ characteristics are displayed in Table 4. These values were obtained from the fitted structured Bayesian spatial logistic regression model implemented in INLA.

Most of the covariates included in the study were significant providing insights into the factors associated with HIV prevalence. Covariates levels with 95% credible intervals including 1, were not statistically significant, and as a result, we did not consider them as predictors of HIV prevalence in our study.

The findings revealed that the odds of HIV prevalence for participants in the age groups 20–24, 25–29 and 30–34 were 2.3373 (OR = 2.3373, 95% CI: 1.7914–3.0526), 4.7446 (OR = 4.7446, 95% CI: 3.6111–6.2339) and 9.1981 (OR = 9.1981, 95% CI: 2.8826–12.2926) times higher than that of age group 15–19, respectively.

Considering education, participants with incomplete secondary were 1.4049 (OR = 1.4049, 95% CI: 1.1948–1.6520) times more likely to be HIV infected compared to those with complete secondary. Participants with no schooling were 1.7177 (OR = 1.7177, 95% CI: 1.0650–2.7732) times more likely to be HIV infected compared to participants with complete secondary. Also, participants with primary education were 2.6117 (OR = 2.6117, 95% CI: 1.5968–4.2759) times more likely to be HIV infected compared to participants with complete secondary. Importantly, participants with tertiary education were 0.5337 (OR = 0.5337, 95% CI: 0.3910–0.7276) times less likely to be HIV infected compared to those with complete secondary.

Results based on main income covariate revealed that participants with salary and or wage had a reduced risk of getting infected with HIV (OR = 0.7061, 95% CI: 0.5215–0.9560), compared to those with no source of income.

We found that individuals who were legally married had a reduced risk of getting infected with HIV (OR = 0.3708, 95% CI: 0.1497–0.9185), compared to those who were divorced. The results also revealed that there was a higher likelihood of being infected by HIV among individuals who were diagnosed with TB (OR = 1.7986, 95% CI: 1.2473–2.5935), compared to those who never suffered from TB. We also discovered that there was a higher likelihood of getting HIV infection among participants who were diagnosed with STIs (OR = 1.6938, 95% CI: 1.2448–2.3025), compared to those who were not diagnosed with STIs.

Considering number of sexual partners, there was a higher likelihood of being HIV infected among participants who had 3 or more sexual partners (OR = 1.7647, 95% CI: 1.2751–2.4449), compared to those who had one partner. Results based on alcohol consumption showed that individuals who consumed alcohol had odds of HIV prevalence that was 1.6438 (OR = 1.6438, 95% CI: 1.3100–2.1684) times higher than those who were not consuming alcohol. Lastly, we found that using condoms as a prevention method, reduced the risk of being HIV infected (OR = 0.5516, 95% CI: 0.3482–0.8737), compared to not using condoms.

The results above indicate that age group, education levels, source of income and marital status, along with behaviours like alcohol use, condom use and having multiple sexual partners, are the key predictors of HIV prevalence. Also, being diagnosed with sexually transmitted infections (STIs) and TB increases the chances of getting infected with HIV.

After fitting the model, the smoothed HIV prevalence rates were calculated and are displayed in Figure 4.

Comparing the HIV prevalence intervals in Figure 2 (unsmoothed prevalence rates) and Figure 4 (smoothed prevalence rates), we observe that the intervals differ. However, areas with high HIV prevalence in the unsmoothed data remain high-prevalence areas in the smoothed data, indicating consistency in spatial patterns.

Model performance was assessed using residuals plot and normal Q-Q plot for model adequacy and Moran’s I, Geary’s C statistic, and the variogram plot to evaluate spatial autocorrelation in residuals. Figure 5 below displays the residuals plot.

Figure 5 displays residuals that are symmetrically distributed around zero, showing no clear pattern, and having constant variance. The plot suggests that there is no systematic bias in the model's predictions implying that the model fits the data well. Figure 6 below displays the Q-Q plot of the residuals.

The Q-Q plot in Figure 6 shows an S-shaped pattern, indicating deviation from normality with heavier tails. However, in spatial modelling, residuals are not always expected to be normally distributed due to inherent spatial dependencies. This characteristic is well-documented in spatial statistics [31,37,43,47,48].

The global Moran’s I statistic for residuals was 0.0009971 (p = 0.4549), suggesting no significant spatial autocorrelation. This indicates that the structured model has adequately captured the spatial structure in the data.

Similarly, Geary’s C statistic was 1.0010397 (p = 0.5349), further confirming that residuals are not spatially autocorrelated. Since both Moran’s I and Geary’s C suggest no significant spatial dependence, the model appears to fit well.

Additionally, the variogram plot in Figure 7 provides strong evidence that residuals are spatially uncorrelated, further supporting the model’s adequacy.

The plot shows a flat semi variance around 0.20, indicating no spatial autocorrelation in the residuals. This suggests that the residuals are spatially independent. If spatial dependence were present, the semi variance would increase with distance, which was not observed.

The posterior density plots for statistically significant regression parameters in Figure 8 display smooth curves with single peaks, indicating stability and proper model convergence.

Based on the spatial autocorrelation tests applied, the results revealed that the residuals were not spatially autocorrelated, implying that the structured model was appropriate and had captured the spatial structure in our data. This is also supported by smooth and unimodal plots displayed by the posterior density plots.

4. Discussion

This study employed a Bayesian spatial logistic regression approach to identify the prevalence and risk factors associated with HIV/AIDS among female youth in KwaZulu-Natal, South Africa, using a structured additive model. The findings indicate significant spatial clustering of HIV prevalence, with sociodemographic, behavioural, and biological factors influencing the risk of infection.

The structured additive model revealed significant spatial variations in HIV prevalence, underscoring the influence of geographic location on HIV risk among female youth. These findings align with existing literature highlighting the clustering of HIV infections in areas with limited access to healthcare services, high population densities, and socio-economic disparities [49]. Addressing spatial inequalities requires targeted interventions in high-risk areas, such as rural and peri-urban settings, where female youth may face barriers to accessing sexual and reproductive health services.

The study identified age as a key determinant of HIV prevalence. Participants aged 20–24, 25–29, and 30–34 had significantly higher odds of HIV infection compared to those aged 15–19. This reflects the high burden of HIV among young adult women, often driven by power imbalances in relationships, transactional sex, and limited access to preventive services [8,50,51].

Education emerged as another critical factor. Lower educational attainment was strongly associated with higher HIV prevalence, while tertiary education was protective. These findings highlight the role of education in empowering young women with knowledge about HIV prevention and increasing their ability to make informed decisions about their sexual health [52,53,54]. Socioeconomic factors, such as income source, were also significant. Female youth earning a salary or wage were less likely to be HIV positive, emphasizing the protective role of financial independence and economic empowerment [55,56].

Risky sexual behaviours, including multiple sexual partners and inconsistent condom use, were significant predictors of HIV infection. These findings align with studies showing that such behaviours amplify the risk of HIV transmission in high-prevalence settings [57,58]. Alcohol use was also associated with higher odds of HIV infection, consistent with evidence that alcohol impairs judgment and increases engagement in risky sexual behaviours [59,60,61].

Health factors, including co-infections with TB and STIs, were strongly associated with HIV infection. These co-morbidities exacerbate vulnerability to HIV, emphasizing the need for integrated healthcare approaches addressing HIV and related infections [62,63,64].

Legal marriage and condom use were protective against HIV infection. Condom use remains one of the most effective strategies to prevent HIV transmission [56,65]. The protective effect of marriage may be due to reduced exposure to high-risk sexual networks, although this depends on the stability and fidelity of the marital relationship [66,67].

The identification of two HIV hotspots in the study area around Greater Endendale and areas around Nadi and KwaMbanjwa is crucial for understanding spatial disparities in HIV prevalence. These hotspots highlight the geographic concentration of HIV prevalence, suggesting that local factors such as socio-economic conditions, healthcare access, and behavioural risk factors may contribute to higher transmission rates in specific areas [54,58,61]. Identifying these clusters is vital for targeted public health interventions, such as focused awareness programs, HIV testing, and prevention efforts, particularly in regions with higher risk [63,64]. By understanding the spatial distribution of HIV, resources can be more effectively allocated, ultimately reducing the burden of the epidemic in high-prevalence regions [67].

5. Conclusions

This study highlights the significant spatial variation in HIV prevalence among female youth in KwaZulu-Natal, South Africa, and identifies key demographic, behavioural, and health factors associated with increased risk. The Bayesian spatial logistic regression approach allowed for the integration of spatial effects and covariates, providing a nuanced understanding of HIV risk in this population. Key findings include the high odds of HIV prevalence among young adults aged 20–34, particularly those with lower educational attainment and limited economic opportunities. Risky behaviours, such as alcohol use, multiple sexual partnerships, and inconsistent condom use, further compound the vulnerability of female youth. Additionally, co-infections with TB and STIs significantly increase the likelihood of HIV infection. The protective role of tertiary education, salaried employment, condom use, and legal marriage emphasizes the importance of multi-sectoral approaches to HIV prevention. Strategies should focus on improving education, promoting economic empowerment, increasing access to HIV prevention services, and addressing co-infections in integrated healthcare programs. Given the spatial heterogeneity observed, interventions must be geographically targeted to address localized drivers of HIV prevalence. Future studies should build on these findings by incorporating longitudinal data to assess the causal pathways between socio-demographic factors, spatial effects, and HIV prevalence. By tailoring interventions to the unique needs of female youth in specific locations, policymakers can enhance the effectiveness of HIV prevention efforts and reduce the burden of the epidemic in KwaZulu-Natal.