Analysis of HIV/AIDS Epidemic and Socioeconomic Factors in Sub-Saharan Africa

Abstract

Sub-Saharan Africa has been the epicenter of the outbreak since the spread of acquired immunodeficiency syndrome (AIDS) began to be prevalent. This article proposes several regression models to investigate the relationships between the HIV/AIDS epidemic and socioeconomic factors (the gross domestic product per capita, and population density) in ten countries of Sub-Saharan Africa, for 2011–2016. The maximum likelihood method was used to estimate the unknown parameters of these models along with the Newton–Raphson procedure and Fisher scoring algorithm. Comparing these regression models, there exist significant spatiotemporal non-stationarity and auto-correlations between the HIV/AIDS epidemic and two socioeconomic factors. Based on the empirical results, we suggest that the geographically and temporally weighted Poisson autoregressive (GTWPAR) model is more suitable than other models, and has the better fitting results.

1. Introduction

Acquired immunodeficiency syndrome (AIDS) is a malignant infectious disease with a high fatality rate caused by human immunodeficiency virus (HIV). The HIV/AIDS epidemic has been one of the greatest global public health and social development problems since 1981, particularly in Sub-Saharan Africa. As of 31 December 2016, over 30 million people had died from the disease [1]. More than 70% of the 35 million people are infected with the HIV/AIDS disease in Sub-Saharan Africa. Thus, the HIV/AIDS epidemic of Sub-Saharan Africa has attracted extensive attention from researchers around the world [2,3,4].

In earlier studies, Janet et al. [5] and Hallman et al. [6] demonstrated the relationship between the disease and socioeconomic status. Chris et al. [7] indicated socioeconomic factors to explain this disease outperformed cultural ones in South Africa. Mathematical models always play an important role in evaluating the trends of the HIV/AIDS epidemic [8]. For example, regression models have been widely used in the study of the relationship between this disease and influencing factors. Shiboski et al. [9] considered a generalized linear model to obtain the statistical analysis of the HIV/AIDS disease. A mixed-effects linear regression model was used to analyze the correlation between national population and antenatal care [10]. Laurence et al. [11] applied a spatial regression model to show that the epidemic had substantial geographic variance across Sub-Saharan Africa.

This paper proposes several regressive models to investigate the relationships between the HIV/AIDS epidemic, the gross domestic product (GDP) per capita and the population density in ten countries of Sub-Saharan Africa. The Poisson regression model is introduced in Section 2.1. Section 2.2 and Section 2.3 describe two spatial models, respectively. A spatiotemporal autoregressive model is proposed in Section 2.4. The maximum likelihood method is used to obtain the iterative formulas of coefficient estimations in Section 3. The main results are shown in Section 4, followed by discussion in Section 5.

2. Methodologies

2.1. Poisson Regression Model

Regression models are a set of statistical processes for estimating the relationships between response and explanatory variables. The classical model is a linear regression. Nelder and Wedderburn [12] extended the linear model to a generalized linear regression for solving the discrete data problem. This kind of models are very important in ecology, medicine and economics [13,14,15]. Suppose that $Y=(Y_1,Y_2,\dots ,Y_n)$ is the response variable, where $Y_i(i=1,\dots ,n)$ are independent. The density function is

$$\begin{array}{c}\hfill f(y_i;{\theta }_i,{\varphi }_i)=exp\left(\frac{y_i{\theta }_i-b({\theta }_i)}{a({\varphi }_i)}+c(y_i,{\varphi }_i)\right),\end{array}$$

where $a(·),b(·),c(·,·)$ are known functions, and ${\theta }_i,{\varphi }_i$ are unknown parameters for $i=1,2,\dots ,n$. Denote ${\mu }_i=E(Y_i)$, and $g({\mu }_i)=ln({\mu }_i)$ is a link function. Let $X_{ij}$ be explanatory variables for the ith observation in the jth variable. Then, the Poisson regression (PR) model is given by

$$\begin{array}{c}\hfill g({\mu }_i)\triangleq {\eta }_i=\sum _{j=1}^p{\beta }_j{X}_{ij},\end{array}$$

(1)

where $i=1,2,\dots ,n$, and ${\beta }_j(j=1,2,\dots ,p)$ are unknown parameters.

2.2. Geographically Weighted Poisson Regression Model

With in-depth study, regression models have been frequently applied in epidemiology and health geography for trying to investigate the persistent geographical variations in disease [16]. Based on the generalized linear regression, Brunsdon et al. [17] proposed the geographically weighted regression model to analyze the spatial non-stationary processes of discrete data. The disease maps arising from this process are considered through the establishment of the geographically weighted Poisson regression (GWPR) model [18,19,20] below

$$\begin{array}{c}\hfill g({\mu }_i)\triangleq {\eta }_i=\sum _{j=1}^p{\beta }_j(u_i,v_i)X_{ij},\end{array}$$

(2)

where $(u_i,v_i)(i=1,2,\dots ,n)$ are the geographical locations, and ${\beta }_j(u_i,v_i)(j=1,2,\dots ,p)$ are unknown parameters at the position $(u_i,v_i)$.

2.3. Geographically Weighted Poisson Autoregressive Model

Another issue deserving of special attention is whether there exists an interaction between different regions in terms of spatial data. Previous studies [21,22,23,24] showed that spatial data has not only spatial non-stationarity but also correlation. Zhang [25] proposed the geographically weighted Poisson autoregressive (GWPAR) model as follows:

$$\begin{array}{c}\hfill g({\mu }_i)\triangleq {\eta }_i=\rho \sum _{k=1}^n{c}_{ik}{\eta }_k+\sum _{j=1}^p{\beta }_j(u_i,v_i)X_{ij},\end{array}$$

(3)

where $\rho$ is a scalar autoregressive parameter, and $c_{ik}(i,k=1,2,\dots ,n)$ is the adjacency relation between the ith and kth locations. Let $c^i$ be the number of regions adjacent to the ith position. If the kth position is next to the ith’s, then $c_{ik}=1/c^i$. Otherwise, $c_{ik}=0$.

2.4. Geographically and Temporally Weighted Poisson Autoregressive Model

Recently, many spatiotemporal models have been proposed to describe the spatiotemporal variations in the relationships of response and explanatory variables [26,27]. Concerning the modeling of spatiotemporal data, there are two important properties: non-stationarity and auto-correlation. The non-stationarity indicates that there exists more than one linear relation between response and explanatory variables. It can be used to identify where interesting relationships are likely to occur or where detailed investigation is necessary in the study areas [28]. Spatiotemporal auto-correlation is an important factor to determine the temporal correlations of observations [29]. These two problems always appeared together [30]. A geographically and temporally weighted autoregressive (GTWPAR) model can be applied to account for non-stationary and auto-correlated effects simultaneously.

Let Y be the response variable, and $Y_{ik}(i=1,2,\dots ,n_k,k=1,2,\dots ,T)$ be the independent variables of Y in the ith position and the kth time. The density function can be defined as follows:

$$\begin{array}{c}\hfill f(y_{ik};{\theta }_{ik},{\varphi }_{ik})=exp\left(\frac{y_{ik}{\theta }_{ik}-b({\theta }_{ik})}{a({\varphi }_{ik})}+c(y_{ik},{\varphi }_{ik})\right),\end{array}$$

where the parameters are similar to Section 2.1. Denote ${\mu }_{ik}=E(Y_{ik})$, and $g({\mu }_{ik})=ln({\mu }_{ik})$. Let $X_{ijk}(j=1,2,\dots ,p)$ be the jth explanatory variable. The GTWPAR model is expressed by

$$\begin{array}{c}\hfill g({\mu }_{ik})\triangleq {\eta }_{ik}=\rho \sum _{m=1}^T\sum _{l=1}^{n_k}{c}_{lm}^{(ik)}{\eta }_{lm}+\sum _{j=1}^p\sum _{k=1}^T{\beta }_{jk}(u_{ik},v_{ik},t_k)X_{ijk},\end{array}$$

(4)

where $\{{\beta }_{jk}(u_{ik},v_{ik},t_k)$} is a set of unknown parameters at the ith position in the kth time, and $c_{lm}^{(ik)}$ is the adjacent relation between the location $(u_{ik},v_{ik},t_k)$ and $(u_{lm},v_{lm},t_m)$. Following the work of [31], the spatiotemporal distance between the locations $(u_{ik},v_{ik},t_k)$ and $(u_{lm},v_{lm},t_m)$ can be defined as

$$d_{lm}^{(ik)}=\sqrt{\lambda [{(u_{ik}-u_{lm})}^2+{(v_{ik}-v_{lm})}^2]+\mu {(t_k-t_m)}^2,}$$

where $\mu$ and $\lambda$ are used to balance spatiotemporal distances. Suppose that

$$c_{lm}^{(ik)}=\left\{\begin{array}{cc}1/c^{ik},\hfill & 0<d_{lm}^{(ik)}<d,\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

where d is a constant and satisfies $min\left\{d_{lm}^{(ik)}\right\}<d<max\left\{d_{lm}^{(ik)}\right\}$.

Next, we rewrite the model (4) in a matrix form

$$\begin{array}{c}\hfill \eta =\rho C\eta +{\mathbf{B}}^{\prime }{X}^{\prime },\end{array}$$

where $\eta ={({\eta }_{11},\cdots ,{\eta }_{n_{1}1},{\eta }_{12},\cdots ,{\eta }_{n_{2}2},\cdots ,{\eta }_{1T},\cdots ,{\eta }_{n_{T}T})}^{{}^{\prime }}$, $C=(c_{lm}^{(ik)})$, $X=(X_{ijk})$ and $\mathbf{B}=({\beta }_{jk}(u_{ik},v_{ik},t_k))$. For convenience, define ${\eta }_K$ as the Kth element of $\eta$; $C_{IK}$ and $X_{IK}$ are the Ith row and the Kth column of the matrices C and X, respectively. The detailed expressions of C, X and $\mathbf{B}$ are given in Appendix A.1.

Remark 1.

For the GTWPAR model (4), if $\rho =0$ and ${\beta }_{jk}(u_{ik},v_{ik},t_k)$ is independent of the spatiotemporal effect, the model is a PR model. If $\rho =0$ and ${\beta }_{jk}(u_{ik},v_{ik},t_k)$ is dependent on spatial effect but independent of temporal effect, the model becomes GWPR model. If $\rho \ne 0$ and ${\beta }_{jk}(u_{ik},v_{ik},t_k)$ is independent of temporal effect, it is the GWPAR model. Thus, PR, GWPR and GWPAR models are the special cases of the GTWPAR model.

3. Coefficient Estimation

In this section, we only provide the estimation method of the GTWPAR model since the PR, GWPR and GWPAR models are its special cases (Remark 1). Let $(u_{ik},v_{ik},t_k)(i=1,2,\dots ,n_k,k=1,2,\dots ,T)$ be any point in the studied spatiotemporal region. We fix a point $(u_{00},v_{00},t_0)$ and assume that ${\beta }_{jk}(u_{ik},v_{ik},t_k)\approx {\beta }_{j0}(u_{00},v_{00},t_0)(j=1,2,\dots ,p)$. Then, the model (4) can be rewritten by

$$\begin{array}{c}\hfill {\eta }_{ik}=g({\mu }_{ik})=\rho \sum _{m=1}^T\sum _{l=1}^{n_k}{c}_{lm}^{(ik)}{\eta }_{lm}+\sum _{j=1}^p\sum _{k=1}^T{\beta }_{j0}(u_{00},v_{00},t_0)X_{ijk}.\end{array}$$

(5)

Denote $\beta (u_{00},v_{00},t_0)={({\beta }_{10},\dots ,{\beta }_{p0})}^{\prime }$, $\mathbf{X}=diag(X_{i.})$ and $X_{i·}=(X_{i1},\dots ,X_{ip})$. The corresponding matrix form can be represented as $\eta =\rho C\eta +{\beta }^{\prime }(u_{00},v_{00},t_0){\mathbf{X}}^{\prime }$.

3.1. Estimation of Parameter Vector $\beta$

For the fixed point $(u_{00},v_{00},t_0)$, we define a spatiotemporal distance $d_{ik}^{(0)}$ from this point to $(u_{ik},v_{ik},t_k)$ as $d_{ik}^{(0)}=\sqrt{\lambda [{(u_{00}-u_{ik})}^2+{(v_{00}-v_{ik})}^2]+\mu {(t_0-t_k)}^2}.$ The Gauss kernel function of these two points can be written by

$$\begin{array}{ccc}\hfill w_{ik}(u_{00},v_{00},t_0)& =& \frac{1}{\sqrt{2\pi }}exp\left\{-\frac{1}{2}{\left(\frac{d_{ik}^{(0)}}{h_{ST}}\right)}^2\right\}\hfill \\ & =& \frac{1}{\sqrt{2\pi }}exp\left\{-\frac{1}{2}\frac{\lambda [{(u_{00}-u_{ik})}^2+{(v_{00}-v_{ik})}^2]+\mu {(t_0-t_k)}^2}{h_{ST}^2}\right\}\hfill \\ & =& \frac{1}{\sqrt{2\pi }}exp\left\{-\frac{1}{2}\left(\frac{{(u_{00}-u_{ik})}^2+{(v_{00}-v_{ik})}^2}{h_S^2}+\frac{{(t_0-t_k)}^2}{\tau h_S^2}\right)\right\},\hfill \end{array}$$

where $h_{ST}$ and $h_S$ are the space-time bandwidth and space bandwidth, respectively. Meanwhile, we have $h_{ST}^2=\lambda h_S^2$, and $\tau =\lambda /\mu$ is a spatiotemporal factor. Without loss of generality, let $\lambda =1$. Then, the weighted maximum likelihood of $Y_{ik}(i=1,2,\dots ,n_k,k=1,2,\dots ,T)$ at the point $(u_{00},v_{00},t_0)$ is

$$\begin{array}{ccc}& & L({\beta }_{10},{\beta }_{20},\cdots ,{\beta }_{p0})=\prod _{k=1}^T\prod _{i=1}^{n_k}f(y_{ik};{\theta }_{ik},{\varphi }_{ik})w_{ik}(u_{00},v_{00},t_0),\hfill \end{array}$$

where $f(y_{ik};{\theta }_{ik},{\varphi }_{ik})$ is the density function. The log-likelihood can be obtained as follows:

$$\begin{array}{c}\hfill {\mathbf{L}}_1(\beta (u_{00},v_{00},t_0))=\sum _{k=1}^T\sum _{i=1}^{n_k}\left(\frac{y_{ik}{\theta }_{ik}-b({\theta }_{ik})}{a({\varphi }_{ik})}+c(y_{ik},{\varphi }_{ik})\right)w_{ik}(u_{00},v_{00},t_0).\end{array}$$

Note that $c(y_{ik},{\varphi }_{ik})=-ln(y_{ik}!)$, $b({\theta }_{ik})={\mu }_{ik}=exp({\theta }_{ik})$, and $a({\varphi }_{ik})={\varphi }_{ik}=1$. Thus, $E(Y_{ik})=b^{\prime }({\theta }_{ik})=exp({\theta }_{ik})={\mu }_{ik},Var(Y_{ik})=b^{″}({\theta }_{ik})a({\varphi }_{ik})=exp({\theta }_{ik})={\mu }_{ik}.$ Differentiating ${\mathbf{L}}_1$ with respect to $\beta (u_{00},v_{00},t_0)$ yields

$$\begin{array}{c}\hfill \frac{\partial {\mathbf{L}}_1}{\partial {\beta }_{r0}}=\sum _{k=1}^T\sum _{i=1}^{n_k}\left(\frac{y_{ik}-{\mu }_{ik}}{a_{ik}\varphi }\frac{\partial {\theta }_{ik}}{\partial {\beta }_{r0}}\right)w_{ik}(u_{00},v_{00},t_0)=0,\end{array}$$

(6)

where ${\beta }_{r0}={\beta }_r(u_{00},v_{00},t_0)(r=1,2,\cdots ,p)$, and

$$\begin{array}{c}\hfill \frac{\partial {\theta }_{ik}}{\partial {\beta }_{r0}}={\left(\frac{\partial {\mu }_{ik}}{\partial {\theta }_{ik}}\right)}^{-1}\frac{\partial {\mu }_{ik}}{\partial g({\mu }_{ik})}\frac{\partial g({\mu }_{ik})}{\partial {\beta }_{r0}}=\frac{1}{b^{″}({\theta }_{ik})}\frac{1}{g^{\prime }({\mu }_{ik})}\frac{\partial {\eta }_{ik}}{\partial {\beta }_{r0}}.\end{array}$$

For convenience, let $N={\sum }_{k=1}^N{n}_k$ and $W={(w_{ik}(u_{00},v_{00},t_0))}_{N\times N}$. Denote $A={(I_N-\rho C)}^{-1}$, $Y={(Y_{11},\cdots ,Y_{n_{1}1},\cdots ,Y_{1T},\cdots ,Y_{n_{T}T})}^{{}^{\prime }}$, $\mu ={({\mu }_{11},\cdots ,{\mu }_{n_{1}1},\cdots ,{\mu }_{1T},\cdots ,{\mu }_{n_{T}T})}^{{}^{\prime }}$, $\theta ={({\theta }_{11},\cdots ,{\theta }_{n_{1}1},\cdots ,{\theta }_{1T},\cdots ,{\theta }_{n_{T}T})}^{{}^{\prime }}$, $\varphi ={({\varphi }_{11},\cdots ,{\varphi }_{n_{1}1},\cdots ,{\varphi }_{1T},\cdots ,{\varphi }_{n_{T}T})}^{{}^{\prime }}$. Suppose that $Y_K$, ${\mu }_K$, ${\theta }_K$ and ${\varphi }_K$ are the Kth elements of Y, $\mu$, $\theta$ and $\varphi$, respectively. Then, we take the derivative of the model (5) with respect to ${\beta }_{r0}$, and obtain

$$\begin{array}{c}\hfill \frac{\partial {\eta }_l}{\partial {\beta }_{r0}}=\sum _{h=1}^N{A}_{lh}{X}_{hr}=A_{l·}{X}_{·r},l=1,2,\dots ,N.\end{array}$$

The calculation process is given in Appendix A.2. Thus, the Equation (6) can be rewritten as

$$\begin{array}{c}\hfill \frac{\partial {\mathbf{L}}_1}{\partial {\beta }_{r0}}=\frac{1}{\varphi }\sum _{l=1}^N{T}_l{A}_{l·}{X}_{·r}(Y_l-{\mu }_l)g^{\prime }({\mu }_l)W_l(u_{00},v_{00},t_0)=0.\end{array}$$

However, there is not a close-form solution for $\beta (u_{00},v_{00},t_0)$. The Newton–Raphson procedure and Fisher scoring algorithm are used to get the estimation of $\beta$. The iterative formula is expressed as

$$\begin{array}{ccc}\hfill {\widehat{\beta }}^{(m+1)}(u_{00},v_{00},t_0)& =& {\widehat{\beta }}^{(m)}(u_{00},v_{00},t_0)+{\mathbf{I}}^{-1}({\widehat{\beta }}^{(m)}(u_{00},v_{00},t_0))S({\widehat{\beta }}^{(m)}(u_{00},v_{00},t_0))\hfill \\ & =& {({(A^{(m)}X)}^{\prime }{T}^{(m)}W(u_{00},v_{00},t_0)(A^{(m)}X))}^{-1}\hfill \\ & & \times {(A^{(m)}X)}^{\prime }{T}^{(m)}W(u_{00},v_{00},t_0)Z^{(m)},\hfill \end{array}$$

(7)

where the Fisher information matrix $\mathbf{I}(\beta )=E(I(\beta ))$, and

$$S({\widehat{\beta }}^{(m)}(u_{00},v_{00},t_0))={\left(\frac{\partial {\mathbf{L}}_1}{\partial {\beta }_{10}},\frac{\partial {\mathbf{L}}_1}{\partial {\beta }_{20}},\cdots ,\frac{\partial {\mathbf{L}}_1}{\partial {\beta }_{p0}}\right)}^{\prime }$$

is the scalar vector. The detail process is provided in Appendix A.2. For the fixed point $(u_{ik},v_{ik},t_k)(i=1,2,\dots ,n_k;k=1,2,\dots ,T)$, ${\widehat{\beta }}_{jk}(u_{ik},v_{ik},t_k)$ can be obtained by (7).

Remark 2.

The estimations $\widehat{\beta }(u_{ik},v_{ik},t_k)(i,l=1,2,\dots ,n_k,k,m=1,2,\dots ,T)$ are related to the temporal and spatial effects in the GTWPAR model. If $m\ne k$, $w_{ik}(u_{lm},v_{lm},t_m)=0$ and $c_{lm}^{(ik)}=0$, then $\widehat{\beta }(u_{ik},v_{ik},t_k)=\widehat{\beta }(u_i,v_i)$ correspond to the parameter estimations of the GWPAR model. If $w_{ik}(u_{lm},v_{lm},t_m)=0(m\ne k)$ and $C=\mathbf{0}$, they are the estimations of the GWPR model. If $W=\mathbf{0}$ and $C=\mathbf{0}$, then $\widehat{\beta }(u_{ik},v_{ik},t_k)=\widehat{\beta }$ are the global estimation values of the PR model.

3.2. Estimation of Parameter $\rho$

Based on the density function, the log-likelihood function of $\rho$ is

$$\begin{array}{c}\hfill {\mathbf{L}}_2(\rho )=\sum _{k=1}^T\sum _{i=1}^{n_k}\left(\frac{y_{ik}{\theta }_{ik}-b({\theta }_{ik})}{a({\varphi }_{ik})}+c(y_{ik},{\varphi }_{ik})\right).\end{array}$$

Differentiating ${\mathbf{L}}_2(\rho )$ with respect to $\rho$, we have

$$\begin{array}{c}\hfill \frac{\partial {\mathbf{L}}_2}{\partial \rho }=\sum _{k=1}^T\sum _{i=1}^{n_k}\left(\frac{y_{ik}-{\mu }_{ik}}{a_{ik}\varphi }\frac{\partial {\theta }_{ik}}{\partial \rho }\right)=0,\end{array}$$

(8)

where $\frac{d{\theta }_{ik}}{d\rho }=\frac{1}{b^{″}({\theta }_{ik})g^{\prime }({\mu }_{ik})}\frac{d{\eta }_{ik}}{d\rho }.$ Then, we take the derivative of the model (5) with respect to $\rho$ as follows:

$$\begin{array}{c}\hfill \frac{d{\eta }_l}{d\rho }=\frac{dg({\mu }_l)}{d\rho }=\sum _{h=1}^N{A}_{l·}{C}_{·h}{\eta }_h.\end{array}$$

The detail calculation is given in Appendix A.3. Then, the Equation (8) can be rewritten in the following nonlinear form

$$\frac{d{\mathbf{L}}_2}{d\rho }=\sum _{l=1}^N\frac{(Y_l-{\mu }_l){\sum }_{h=1}^N{A}_{l·}{C}_{·h}{\eta }_h}{a_l\varphi V({\mu }_l)g^{\prime }({\mu }_l)}=0.$$

According to the Newton–Raphson procedure and Fisher scoring algorithm, the iterative formula of ${\widehat{\rho }}^{(m+1)}$ is

$$\begin{array}{ccc}\hfill {\widehat{\rho }}^{(m+1)}& =& {\widehat{\rho }}^{(m)}+{\mathcal{I}}^{-1}({\widehat{\rho }}^{(m)})S({\widehat{\rho }}^{(m)})\hfill \\ & =& {\displaystyle {\widehat{\rho }}^{(m)}+{({(A^{(m)}C{\eta }^{(m)})}^{\prime }{T}^{(m)}(A^{(m)}C{\eta }^{(m)}))}^{-1}}\hfill \\ & & \times {(A^{(m)}C{\eta }^{(m)})}^{\prime }{T}^{(m)}(Z^{(m)}-{\eta }^{(m)}),\hfill \end{array}$$

(9)

where the scalar vector $S({\widehat{\rho }}^{(m)})=\frac{1}{\varphi }{(AC\eta )}^{\prime }T(Z-\eta )$ and the Fisher information matrix $\mathcal{I}(\rho )=\frac{1}{\varphi }{(AC\eta )}^{\prime }T(AC\eta )$. The calculation process of the scalar vector $S({\widehat{\rho }}^{(m)})$ and the information matrix $\mathcal{I}$ is given in Appendix A.3.

4. Main Results

In this section, we apply the PR, GWPR, GWPAR and GTWPAR models to analyze the relationships between the HIV/AIDS epidemic, the GDP per capita and population density in ten countries of Sub-Saharan Africa from 2011 to 2016. The ten countries are Angola, Botswana, Lesotho, Malawi, Mozambique, Namibia, South Africa, Swaziland, Zimbabwe and Zambia. The parameters of these four models are estimated by the Newton–Raphson procedure and Fisher scoring algorithm. The coefficient of determination R${}^2$, the corrected Akaike information criterion (AICc), the deviation (D) and mean-square error (MSE) are used to compare the performances of the four models [18].

4.1. The HIV/AIDS Epidemic Models

The data of HIV/AIDS incidence, GDP per capita and population density were derived from http://data.cnki.net/InternationalData/Report. Readers should note that authorization is required to access the database on this website. Figure 1 describes the HIV/AIDS incidence in ten countries from 2011 to 2016. It shows that the incidence varies significantly in different regions. Angola has a minimum incidence of less than 5%, while Botswana and Swaziland have higher incidences of more than 20% every years. Therefore, it may be necessary to consider the temporal and spatial factors in analyzing the HIV/AIDS epidemic.

The distributions of HIV/AIDS cases, GDP per capita and population density are displayed in Figure 2. The Pearson correlation coefficients between these cases and GDP per capita and population density are 0.2739 and $-0.1179$, respectively. Meanwhile, the two socioeconomic factors have different effects on the HIV/AIDS cases at the spatiotemporal locations. These reflect a spatiotemporal non-stationarity between the cases and two factors in ten countries from 2011 to 2016.

Table 1. p-values of the spatial and temporal autocorrelation analysis.

Time	2011	2012	2013	2014	2015	2016
p-value	0.0024	0.0028	0.0021	0.0019	0.0017	0.0022
Regions	Angola	Botswana	Lesotho	Malawi	Mozambique	Namibia
p-value	0.9094	0.8807	0.0092	0.5300	0.1289	0.8267
	South Africa	Swaziland	Zimbabwe	Zambia
	0.0045	0.1491	0.5688	0.2231

lists the p-values of the first-order autocorrelation of HIV/AIDS cases in the different years of the same region or the different regions of the same year. Each region has a significant spatial autocorrelation (p-value < 0.01) each year. Lesotho and South Africa had temporal autocorrelation during 2011 to 2016. Thus, the spatial and temporal autocorrelation should not be ignored.

Next, we standardized the two socioeconomic factors. The multiplex collinear test [32] was performed by the condition number $k=\sqrt{{\lambda }_{max}/{\lambda }_{min}}=1.804(\le 15)$ ($\lambda$ is the eigenvalue of explanatory variable matrix). If $k>15$, then the data have collinearity. Otherwise, there is no collinearity. Thus, there is no collinearity between the two factors. Let ${\mu }_{ik}$, $r_{ik}$ and $P_{ik}$ be the annual HIV/AIDS cases (Unit: 1/1000 people), incidence (Unit: 1/100) and total population (unit: 100,000 people) in the kth year of the ith region, respectively. Denote $g({\mu }_{ik})={\eta }_{ik}=ln{\mu }_{ik}=ln{r}_{ik}+ln{P}_{ik}({\mu }_{ik}=r_{ik}{P}_{ik},i=1,2,\dots ,10$ and $k=1,2,\dots ,6$). Let $X_{i1k}$ and $X_{i2k}$ be the GDP per capita and population density in the ith region at the kth year, respectively. The PR model is written by

$$\begin{array}{c}\hfill g({\mu }_{ik})={\beta }_0+{\beta }_1{X}_{i1k}+{\beta }_2{X}_{i2k},i=1,2,\dots ,10,k=1,2,\dots ,6,\end{array}$$

(10)

where ${\beta }_j(j=0,1,2)$ are unknown constants. The GWPR model is introduced as

$$\begin{array}{c}\hfill g({\mu }_{ik})={\beta }_0(u_{ik},v_{ik})+{\beta }_1(u_{ik},v_{ik})X_{i1k}+{\beta }_2(u_{ik},v_{ik})X_{i2k},i=1,2,\dots ,10,\end{array}$$

(11)

where k is a fixed constant taken from $\{1,2,\dots ,6\}$, and ${\beta }_j(u_{ik},v_{ik})$ are unknown spatial parameters for the ith country $(u_{ik},v_{ik})$ in the kth year. Let $\rho$ be a scalar autoregressive parameter, and $c_{il}$ be a constant that represents an adjacency relation. The GWPAR model is

$$\begin{array}{c}\hfill g({\mu }_{ik})=\rho \sum _{l=1}^n{c}_{il}{\eta }_{ik}+{\beta }_0(u_{ik},v_{ik})+{\beta }_1(u_{ik},v_{ik})X_{i1k}+{\beta }_2(u_{ik},v_{ik})X_{i2k},\end{array}$$

(12)

where $n=10$, k is a fixed constant, and ${\beta }_j(u_{ik},v_{ik})$ are defined as above. Let $c_{lm}^{(ik)}$ be a spatiotemporal adjacency relation, and ${\beta }_{jk}(u_{ik},v_{ik},t_k)(k=1,2,\dots ,6)$ be unknown spatiotemporal parameters in the ith country $(u_{ik},v_{ik})$ in the kth year. The GTWPAR model is established as follows:

$$\begin{array}{ccc}\hfill g({\mu }_{ik})& =& \rho \sum _{m=1}^T\sum _{l=1}^{n_k}{c}_{lm}^{(ik)}{\eta }_{lm}+{\beta }_{0k}(u_{ik},v_{ik},t_k)\hfill \\ & & +{\beta }_{1k}(u_{ik},v_{ik},t_k)X_{i1k}+{\beta }_{2k}(u_{ik},v_{ik},t_k)X_{i2k},\hfill \end{array}$$

(13)

where $T=6$; $n_k=10$ for every k years; and $\rho$ is defined as above.

Algorithms I, II, III and IV of PR, GWPR, GWPAR and GTWPAR models are provided in Appendix A.4, respectively.

4.2. Statistical Analysis

For the PR model, we get the estimated values of unknown parameters by Algorithm I. Then, the best space bandwidth is chosen by the cross-validation method. Following Huang et al. [28], the range [0.09, 2.49] of the space bandwidth is selected according to the minimum and maximum distance of the geographical positions. In the GWPR model, the best space bandwidth is $h=0.62,0.59,0.62,0.61,0.60,0.60$, and the estimations of coefficient functions are given by Algorithm II. The optimal space bandwidth of the GWPAR model is selected as $h=1.2895,1.1316,1.1316,1.0526,1.0526,1.0526.$ Based on Algorithm III, we can get the estimations of coefficient functions and the scalar autoregressive parameter $\widehat{\rho }=0.267,0.269,0.263,0.264,0.264,0.264.$ For the GTWPAR model, we chose $h_s=1.1316,0.9737,0.8947,0.9211,0.7842,0.8789$ and $\tau =0.1$, where $\tau (>0)$ is a balanced parameter. The coefficient estimations and scalar autoregressive parameter $\widehat{\rho }=0.126$ can be obtained by Algorithm IV. The quantile and mean values of coefficient estimations and response variables are shown in

Table 2. The quantile and mean values of coefficient estimations and response variables.

Model	Coefficient	Min	1st Qu	Median	3rd Qu	Max	Mean
True	$\eta$	5.293	5.644	6.455	7.288	8.871	6.559
PR	${\widehat{\beta }}_1$	0.581	0.581	0.581	0.581	0.581	0.581
	${\widehat{\beta }}_2$	0.385	0.385	0.385	0.385	0.385	0.385
	$\widehat{\eta }$	6.344	6.527	7.185	7.447	8.162	7.078
GWPR	${\widehat{\beta }}_1$	−0.519	−0.249	0.009	0.641	1.774	0.288
	${\widehat{\beta }}_2$	−0.672	−0.139	0.068	2.044	3.803	0.747
	$\widehat{\eta }$	5.618	6.422	7.022	7.202	8.813	6.932
GWPAR	${\widehat{\beta }}_1$	−6.286	−1.258	−0.910	0.240	3.036	−0.556
	${\widehat{\beta }}_2$	−2.089	−0.669	0.019	2.281	9.152	0.912
	$\widehat{\eta }$	5.176	5.625	6.435	7.219	8.940	6.539
GTWPAR	${\widehat{\beta }}_1$	−5.642	−0.956	−0.677	0.186	2.992	−0.443
	${\widehat{\beta }}_2$	−1.865	−0.464	0.215	2.134	9.104	0.988
	$\widehat{\eta }$	5.200	5.589	6.387	7.215	8.803	6.489

. We note that the GWPR, GWPAR and GTWPAR models can reflect the non-stationarity property of the influencing factors; the PR model cannot. Moreover, the GTWPAR model has a better performance than other models by comparing the true and fitted values.

The average estimated coefficients are visualized in Figure 3. For the PR model, the GDP per capita and population density had the same effect on the HIV/AIDS epidemic for ten countries in six years. However, there exist significant spatial non-stationarity and auto-correlation for different countries under the GWPR, GWPAR and GTWPAR models. Figure 4 shows the spatial distribution of the average MSE of their response variables. The lighter the color, the smaller the average error is. Thus, the GWPAR and GTWPAR models have the better fitting results.

These four indicators can effectively compare the performances of the proposed models (

Table 3. The comparison of the four models.

Model	R${}^2$	AICc	D	MSE
PR	0.1291	42,504.40	42,624.07	1.6488
GWPR	0.6139	6495.12	6613.02	0.4326
GWPAR	0.9940	155.46	236.02	0.0067
GTWPAR	0.9957	115.25	190.05	0.0048

). The calculation formulas of R${}^2$, AICc, D and MSE are given in Appendix A.5. The coefficient of determination R${}^2$ gradually increases from 12.91% of the PR model to 99.57% of the GTWPAR model. The MSE, AICc and D values of the GTWPAR model are smaller than those of other models. Therefore, the GTWPAR model is more suitable to investigate the spatiotemporal HIV/AIDS epidemic.

Based on the GTWPAR model, the mean values and 95% confidence intervals of the coefficient estimations are shown in Figure 5. The mean estimations are represented by the dot, and the 95% confidence intervals are given by the upper and lower lines. Note that the GDP per capita in Botswana, Namibia and South Africa has a positive effect on the HIV/AIDS cases. Six other countries (except Lesotho) had the opposite results. The population density for five countries had a positive effect on the HIV/AIDS cases—Angola, Botswana, Namibia, South Africa and Zambia. The population density of other five countries had the negative effect. Moreover, the impact of the GDP per capita on HIV/AIDS epidemic had a strong spatiotemporal non-stationarity in Lesotho, Malawi and Zimbabwe, while the population density had a strong spatiotemporal non-stationarity in Angola.

5. Conclusions

In this paper, we propose four regression models, including the PR, GWPR, GWPAR and GTWPAR, to investigate the non-stationary and auto-correlation properties. The relationships between the HIV/AIDS epidemic, GDP per capita and population density were analyzed in ten countries of Sub-Saharan Africa from 2011 to 2016. The unknown parameters of these models can be estimated by the Newton–Raphson procedure and Fisher scoring algorithm.

The PR model is a classical generalized model, which considers the global relationships between the response and explanatory variables. The GWPR and GWPAR models have been introduced to determine the spatial non-stationarity or auto-correlation. The GTWPAR model proposed by this article can be used to investigate not only spatiotemporal non-stationary but also auto-correlation. Thus, the PR, GWPR and GWPAR models are several special cases of the GTWPAR model (see Remark 1 and Remark 2). The performances of these models were evaluated by analyzing the correlations between the HIV/AIDS epidemic and two socioeconomic factors. The parameter estimations of the models can be obtained by Algorithms I, II, III and IV in Appendix A.4.

The results show that the impacts of GDP per capita and population density on HIV /AIDS cases had significant spatiotemporal non-stationarity and auto-correlation. The GWPR, GWPAR and GTWPAR models can reflect the strong spatial or spatiotemporal non-stationarity. The auto-correlation can be reflected in the GWPAR and GTWPAR models. Compared with other models, the GTWPAR model is more effective in terms of four comparison indicators. Thus, we suggest that the GTWPAR model can be used to analyze the spatiotemporal characteristics of the HIV/AIDS epidemic and the influences of the GDP per capita and population density.

Further work also exists in our study. For example, we observed that the effects of the GDP per capita for Lesotho, Malawi and Zimbabwe and the population density for Angola on HIV/AIDS had strong spatiotemporal non-stationarity. These may be the result of local environmental or political factors. Whether the fitting results of these regions will perform better if explanatory variables such as local unique environmental or political factors are added needs to be further investigated.