**3.1 Trials**

The efficacy and safety of TAS-102 were a subject of a double-blind randomized phase II trial including Japanese patients with metastatic colorectal cancer. The patients were randomized in a 2:1 ratio to either TAS-102 or placebo arm. Median OS in the TAS-102 arm turned out to be longer than in the placebo arm (9.0 vs. 6.6 months). The most common AEs were hematologic toxicities, with grade 3 or 4 of neutropenia found in 50% of TAS-102-treated patients [28].

In a randomized double-blind phase III registration trial, TERRA, Asian patients with metastatic colorectal cancer who previously received at least two lines of systemic treatment were randomized in a 2:1 ratio to either TAS-102 or placebo group. The study confirmed the efficacy of TAS-102 in Asian population. Median OS in the TAS-102 arm was significantly longer than in the placebo group (7.8 vs. 7.1 months). Moreover, the TAS-102-treated patients had significantly lower mortality risk and significantly longer median survival follow-up time [29].

In a double-blind phase III registration trial, RECOURSE, 800 patients from Europe, North America, and Asia were randomly assigned to receive TAS-102 or placebo. PFS for the TAS-102 and placebo arm was 2.0 and 1.7 months, respectively. Patients from the TAS-102 arm had significantly longer median OS and significantly longer median time to worsening performance status (7.1 and 5.7 months, respectively) than individuals from the placebo arm (5.3 and 4.0 months, respectively). One-year OS rates for the TAS-102 and placebo arm were 27 and 18%, respectively [30]. Those promising results were further confirmed on a subgroup analysis; median OS for TAS-102-treated patients from the USA, Japan, EU, and Spain ranged between 6.5 and 7.8 months as compared with 4.3–6.7 months for the respective placebo arms, whereas median PFS in the TAS-102 and placebo groups amounted to 2.0–2.8 and 1.7–1.8 months, respectively [31, 32].

## **3.2 Predictors**

Genetic polymorphisms in homologous recombination pathway seem to be a predictor of therapeutic response in metastatic colorectal cancer patients treated with TAS-102. According to Suenaga et al., TAS-102-treated patients with a combination of ENT1 rs760370, MATE1 rs2289669, and OCT2 rs316019 single-nucleotide

**83**

**3.3 Toxicity**

*Further Therapeutic Options in Heavily Pretreated Colorectal Cancer Patients*

significantly longer OS and PFS in TAS-102-treated patients [34].

G3 or G4 neutropenia after introduction of this agent [37].

confirmed in a meta-analysis conducted by Chen et al. [39].

the KRAS status, as well as with patients' age and ethnicity [31].

factors turned out to be a predictor of therapeutic response [30].

polymorphisms of genes involved in trifluridine metabolism and thymidine phosphorylase inhibitor excretion had significantly longer PFS and OS [33]. Also, single-nucleotide polymorphisms of genes involved in homologous recombination, such as ATM and XRCC3, could be predictive and prognostic markers in metastatic colorectal cancer treated with TAS-102. The TAS-102-treated patients who carried any G allele in ATM rs609429 polymorphism had significantly longer OS and nonsignificantly longer PFS than carriers of the C/C variant. Also, the presence of any A allele in XRCC3 rs861539 polymorphism was shown to be associated with

Patients with ≥grade 2 chemotherapy-induced neutropenia (absolute neutrophil

median PFS (3.0 vs. 2.4 months) and OS (14.0 vs. 5.6 months) than those without the neutropenia. The presence of neutropenia at 1 month of TAS-102 treatment and higher baseline CEA levels were identified as independent predictors of OS [35]. In another study, grade 3 or 4 CIN-1-month was shown to be associated with longer PFS than grade 0-2 CIN-1-month (4.3 vs. 2.0 months). Moreover, G3 or G4 neutropenia during the first cycle of TAS-102 therapy turned out to be associated with significantly higher DCR (72.2%) than grade 0–2 neutropenia (72.2 vs. 30.9%) [36]. Finally, a significant association was found between baseline creatinine clearance rate of less than 57.1 mL/min prior to TAS-102 administration and the incidence of

Also, longer time elapsed since the onset of the first-line therapy to disease progression (more than 18 months) seems to be a predictor of better response to TAS-102 therapy. In one study, median PFS in TAS-102-treated patients with the time to progression exceeding 18 months was longer than in those who received the first-line therapy ≤18 months before (7 vs. 5 months) [38]. These findings were later

Kwakman et al. identified KRAS-wild type tumor, good performance status (PS0 or PS1), and normal serum levels of lactate dehydrogenase and alkaline phosphatase as independent predictors of better response to TAS-102 treatment. All these factors correlated positively with longer OS. Patients with KRAS-wild type tumors had longer median OS than those with KRAS-mutated malignancies (6.9 vs. 4.9 months), and median OS in persons with ECOG PS0–1 turned out to be longer than in individuals with worse performance status (5.9 vs. 3.2 months) [40]. However, in a meta-analysis involving the data of 1318 patients who received TAS-102, OS was significantly longer than in the study mentioned above, regardless the KRAS mutation status. Furthermore, the treatment response was not influenced by the number of metastatic sites (1, 2, or more) [39]. Also, in a subgroup analysis of patients participating in the RECOURSE trial, neither OS nor PFS correlated with

ECOG performance status (PS0 or PS1), the number of metastatic sites (1 or 2), and the time elapsed since the diagnosis of the first metastasis (18 months or longer) were identified as prognostic factors in the RECOURSE trial; however, none of those

Hematologic toxicities, including leukopenia, neutropenia, and anemia, and gastrointestinal toxicities, such as nausea/vomiting, diarrhea, and the loss of appetite, as well as fatigue of various grade, were the most frequent side effects observed in patients treated with trifluridine/tipiracil [40–45]. The most common grade 3/4 toxicity was myelosuppression (neutropenia, anemia, febrile neutropenia) [46]. In turn, cardiac ischemia seems to be one of the rarest AEs observed during the

) at 1 month (CIN-1-month) of TAS-102 therapy had longer

*DOI: http://dx.doi.org/10.5772/intechopen.85027*

count <1500/mm3

## *Further Therapeutic Options in Heavily Pretreated Colorectal Cancer Patients DOI: http://dx.doi.org/10.5772/intechopen.85027*

polymorphisms of genes involved in trifluridine metabolism and thymidine phosphorylase inhibitor excretion had significantly longer PFS and OS [33]. Also, single-nucleotide polymorphisms of genes involved in homologous recombination, such as ATM and XRCC3, could be predictive and prognostic markers in metastatic colorectal cancer treated with TAS-102. The TAS-102-treated patients who carried any G allele in ATM rs609429 polymorphism had significantly longer OS and nonsignificantly longer PFS than carriers of the C/C variant. Also, the presence of any A allele in XRCC3 rs861539 polymorphism was shown to be associated with significantly longer OS and PFS in TAS-102-treated patients [34].

Patients with ≥grade 2 chemotherapy-induced neutropenia (absolute neutrophil count <1500/mm3 ) at 1 month (CIN-1-month) of TAS-102 therapy had longer median PFS (3.0 vs. 2.4 months) and OS (14.0 vs. 5.6 months) than those without the neutropenia. The presence of neutropenia at 1 month of TAS-102 treatment and higher baseline CEA levels were identified as independent predictors of OS [35]. In another study, grade 3 or 4 CIN-1-month was shown to be associated with longer PFS than grade 0-2 CIN-1-month (4.3 vs. 2.0 months). Moreover, G3 or G4 neutropenia during the first cycle of TAS-102 therapy turned out to be associated with significantly higher DCR (72.2%) than grade 0–2 neutropenia (72.2 vs. 30.9%) [36]. Finally, a significant association was found between baseline creatinine clearance rate of less than 57.1 mL/min prior to TAS-102 administration and the incidence of G3 or G4 neutropenia after introduction of this agent [37].

Also, longer time elapsed since the onset of the first-line therapy to disease progression (more than 18 months) seems to be a predictor of better response to TAS-102 therapy. In one study, median PFS in TAS-102-treated patients with the time to progression exceeding 18 months was longer than in those who received the first-line therapy ≤18 months before (7 vs. 5 months) [38]. These findings were later confirmed in a meta-analysis conducted by Chen et al. [39].

Kwakman et al. identified KRAS-wild type tumor, good performance status (PS0 or PS1), and normal serum levels of lactate dehydrogenase and alkaline phosphatase as independent predictors of better response to TAS-102 treatment. All these factors correlated positively with longer OS. Patients with KRAS-wild type tumors had longer median OS than those with KRAS-mutated malignancies (6.9 vs. 4.9 months), and median OS in persons with ECOG PS0–1 turned out to be longer than in individuals with worse performance status (5.9 vs. 3.2 months) [40]. However, in a meta-analysis involving the data of 1318 patients who received TAS-102, OS was significantly longer than in the study mentioned above, regardless the KRAS mutation status. Furthermore, the treatment response was not influenced by the number of metastatic sites (1, 2, or more) [39]. Also, in a subgroup analysis of patients participating in the RECOURSE trial, neither OS nor PFS correlated with the KRAS status, as well as with patients' age and ethnicity [31].

ECOG performance status (PS0 or PS1), the number of metastatic sites (1 or 2), and the time elapsed since the diagnosis of the first metastasis (18 months or longer) were identified as prognostic factors in the RECOURSE trial; however, none of those factors turned out to be a predictor of therapeutic response [30].


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

## **Shuman Sun, Zhiming Li, Huiguo Zhang , Haijun Jiang and Xijian Hu \***

College of Mathematics and System Science, Xinjiang University, Urumqi 830046, China<EMAIL_ADDRESS>(S.S<EMAIL_ADDRESS>(Z.L<EMAIL_ADDRESS>(H.Z<EMAIL_ADDRESS>(H.J.)

**\*** Correspondence<EMAIL_ADDRESS>
Received: 21 September 2020; Accepted: 26 October 2020; Published: 29 October 2020

**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.

**Keywords:** HIV/AIDS epidemic; regression model; Newton–Raphson procedure; Fisher scoring algorithm; time series

#### **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–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. Sections 2.2 and 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–15]. Suppose that *Y* = (*Y*1,*Y*2, ... ,*Yn*) is the response variable, where *Yi*(*i* = 1, . . . , *n*) are independent. The density function is

$$f(y\_i; \theta\_i, \phi\_i) = \exp\left(\frac{y\_i \theta\_i - b(\theta\_i)}{a(\phi\_i)} + c(y\_i, \phi\_i)\right),$$

where *a*(·), *b*(·), *c*(·, ·) are known functions, and *θi*, *φi* are unknown parameters for *i* = 1, 2, ... , *n*. Denote *μi* = *E*(*Yi*), and *g*(*μi*) = ln(*μi*) is a link function. Let *Xij* be explanatory variables for the *i*th observation in the *j*th variable. Then, the Poisson regression (PR) model is given by

$$\log(\mu\_i) \triangleq \eta\_i = \sum\_{j=1}^p \pounds\_j X\_{ij\prime} \tag{1}$$

where *i* = 1, 2, . . . , *n*, and *βj*(*j* = 1, 2, . . . , *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–20] below

$$\log(\mu\_i) \triangleq \eta\_i = \sum\_{j=1}^p \beta\_j (u\_{i\prime} v\_{i\prime}) X\_{ij\prime} \tag{2}$$

where (*ui*, *vi*)(*i* = 1, 2, ... , *n*) are the geographical locations, and *βj*(*ui*, *vi*)(*j* = 1, 2, ... , *p*) are unknown parameters at the position (*ui*, *vi*).

#### *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–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:

$$\log(\mu\_i) \triangleq \eta\_i = \rho \sum\_{k=1}^n c\_{ik} \eta\_k + \sum\_{j=1}^p \beta\_j (\mu\_{i\prime} \upsilon\_i) X\_{ij\prime} \tag{3}$$

where *ρ* is a scalar autoregressive parameter, and *cik*(*i*, *k* = 1, 2, ... , *n*) is the adjacency relation between the *i*th and *k*th locations. Let *ci* be the number of regions adjacent to the *i*th position. If the *k*th position is next to the *i*th's, then *cik* = 1/*ci* . Otherwise, *cik* = 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 *Yik*(*i* = 1, 2, ... , *nk*, *k* = 1, 2, ... , *T*) be the independent variables of *Y* in the *i*th position and the *k*th time. The density function can be defined as follows:

$$f(y\_{ik}; \theta\_{ik}, \phi\_{ik}) = \exp\left(\frac{y\_{ik}\theta\_{ik} - b(\theta\_{ik})}{a(\phi\_{ik})} + \mathfrak{c}(y\_{ik}, \phi\_{ik})\right),$$

where the parameters are similar to Section 2.1. Denote *μik* = *E*(*Yik*), and *g*(*μik*) = ln(*μik*). Let *Xijk*(*j* = 1, 2, . . . , *p*) be the *j*th explanatory variable. The GTWPAR model is expressed by

$$\log(\mu\_{ik}) \triangleq \eta\_{ik} = \rho \sum\_{m=1}^{T} \sum\_{l=1}^{n\_k} c\_{lm}^{(\bar{k})} \eta\_{lm} + \sum\_{j=1}^{p} \sum\_{k=1}^{T} \beta\_{jk} (\mu\_{i\bar{k}\lambda}, \upsilon\_{i\bar{k}\lambda}, t\_k) \chi\_{i\bar{j}\lambda}, \tag{4}$$

where {*βjk*(*uik*, *vik*, *tk*)} is a set of unknown parameters at the *i*th position in the *k*th time, and *c* (*ik*) *lm* is the adjacent relation between the location (*uik*, *vik*, *tk*) and (*ulm*, *vlm*, *tm*). Following the work of [31], the spatiotemporal distance between the locations (*uik*, *vik*, *tk*) and (*ulm*, *vlm*, *tm*) can be defined as

$$d\_{lm}^{(ik)} = \sqrt{\lambda[(\mu\_{ik} - \mu\_{lm})^2 + (\upsilon\_{ik} - \upsilon\_{lm})^2] + \mu(t\_k - t\_m)^2},$$

where *μ* and *λ* are used to balance spatiotemporal distances. Suppose that

$$
\mathcal{c}\_{lm}^{(ik)} = \begin{cases}
1/c^{ik}, & 0 < d\_{lm}^{(ik)} < d\_{\prime i} \\
0, & \text{otherwise},
\end{cases}
$$

where *d* is a constant and satisfies min{*d* (*ik*) *lm* } < *d* < max{*d* (*ik*) *lm* }.

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

$$
\eta = \rho \mathbf{C} \eta + \mathbf{B}' \mathbf{X}',
$$

where *η* = (*η*11, ··· , *ηn*11, *η*12, ··· , *ηn*22, ··· , *η*1*T*, ··· , *ηnTT*) , *C* = (*c* (*ik*) *lm* ), *X* = (*Xijk*) and **B** = (*βjk*(*uik*, *vik*, *tk*)). For convenience, define *ηK* as the *K*th element of *η*; *CIK* and *XIK* are the *I*th row and the *K*th column of the matrices *C* and *X*, respectively. The detailed expressions of *C*, *X* and **B** are given in Appendix A.1.

**Remark 1.** *For the GTWPAR model (4), if ρ* = 0 *and βjk*(*uik*, *vik*, *tk*) *is independent of the spatiotemporal effect, the model is a PR model. If ρ* = 0 *and βjk*(*uik*, *vik*, *tk*) *is dependent on spatial effect but independent of temporal effect, the model becomes GWPR model. If ρ* = 0 *and βjk*(*uik*, *vik*, *tk*) *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 (*uik*, *vik*, *tk*)(*i* = 1, 2, ... , *nk*, *k* = 1, 2, ... , *T*) be any point in the studied spatiotemporal region. We fix a point (*u*00, *v*00, *t*0) and assume that *βjk*(*uik*, *vik*, *tk*) ≈ *βj*0(*u*00, *v*00, *t*0)(*j* = 1, 2, . . . , *p*). Then, the model (4) can be rewritten by

$$\eta\_{ik} = \mathcal{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} \mathcal{\beta}\_{j0}(\mu\_{00}, v\_{00}, t\_0) X\_{ijk} \tag{5}$$

Denote *β*(*u*00, *v*00, *t*0)=(*β*10, ... , *βp*0) , **X** = diag(*Xi*.) and *Xi*· = (*Xi*1, ... , *Xip*). The corresponding matrix form can be represented as *η* = *ρCη* + *β* (*u*00, *v*00, *t*0)**X** .

#### *3.1. Estimation of Parameter Vector β*

For the fixed point (*u*00, *v*00, *t*0), we define a spatiotemporal distance *d* (0) *ik* from this point to (*uik*, *vik*, *tk*) as *d* (0) *ik* = *λ*[(*u*00 − *uik*)2 + (*v*00 − *vik*)2] + *μ*(*t*0 − *tk*)2. The Gauss kernel function of these two points can be written by

$$\begin{split} w\_{ik}(\boldsymbol{\mu\_{00}}, \boldsymbol{\upsilon\_{00}}, t\_{0}) &= \quad \frac{1}{\sqrt{2\pi}} \exp\left\{-\frac{1}{2} \left(\frac{d\_{ik}^{(0)}}{h\_{ST}}\right)^{2}\right\} \\ &= \quad \frac{1}{\sqrt{2\pi}} \exp\left\{-\frac{1}{2} \frac{\lambda\left[(\boldsymbol{\mu\_{00}} - \boldsymbol{\mu\_{ik}})^{2} + (\boldsymbol{\upsilon\_{00}} - \boldsymbol{\upsilon\_{ik}})^{2}\right] + \mu(t\_{0} - t\_{k})^{2}}{h\_{ST}^{2}}\right\} \\ &= \quad \frac{1}{\sqrt{2\pi}} \exp\left\{-\frac{1}{2} \left(\frac{(\boldsymbol{\mu\_{00}} - \boldsymbol{\mu\_{ik}})^{2} + (\boldsymbol{\upsilon\_{00}} - \boldsymbol{\upsilon\_{ik}})^{2}}{h\_{S}^{2}} + \frac{(t\_{0} - t\_{k})^{2}}{\tau h\_{S}^{2}}\right)\right\}, \end{split}$$

where *hST* and *hS* are the space-time bandwidth and space bandwidth, respectively. Meanwhile, we have *h*2 *ST* = *λh*2 *S*, and *τ* = *λ*/*μ* is a spatiotemporal factor. Without loss of generality, let *λ* = 1. Then, the weighted maximum likelihood of *Yik*(*i* = 1, 2, . . . , *nk*, *k* = 1, 2, . . . , *T*) at the point (*u*00, *v*00, *t*0) is

$$L(\boldsymbol{\beta}\_{10\prime}, \boldsymbol{\beta}\_{20\prime}, \dots, \boldsymbol{\beta}\_{p0}) = \prod\_{k=1}^{T} \prod\_{i=1}^{n\_k} f(\boldsymbol{y}\_{ik\prime} \boldsymbol{\theta}\_{ik\prime} \boldsymbol{\phi}\_{ik}) w\_{ik}(\boldsymbol{u}\_{00\prime}, \boldsymbol{v}\_{00\prime}, t\_0),$$

where *f*(*yik*; *θik*, *φik*) is the density function. The log-likelihood can be obtained as follows:

$$\mathbf{L}\_1(\boldsymbol{\beta}(\boldsymbol{\mu\_{00}}, \boldsymbol{\nu\_{00}}, t\_0)) = \sum\_{k=1}^{T} \sum\_{l=1}^{n\_k} \left( \frac{y\_{lk}\theta\_{lk} - b(\theta\_{lk})}{a(\phi\_{lk})} + c(y\_{lk}, \phi\_{lk}) \right) w\_{lk}(\boldsymbol{\mu\_{00}}, \boldsymbol{\nu\_{00}}, t\_0) \dots$$

Note that *c*(*yik*, *φik*) = − ln(*yik*!), *b*(*θik*) = *μik* = exp(*θik*), and *a*(*φik*) = *φik* = 1. Thus, *E*(*Yik*) = *b* (*θik*) = exp(*θik*) = *μik*, *Var*(*Yik*) = *b*(*θik*)*a*(*φik*) = exp(*θik*) = *μik*. Differentiating **L**1 with respect to *β*(*u*00, *v*00, *t*0) yields

$$\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}\phi} \frac{\partial \theta\_{lk}}{\partial \beta\_{r0}} \right) w\_{ik} (u\_{00\prime} v\_{00\prime} t\_0) = 0,\tag{6}$$

where *βr*0 = *βr*(*u*00, *v*00, *t*0)(*r* = 1, 2, ··· , *p*), and

$$\frac{\partial \theta\_{ik}}{\partial \beta\_{r0}} = \left(\frac{\partial \mu\_{ik}}{\partial \theta\_{ik}}\right)^{-1} \frac{\partial \mu\_{ik}}{\partial \mathcal{S}(\mu\_{ik})} \frac{\partial \mathcal{g}(\mu\_{ik})}{\partial \beta\_{r0}} = \frac{1}{b''(\theta\_{ik})} \frac{1}{\mathcal{g}'(\mu\_{ik})} \frac{\partial \eta\_{ik}}{\partial \beta\_{r0}}.$$

For convenience, let *N* = ∑*N k*=1 *nk* and *W* = (*wik*(*u*00, *v*00, *t*0))*N*×*N*. Denote *A* = (*IN* − *ρC*)−1, *Y* = (*Y*11, ··· ,*Yn*11, ··· ,*Y*1*T*, ··· ,*YnTT*) , *μ* = (*μ*11, ··· , *μn*11, ··· , *μ*1*T*, ··· , *μnTT*) , *θ* = (*θ*11, ··· , *θn*11, ··· , *θ*1*T*, ··· , *θnTT*) , *φ* = (*φ*11, ··· , *φn*11, ··· , *φ*1*T*, ··· , *φnTT*) . Suppose that *YK*, *μK*, *θK* and *φK* are the *K*th elements of *Y*, *μ*, *θ* and *φ*, respectively. Then, we take the derivative of the model (5) with respect to *βr*0, and obtain

$$\frac{\partial \eta\_l}{\partial \beta\_{r0}} = \sum\_{h=1}^{N} A\_{lh} X\_{hr} = A\_l X\_{r\prime} \quad l = 1, 2, \dots, N.$$

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

$$\frac{\partial \mathbf{L}\_1}{\partial \beta\_{r0}} = \frac{1}{\Phi} \sum\_{l=1}^{N} T\_l A\_{l\cdot\cdot} X\_{\cdot r} (\mathbf{Y}\_l - \mu\_l) \mathbf{g}^{\prime}(\mu\_l) \mathcal{W}\_l(\mu\_{00\prime}, \upsilon\_{00\prime}, t\_0) = 0.$$

However, there is not a close-form solution for *β*(*u*00, *v*00, *t*0). The Newton–Raphson procedure and Fisher scoring algorithm are used to get the estimation of *β*. The iterative formula is expressed as

$$\begin{array}{rcl}\beta^{(m+1)}(\boldsymbol{u}\_{00},\boldsymbol{v}\_{00},t\_{0})&=&\beta^{(m)}(\boldsymbol{u}\_{00},\boldsymbol{v}\_{00},t\_{0})+\mathbf{I}^{-1}(\beta^{(m)}(\boldsymbol{u}\_{00},\boldsymbol{v}\_{00},t\_{0}))S(\beta^{(m)}(\boldsymbol{u}\_{00},\boldsymbol{v}\_{00},t\_{0})) \\&=&((\boldsymbol{A}^{(m)}\boldsymbol{X})^{\prime}\boldsymbol{T}^{(m)}\boldsymbol{W}(\boldsymbol{u}\_{00},\boldsymbol{v}\_{00},t\_{0})(\boldsymbol{A}^{(m)}\boldsymbol{X}))^{-1} \\&\times(\boldsymbol{A}^{(m)}\boldsymbol{X})^{\prime}\boldsymbol{T}^{(m)}\boldsymbol{W}(\boldsymbol{u}\_{00},\boldsymbol{v}\_{00},t\_{0})\boldsymbol{Z}^{(m)}),\end{array}\tag{7}$$

where the Fisher information matrix **I**(*β*) = *E*(*I*(*β*)), and

$$S(\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}}, \dots, \frac{\partial \mathbf{L}\_1}{\partial \beta\_{p0}}\right)'$$

is the scalar vector. The detail process is provided in Appendix A.2. For the fixed point (*uik*, *vik*, *tk*)(*i* = 1, 2, . . . , *nk*; *k* = 1, 2, . . . , *T*), *β*ˆ *jk*(*uik*, *vik*, *tk*) can be obtained by (7).

**Remark 2.** *The estimations β*ˆ(*uik*, *vik*, *tk*)(*i*, *l* = 1, 2, ... , *nk*, *k*, *m* = 1, 2, ... , *T*) *are related to the temporal and spatial effects in the GTWPAR model. If m* = *k, wik*(*ulm*, *vlm*, *tm*) = 0 *and c* (*ik*) *lm* = 0*, then β*ˆ(*uik*, *vik*, *tk*) = *β*ˆ(*ui*, *vi*) *correspond to the parameter estimations of the GWPAR model. If wik*(*ulm*, *vlm*, *tm*) = 0(*m* = *k*) *and C* = **0***, they are the estimations of the GWPR model. If W* = **0** *and C* = **0***, then β*ˆ(*uik*, *vik*, *tk*) = *β*ˆ *are the global estimation values of the PR model.*

#### *3.2. Estimation of Parameter ρ*

Based on the density function, the log-likelihood function of *ρ* is

$$\mathbf{L}\_2(\rho) = \sum\_{k=1}^{T} \sum\_{i=1}^{n\_k} \left( \frac{y\_{ik}\theta\_{ik} - b(\theta\_{ik})}{a(\phi\_{ik})} + c(y\_{ik}, \phi\_{ik}) \right).$$

Differentiating **L**2(*ρ*) with respect to *ρ*, we have

$$\frac{\partial \mathbf{L}\_2}{\partial \rho} = \sum\_{k=1}^{T} \sum\_{i=1}^{n\_k} \left( \frac{y\_{ik} - \mu\_{ik}}{a\_{ik} \phi} \frac{\partial \theta\_{ik}}{\partial \rho} \right) = 0,\tag{8}$$

where *dθik dρ* = 1 *b*(*θik* )*g*(*μik* ) *dηik dρ* . Then, we take the derivative of the model (5) with respect to *ρ* as follows:

$$\frac{d\eta\_l}{d\rho} = \frac{d\lg(\mu\_l)}{d\rho} = \sum\_{h=1}^N A\_h \mathcal{C}\_h \eta\_h.$$

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{(\mathbf{Y}\_l - \mu\_l)\sum\_{l=1}^{N} A\_{l\cdot} \mathbf{C}\_{\cdot l\cdot} \eta\_{l\cdot}}{a\_l \phi V(\mu\_l) \mathbf{g}'(\mu\_l)} = 0.1$$

According to the Newton–Raphson procedure and Fisher scoring algorithm, the iterative formula of *ρ*ˆ(*m*+1) is

$$\begin{array}{rcl} \boldsymbol{\beta}^{(m+1)} &=& \boldsymbol{\beta}^{(m)} + \boldsymbol{\mathcal{Z}}^{-1}(\boldsymbol{\beta}^{(m)}) \mathbf{S}(\boldsymbol{\beta}^{(m)}) \\ &=& \boldsymbol{\beta}^{(m)} + ((\boldsymbol{\mathcal{A}}^{(m)} \mathbf{C} \boldsymbol{\eta}^{(m)})' \boldsymbol{T}^{(m)} (\boldsymbol{\mathcal{A}}^{(m)} \mathbf{C} \boldsymbol{\eta}^{(m)}))^{-1} \\ & \times (\boldsymbol{\mathcal{A}}^{(m)} \mathbf{C} \boldsymbol{\eta}^{(m)})' \boldsymbol{T}^{(m)} (\boldsymbol{\mathcal{Z}}^{(m)} - \boldsymbol{\eta}^{(m)}), \end{array} \tag{9}$$

where the scalar vector *S*(*ρ*ˆ(*m*)) = 1 *φ* (*ACη*) *T*(*Z* − *η*) and the Fisher information matrix I(*ρ*) = 1 *φ* (*ACη*) *T*(*ACη*). The calculation process of the scalar vector *S*(*ρ*ˆ(*m*)) and the information matrix 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 R2, 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.

**Figure 1.** Spatiotemporal HIV/AIDS incidence of ten countries, 2011–2016.

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 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.

**Figure 2.** Distributions of HIV/AIDS cases, GDP per capita and population density of ten countries, 2011–2016.

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

Next, we standardized the two socioeconomic factors. The multiplex collinear test [32] was performed by the condition number *k* = √*λ*max/*λ*min = 1.804(≤ 15) (*λ* 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 *μik*, *rik* and *Pik* be the annual HIV/AIDS cases (Unit: 1/1000 people), incidence (Unit: 1/100) and total population (unit: 100,000 people) in the *k*th year of the *i*th region, respectively. Denote *g*(*μik*) = *ηik* = ln *μik* = ln *rik* + ln *Pik*(*μik* = *rikPik*, *i* = 1, 2, ... , 10 and *k* = 1, 2, ... , 6). Let *Xi*1*k* and *Xi*2*k* be the GDP per capita and population density in the *i*th region at the *k*th year, respectively. The PR model is written by

$$\lg(\mu\_{ik}) = \beta\_0 + \beta\_1 X\_{i1k} + \beta\_2 X\_{i2k}, \quad i = 1, 2, \dots, 10, k = 1, 2, \dots, 6,\tag{10}$$

where *βj*(*j* = 0, 1, 2) are unknown constants. The GWPR model is introduced as

$$\mathcal{S}\_{\mathcal{S}}(\mu\_{ik}) = \beta\_0(u\_{ik}, v\_{ik}) + \beta\_1(u\_{ik}, v\_{ik})X\_{i1k} + \beta\_2(u\_{ik}, v\_{ik})X\_{i2k}, \quad i = 1, 2, \dots, 10,\tag{11}$$

where *k* is a fixed constant taken from {1, 2, ... , 6}, and *βj*(*uik*, *vik*) are unknown spatial parameters for the *i*th country (*uik*, *vik*) in the *k*th year. Let *ρ* be a scalar autoregressive parameter, and *cil* be a constant that represents an adjacency relation. The GWPAR model is

$$\log(\mu\_{i\bar{k}}) = \rho \sum\_{l=1}^{n} c\_{l\bar{l}} \eta\_{i\bar{k}} + \beta \mathbf{o}(u\_{i\bar{k}\succ i} \upsilon\_{i\bar{k}}) + \beta \mathbf{1}(u\_{i\bar{k}\succ i} \upsilon\_{i\bar{k}}) X\_{i1\bar{k}} + \beta \mathbf{2}(u\_{i\bar{k}\succ i} \upsilon\_{i\bar{k}}) X\_{i2\bar{k}} \tag{12}$$

where *n* = 10, *k* is a fixed constant, and *βj*(*uik*, *vik*) are defined as above. Let *c* (*ik*) *lm* be a spatiotemporal adjacency relation, and *βjk*(*uik*, *vik*, *tk*)(*k* = 1, 2, ... , 6) be unknown spatiotemporal parameters in the *i*th country (*uik*, *vik*) in the *k*th year. The GTWPAR model is established as follows:

$$\begin{aligned} \text{g}(\mu\_{ik}) &= \rho \sum\_{m=1}^{T} \sum\_{l=1}^{n\_k} c\_{lm}^{(ik)} \eta\_{lm} + \beta\_{0k} (\mu\_{ik}, \upsilon\_{ik}, t\_k) \\ &+ \beta\_{1k} (\mu\_{ik}, \upsilon\_{ik}, t\_k) X\_{i1k} + \beta\_{2k} (\mu\_{ik}, \upsilon\_{ik}, t\_k) X\_{i2k} \end{aligned} \tag{13}$$

where *T* = 6; *nk* = 10 for every *k* years; and *ρ* 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 *ρ* = 0.267, 0.269, 0.263, 0.264, 0.264, 0.264. For the GTWPAR model, we chose *hs* = 1.1316, 0.9737, 0.8947, 0.9211, 0.7842, 0.8789 and *τ* = 0.1, where *τ*(> 0) is a balanced parameter. The coefficient estimations and scalar autoregressive parameter *ρ* = 0.126 can be obtained by Algorithm IV. The quantile and mean values of coefficient estimations and response variables are shown in Table 2. 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.

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

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.

*Entropy* **2020**, *22*, 1230

**Figure 3.** The spatial distribution of the average coefficient estimations in four models.

**Figure 4.** The average MSE of response variables.

These four indicators can effectively compare the performances of the proposed models (Table 3). The calculation formulas of R2, AICc, D and MSE are given in Appendix A.5. The coefficient of determination R2 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.

**Table 3.** The comparison of the four models.

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.

**Figure 5.** The mean values and 95% confidence intervals of coefficient estimations.

#### **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.

**Author Contributions:** Conceptualization, Z.L., H.Z. and X.H.; methodology, S.S., H.Z. and X.H.; software, S.S.; supervision, Z.L. and H.J.; visualization, Z.L.; writing—original draft, S.S.; writing—review & editing, S.S., Z.L. and H.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Natural Science Foundation of Xinjiang Uygur Autonomous Region (XJEDU2017M001, 2018Q011), and the National Natural Science Foundation of China (U1703237, 11661076, 12061070).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

#### **Appendix A. Detailed Processes**

*Appendix A.1.    *XnT*1*T XnT*2*T* ··· *XnT pT* ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ ,

where *CIK*, *XIK* are respectively the *I*th row and the *K*th column of the matrices *C* and *X*. Moreover,

**B** = (*β*11(*u*11, *v*11, *t*1), ··· , *βp*1(*u*11, *v*11, *t*1), ··· , *β*11(*un*11, *vn*11, *t*1), ··· , *βp*1(*un*11, *vn*11, *t*1), ··· , *β*1*T*(*u*1*T*, *v*1*T*, *tT*), ··· , *βpT*(*u*1*T*, *v*1*T*, *tT*), ··· , *β*1*T*(*unT*1, *vnT*1, *t*1), ··· , *βpT*(*unT*1, *vnT*1, *t*1)) .

*Appendix A.2. Formula and Information Matrix of β*(*u*, *v*, *t*)

(1) For the matrix form of model (5), we obtain

$$\begin{cases} \begin{aligned} &(1 - \rho \mathbb{C}\_{11})\eta\_1 - \rho \mathbb{C}\_{12}\eta\_2 - \cdots - \rho \mathbb{C}\_{1N}\eta\_N = \sum\_{j=1}^p \beta\_{j0} \mathbb{X}\_{1j}, \\ & - \rho \mathbb{C}\_{21}\eta\_1 + (1 - \rho \mathbb{C}\_{22})\eta\_2 - \cdots - \rho \mathbb{C}\_{2N}\eta\_N = \sum\_{j=1}^p \beta\_{j0} \mathbb{X}\_{2j}, \\ & \vdots \\ & \vdots \\ & - \rho \mathbb{C}\_{N1}\eta\_1 - \rho \mathbb{C}\_{N2}\eta\_2 - \cdots + (1 - \rho \mathbb{C}\_{NN})\eta\_N = \sum\_{j=1}^p \beta\_{j0} \mathbb{X}\_{Nj}. \end{aligned} \end{cases}$$

Differentiating the above equations with *βr*0 yields

$$\begin{cases} (1 - \rho \mathbf{C}\_{11}) \frac{\partial \eta\_{1}}{\partial \beta\_{r0}} - \rho \mathbf{C}\_{12} \frac{\partial \eta\_{2}}{\partial \beta\_{r0}} - \dots - \rho \mathbf{C}\_{1N} \frac{\partial \eta\_{N}}{\partial \beta\_{r0}} = X\_{1r}, \\ -\rho \mathbf{C}\_{21} \frac{\partial \eta\_{1}}{\partial \beta\_{r0}} + (1 - \rho \mathbf{C}\_{22}) \frac{\partial \eta\_{2}}{\partial \beta\_{r0}} - \dots - \rho \mathbf{C}\_{2N} \frac{\partial \eta\_{N}}{\partial \beta\_{r0}} = X\_{2r}, \\ \vdots \\ -\rho \mathbf{C}\_{N1} \frac{\partial \eta\_{1}}{\partial \beta\_{r0}} - \rho \mathbf{C}\_{N2} \frac{\partial \eta\_{2}}{\partial \beta\_{r0}} - \dots + (1 - \rho \mathbf{C}\_{NN}) \frac{\partial \eta\_{N}}{\partial \beta\_{r0}} = X\_{Nr}. \end{cases}$$

Then,

$$
\begin{pmatrix}
1 - \rho \mathbb{C}\_{11} & -\rho \mathbb{C}\_{12} & \cdots & -\rho \mathbb{C}\_{1N} \\
\vdots & \vdots & \vdots & \vdots \\
\end{pmatrix}
\begin{pmatrix}
\frac{\partial \eta\_{1}}{\partial \mathbb{B}\_{r0}} \\
\frac{\partial \eta\_{2}}{\partial \mathbb{B}\_{r0}} \\
\vdots \\
\frac{\partial \eta\_{N}}{\partial \mathbb{B}\_{r0}}
\end{pmatrix} =
\begin{pmatrix}
X\_{1r} \\
X\_{2r} \\
\vdots \\
X\_{Nr}
\end{pmatrix}.
$$

Denote *∂η ∂βr*0 = ( *∂η*1 *∂βr*0 , *∂η*2 *∂βr*0 ,..., *∂ηN ∂βr*0 ) and

$$A = \begin{pmatrix} 1 - \rho \mathbb{C}\_{11} & -\rho \mathbb{C}\_{12} & \cdots & -\rho \mathbb{C}\_{1N} \\ -\rho \mathbb{C}\_{21} & 1 - \rho \mathbb{C}\_{22} & \cdots & -\rho \mathbb{C}\_{2N} \\ \vdots & \vdots & \vdots & \vdots \\ -\rho \mathbb{C}\_{N1} & -\rho \mathbb{C}\_{N2} & \cdots & 1 - \rho \mathbb{C}\_{NN} \end{pmatrix}^{-1},$$

$$X\_{I} = (X\_{1r}, X\_{2r}, \dots, X\_{Nr})', A\_{I} = (A\_{I1}, A\_{I2}, \dots, A\_{IN})'.$$

Thus, *∂η ∂βr*0 = *AX*·*r*, that is

$$\frac{\partial \eta\_l}{\partial \beta\_{r0}} = \sum\_{h=1}^{N} A\_{lh} X\_{hr} = A\_l X\_{,r\_l} \quad l = 1, 2, \dots, N.$$

(2) The element *Irb* of *I*(*β*) satisfies

*Irb*(*β*) = − *∂*2**L**1 *∂βb*0*∂βr*0 = − *∂ ∂βb*0 *N* ∑ *l*=1 *Yl* − *μl alφ Al*·*X*·*r V*(*μl*)*g*(*μl*) *Wl*(*u*00, *v*00, *t*0) = − *N* ∑ *l*=1 *Yl* − *μl alφ ∂ ∂βb*0 *Al*·*X*·*r V*(*μl*)*g*(*μl*) *Wl*(*u*00, *v*00, *t*0) − *N* ∑ *l*=1 *Al*·*X*·*r V*(*μl*)*g*(*μl*) *∂ ∂βb*0 *Yl* − *μl alφ Wl*(*u*00, *v*00, *t*0) = − *N* ∑ *l*=1 *Yl* − *μl alφ ∂ ∂βb*0 *Al*·*X*·*r V*(*μl*)*g*(*μl*) *Wl*(*u*00, *v*00, *t*0) + *N* ∑ *l*=1 *Al*·*X*·*rAl*·*X*·*b alφV*(*μl*)(*g*(*μl*))2 *Wl*(*u*00, *v*00, *t*0).

The Fisher information matrix is

$$\begin{split} \mathbf{I}(\boldsymbol{\beta}) = \mathbf{E}(I(\boldsymbol{\beta})) = \mathbf{E}((I\_{rb}(\boldsymbol{\beta}))\_{p \times p}) &= \sum\_{l=1}^{N} \frac{A\_{l} \mathbf{X}\_{r} A\_{l} \mathbf{X}\_{\cdot b}}{a\_{l} \boldsymbol{\Phi} V(\mu\_{l}) (\mathbf{g}^{r}(\mu\_{l}))^{2}} \mathbf{W}\_{l}(\mu\_{00}, \boldsymbol{\nu}\_{00}, t\_{0}) \\ &= \frac{1}{\boldsymbol{\Phi}} \sum\_{l=1}^{N} T\_{l} A\_{l} \mathbf{X}\_{r} A\_{l} \mathbf{X}\_{b} \mathbf{W}\_{l}(\mu\_{00}, \boldsymbol{\nu}\_{00}, t\_{0}). \end{split}$$

*Appendix A.3. Formula and Information Matrix of ρ*

(1) Differentiating (4) with *ρ* yields

$$\begin{cases} (1 - \rho \mathbb{C}\_{11}) \frac{d\eta\_{1}}{d\rho} - \rho \mathbb{C}\_{12} \frac{d\eta\_{2}}{d\rho} - \dots - \rho \mathbb{C}\_{1N} \frac{d\eta\_{N}}{d\rho} = \mathbb{C}\_{11}\eta\_{1} + \mathbb{C}\_{12}\eta\_{2} + \dots + \mathbb{C}\_{1N}\eta\_{N}, \\ -\rho \mathbb{C}\_{21} \frac{d\eta\_{1}}{d\rho} + (1 - \rho \mathbb{C}\_{22}) \frac{d\eta\_{2}}{d\rho} - \dots - \rho \mathbb{C}\_{2N} \frac{d\eta\_{N}}{d\rho} = \mathbb{C}\_{21}\eta\_{1} + \mathbb{C}\_{22}\eta\_{2} + \dots + \mathbb{C}\_{2N}\eta\_{N}, \\ \vdots \\ -\rho \mathbb{C}\_{N1} \frac{d\eta\_{1}}{d\rho} - \rho \mathbb{C}\_{N2} \frac{d\eta\_{2}}{d\rho} - \dots + (1 - \rho \mathbb{C}\_{NN}) \frac{d\eta\_{N}}{d\rho} = \mathbb{C}\_{N1}\eta\_{1} + \mathbb{C}\_{N2}\eta\_{2} + \dots + \mathbb{C}\_{NN}\eta\_{N}. \end{cases}$$

That is to say

$$
\begin{pmatrix}
1 - \rho \mathbb{C}\_{11} & -\rho \mathbb{C}\_{12} & \cdots & -\rho \mathbb{C}\_{1N} \\
\vdots & \vdots & \vdots & \vdots \\
\end{pmatrix}
\begin{pmatrix}
\frac{d\eta\_{1}}{d\rho} \\
\frac{d\eta\_{2}}{d\rho} \\
\vdots \\
\frac{d\eta\_{N}}{d\rho}
\end{pmatrix} = \begin{pmatrix}
\mathbb{C}\_{11} & \mathbb{C}\_{12} & \cdots & \mathbb{C}\_{1N} \\
\mathbb{C}\_{21} & \mathbb{C}\_{22} & \cdots & \mathbb{C}\_{2N} \\
\vdots & \vdots & \vdots & \vdots \\
\mathbb{C}\_{N1} & \mathbb{C}\_{N2} & \cdots & \mathbb{C}\_{NN}
\end{pmatrix} \begin{pmatrix}
\eta\_{1} \\
\eta\_{2} \\
\vdots \\
\eta\_{N}
\end{pmatrix}.
$$

Let *dη dρ* = ( *dη*1 *dρ* ,..., *dηN dρ* ) . Then, *dη dρ* = *ACη*, where

$$\frac{d\eta\_l}{d\rho} = \frac{d\mathcal{g}(\mu\_l)}{d\rho} = \sum\_{h=1}^N A\_l \mathbb{C}\_{\text{-}h} \eta\_{h\prime} \quad l = 1, 2, \cdots, N.$$

(2) The scalar vector of *ρ* is

$$\begin{aligned} S(\rho) &= \quad \frac{d\mathbf{L}\_2}{d\rho} = \frac{1}{\Phi} \sum\_{l=1}^N T\_l \left( \sum\_{h=1}^N A\_{l\cdot} \mathbb{C}\_{\cdot h} \eta\_h \right) (Y\_l - \mu\_l) \mathbf{g}'(\mu\_l), \\ &= \quad \frac{1}{\Phi} (AC\eta)' T(Z - \eta). \end{aligned}$$

The information matrix is

$$\begin{array}{rcl} I(\rho) &=& -\frac{\partial^2 \mathcal{L}\_2}{\partial \rho^2} = -\frac{d}{d\rho} \sum\_{l=1}^N \left( \frac{\mathcal{Y}\_l - \mu\_l}{a\_l \phi} \right) \left( \frac{\sum\_{h=1}^N A\_l \mathcal{C}\_{\cdot h} \eta\_h}{V(\mu\_l) \mathcal{G}'(\mu\_l)} \right) \\ &=& -\sum\_{l=1}^N \left( \frac{\mathcal{Y}\_l - \mu\_l}{a\_l \phi} \right) \frac{d}{d\rho} \left( \frac{\sum\_{h=1}^N A\_l \mathcal{C}\_{\cdot h} \eta\_h}{V(\mu\_l) \mathcal{G}'(\mu\_l)} \right) \\ & - \sum\_{l=1}^N \left( \frac{\sum\_{h=1}^N A\_l \mathcal{C}\_{\cdot h} \eta\_h}{V(\mu\_l) \mathcal{G}'(\mu\_l)} \right) \frac{d}{d\rho} \left( \frac{\mathcal{Y}\_l - \mu\_l}{a\_l \phi} \right) . \end{array}$$

where

$$\frac{d\mu\_l}{d\rho} = \frac{d\mu\_l}{d\lg(\mu\_l)} \frac{d\lg(\mu\_l)}{d\rho} = \frac{1}{\lg'(\mu\_l)} \sum\_{h=1}^N A\_{l^\cdot} \mathbb{C}\_{h\eta\eta\_h}.$$

Thus,

$$I(\rho) = -\sum\_{l=1}^{N} \left( \frac{\Upsilon\_l - \mu\_l}{a\_l \phi} \right) \frac{d}{d\rho} \left( \frac{\sum\_{h=1}^{N} A\_{l\cdot} \mathcal{C}\_{\cdot h} \eta\_h}{V(\mu\_l) \mathcal{g}'(\mu\_l)} \right) + \sum\_{l=1}^{N} \frac{(\sum\_{h=1}^{N} A\_{l\cdot} \mathcal{C}\_{\cdot h} \eta\_h)^2}{a\_l \phi V(\mu\_l) (\mathcal{g}'(\mu\_l))^2}.$$

The Fisher information matrix is

$$\begin{aligned} \mathcal{Z}(\rho) &= \quad E(I(\rho)) = \sum\_{l=1}^{N} \frac{(\sum\_{h=1}^{N} A\_l \mathbb{C}\_{-h} \eta\_h)^2}{a\_l \Phi V(\mu\_l) (\mathcal{g}'(\mu\_l))^2} \\ &= \quad \frac{1}{\Phi} \sum\_{l=1}^{N} T\_l \left( \sum\_{h=1}^{N} A\_l \mathbb{C}\_{-h} \eta\_h \right)^2 = \frac{1}{\Phi} (AC\eta)' T(AC\eta) .\end{aligned}$$

*Appendix A.4. Algorithms of Coefficient Estimation*

We provide four algorithms to estimate the unknown coefficients of PR, GWPR, GWPAR and GTWPAR models in Section 4.

**Algorithm I:** Estimate the unknown parameters in the PR model. Take the initial values *g*(*μ*(0) *ik* ) = *η*(0) *ik* = ln *μ*(0) *ik* , *yik* = *μ*(0) *ik* , *Z*(0) *ik* = *η*(0) *ik* + *g* (*μ*(0) *ik* )(*yik* − *μ*(0) *ik* ), and *w*(0) *ik* = 1 *aikV*(*μ*(0) *ik* )(*g*(*μ*(0) , *i* = 1, 2, . . . , 10, *k* = 1, 2, . . . , 6.

The iterative formula of *β*ˆ(*m*+1) is

$$\beta^{(m+1)} = \left(X^\prime \mathcal{W}^{(m)} X\right)^{-1} X^\prime \mathcal{W}^{(m)} Z^{(m)}.$$

Repeat the above step until convergence yields. The estimated value *β*ˆ = *β*ˆ(*m*) can be obtained.

*ik* ))2

**Algorithm II:** Estimate the unknown coefficients in the GWPR model. (Note that *k* is a fixed constant taken from {1, 2, ... , 6} and the following steps should be repeated six times independently). Take the initial values *η*(0) *ik* = *g*(*μ*(0) *ik* ), *yik* = *μ*(0) *ik* , *Z*(0) *ik* = *η*(0) *ik* + *g* (*μ*(0) *ik* )(*yik* − *μ*(0) *ik* ), and

$$t\_{ik}^{(0)} = \frac{1}{a\_{ik} V(\mu\_{ik}^{(0)}) (g'(\mu\_{ik}^{(0)}))^2}, i = 1, 2, \dots, 10.$$

Let *Z*(0) = (*Z*(0) 1*k* , *Z*(0) 2*k* , ··· , *Z*(0) 10*k*) and *T*(0) = diag(*T*(0) 1*k* , *T*(0) 2*k* , ··· , *T*(0) 10*k* ). The iterative formula of *β*ˆ(*m*+1) at the location (*u*0*k*, *v*0*k*) is

$$\beta^{(m+1)}(\mathfrak{u}\_{0k}, \mathfrak{v}\_{0k}) \quad = \, \left( A'(\mathfrak{u}\_{0k}, \mathfrak{v}\_{0k}) T^{(m)} \mathcal{W}(\mathfrak{u}\_{0k}, \mathfrak{v}\_{0k}) A(\mathfrak{u}\_{0k}, \mathfrak{v}\_{0k}) \right)^{-1} A'(\mathfrak{u}\_{0k}, \mathfrak{v}\_{0k}) T^{(m)} \mathcal{W}(\mathfrak{u}\_{0k}, \mathfrak{v}\_{0k}) Z^{(m)},$$

where *W*(*u*0*k*, *v*0*k*) = diag(*w*1(*u*0*k*, *v*0*k*), *w*2(*u*0*k*, *v*0*k*), ··· , *w*10(*u*0*k*, *v*0*k*)) and

$$w\_i(u\_{0k}, v\_{0k}) = \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{1}{2} \left(\frac{d\_{ik}^{(0)}}{h}\right)^2\right).$$

Repeat the above step until convergence. When (*u*0*k*, *v*0*k*) takes all the locations (*uik*, *vik*), we will get the estimated value *β*ˆ = *β*ˆ(*m*) in a fixed the *k*th year.

**Algorithm III:** Estimate unknown coefficients in the GWPAR model. (Note that *k* is a fixed constant and the following steps should be repeated six times as in Algorithm II). Take the initial value *β*(0) 1 , *β*(0) 2 from Algorithm II, and *ρ*(0) is the absolute value of spatial Moran's I = 0.4480. The initial values *η*(0) = (*I* − *ρ*(0)*C*)−1*Xβ*(0), *μ*(0) *ik* = *g*−1(*η*(0) *ik* ), and

$$Z\_{ik}^{(0)} = \lg(\mu\_{ik}^{(0)}) + \lg'(\mu\_{ik}^{(0)})(y\_{ik} - \mu\_{ik}^{(0)}),\\t\_{ik}^{(0)} = \frac{1}{a\_{ik}V(\mu\_{ik}^{(0)})(\mathcal{g'}(\mu\_{ik}^{(0)}))^2}, i = 1, 2, \cdots, 10.$$

The (*m* + 1)th iterative estimation *β*ˆ(*m*+1)(*u*0*k*, *v*0*k*) and *ρ*ˆ(*m*+1) is

$$\begin{array}{rcl}\beta^{(m+1)}(\boldsymbol{u}\_{0k},\boldsymbol{v}\_{0k})&=& ( (\boldsymbol{A}^{(m)}\boldsymbol{X})^{\prime}\boldsymbol{T}^{(m)}\boldsymbol{\mathcal{W}}(\boldsymbol{A}^{(m)}\boldsymbol{X}) )^{-1}(\boldsymbol{A}^{(m)}\boldsymbol{X})^{\prime}\boldsymbol{T}^{(m)}\boldsymbol{\mathcal{W}}\boldsymbol{Z}^{(m)},\\\beta^{(m+1)}&=& \boldsymbol{\rho}^{(m)} + ((\boldsymbol{A}^{(m)}\boldsymbol{C}\boldsymbol{\eta}^{(m)})^{\prime}\boldsymbol{T}^{(m)}(\boldsymbol{A}^{(m)}\boldsymbol{C}\boldsymbol{\eta}^{(m)}))^{-1}(\boldsymbol{A}^{(m)}\boldsymbol{C}\boldsymbol{\eta}^{(m)})^{\prime}\boldsymbol{T}^{(m)}(\boldsymbol{Z}^{(m)} - \boldsymbol{\eta}^{(m)}).\end{array}$$

If (*u*0*k*, *v*0*k*) takes all the locations (*uik*, *vik*), the estimate *β*ˆ(*m*+1) can be given. When all estimated values arrive to converge, we will get *β*ˆ = *β*ˆ(*m*)(*uik*, *vik*) and *ρ*ˆ = *ρ*ˆ(*m*) in a fixed the *k*th year.

**Algorithm IV:** Estimate the unknown coefficients in the GTWPAR model. Take the initial values *β*(0) 1*k* (*uik*, *vik*, *tk*), *β*(0) 2*k* (*uik*, *vik*, *tk*), *i* = 1, 2, ... , 10, *k* = 1, 2, ... , 6 from Algorithm III, and *ρ*(0) is the absolute value of spatiotemporal Moran's I= 0.2143. The initial value vector *η*(0) = (*I* − *ρ*(0)*C*)−1*Xβ*(0), *μ*(0) *ik* = *g*−1(*η*(0) *ik* ), *Z*(0) = (*Z*(0) 1 , *Z*(0) 2 , ··· , *Z*(0) 60 ) , *T*(0) = diag(*T*(0) 1 , *T*(0) 2 ,..., *T*(0) 60 ) and

$$Z\_{ik}^{(0)} = \lg(\mu\_{ik}^{(0)}) + \lg'(\mu\_{ik}^{(0)})(y\_{ik} - \mu\_{ik}^{(0)}),\ T\_{ik}^{(0)} = \frac{1}{a\_{ik}V(\mu\_{ik}^{(0)})(\mathcal{g'}(\mu\_{ik}^{(0)}))^2}.$$

The (*m* + 1)th iterative estimations *β*ˆ(*u*00, *v*00, *t*0) and *ρ*ˆ are

$$\begin{array}{rcl}\hat{\beta}^{(\mathfrak{m}+1)}(\mathfrak{u}\_{00},\mathfrak{v}\_{00},t\_{0})&=&((A^{(\mathfrak{m})}X)'T^{(0)}W(A^{(\mathfrak{m})}X))^{-1}\\&\times(A^{(\mathfrak{m})}X)'T^{(\mathfrak{m})}WZ^{(\mathfrak{m})},\\\rho^{(\mathfrak{m}+1)}&=&\rho^{(\mathfrak{m})}+((A^{(\mathfrak{m})}\mathbb{C}\eta^{(\mathfrak{m})})'T^{(\mathfrak{m})}(A^{(\mathfrak{m})}\mathbb{C}\eta^{(\mathfrak{m})}))^{-1}\\&\times(A^{(\mathfrak{m})}\mathbb{C}\eta^{(\mathfrak{m})})'T^{(\mathfrak{m})}(Z^{(\mathfrak{m})}-\eta^{(\mathfrak{m})}).\end{array}$$

where *W* = {*wik*(*u*00, *v*00, *t*0)} and

$$w\_{lk}(\mu\_{00}, v\_{00}, t\_0) = \frac{1}{\sqrt{2\pi}} \exp\left\{ -\frac{1}{2} \left( \frac{(\mu\_{00} - \mu\_{lk})^2 + (v\_{00} - v\_{lk})^2}{h\_S^2} + \frac{(t\_0 - t\_k)^2}{\tau h\_S^2} \right) \right\}.$$

A detailed definition is given in Section 3.1. If (*u*00, *v*00, *t*0) takes all the locations (*uik*, *vik*, *tk*) and all estimations converge, we will get *β*ˆ = *β*ˆ(*m*)(*uik*, *vik*, *tk*) and *ρ*ˆ = *ρ*ˆ(*m*).

It is worth noting that we use the parameter estimates of the previous model as the initial values of the next model to reduce the number of iterations and improve the operational efficiency. For example, the estimations of the GWPR model are selected as the initial values of the GWPAR model.

*Appendix A.5. The Model Comparison Indicators (See Table 3)*

(1) The coefficient of determination is defined by

$$\mathbb{R}^2 = 1 - \frac{\sum (\eta - \hat{\eta})^2}{\sum (\eta - \bar{\eta})^2},$$

where *η* is a set of vectors {*ηik*}, *η*ˆ and *η*¯ are the parameter estimate and the mean value of *η*, respectively. (2) Deviation can be defined as

$$\mathcal{D} = \sum (y \ln(\frac{\hat{\mu}}{y}) + (y - \hat{\mu})),$$

where *y* is a set of response variables, and *μ*ˆ is the estimation of *μ*=E(y).

(3) The corrected Akaike information criterion is

$$\text{AICc} = D + 2P + 2\frac{P(P+1)}{N-P-1}r$$

where *D*, *P* and *N* are the deviation, the number of parameters and the number of samples, respectively. (4) Mean-square error is given by

$$\text{MSE} = \frac{1}{N - P} \sum (\vec{\eta} - \eta)^2 \vec{\eta}$$

where the parameter settings are the same as above.

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).


**1. Background**

*1.1. The New Urbanization "Beyond Mere Growth" in China*

Since 2014, the Chinese central government has promoted new-type urbanization in the *National New-type Urbanization Plan (2014–2020)* [1–3]. The direction of the plan "beyond mere growth" supports green and low-carbon development and requires avoiding urban sprawl, low efficiency in land use, the proliferation of useless urban districts, and waste of urban space in over-dimensioned infrastructures [4,5]. Consistent with the goals of new-type urbanization, the national guidelines on urban planning released in 2016 prohibit expanding cities beyond what their natural resources can support, enforcing urban growth boundaries [6]. In other words, to be sustainable urban development should not expand beyond a planned threshold, as also required by the environmental protection laws issued since 2015. This limit to expansion is also known as the "red line of China's cultivated land", and was first drafted by the National Land Planning Outline issued by the Ministry of Land and Natural Resources in 2009, which also specifies that the impermeabilization of the ground should be diminished [7].

These policies for preserving agricultural land are combined with the regulations for controlling land supply and the real estate market. Since 2003, the central government established year-by-year construction land development quotas for every local government according to the assessment of the local demand–supply situation [8–10].

**Citation:** Chen, J.; Pellegrini, P.; Yang, Z.; Wang, H. Strategies for Sustainable Urban Renewal: Community-Scale GIS-Based Analysis for Densification Decision Making. *Sustainability* **2023**, *15*, 7901. https://doi.org/10.3390/su15107901

Academic Editor: Miguel Amado

Received: 6 April 2023 Revised: 8 May 2023 Accepted: 9 May 2023 Published: 11 May 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Together with the sustainability goals, the increase in the urbanization rate is a key national priority in China. It has been one of the most significant drivers of economic, political, and social progress since the reform was promoted at the end of the 1970s [11–14]. So far, urban population growth has followed economic growth because economic development and city growth are strongly related to Chinese policies. The massive urbanization process has transformed millions of farmers into city dwellers, attracting an essential workforce into the expanding city [15,16].

These recent policies promoted by the Chinese government—on one side the goal of a growing economy and on the other side the limits to urbanization to protect fertile land and natural environments—can produce a contradiction and urban concentration: probable newcomers in growing cities and industrial districts versus a not expanding city.

#### *1.2. The Solution of High Density for Growing Cities and Sustainability*

Multiple actions should be considered for promoting urban sustainability, such as limits to new urbanizations and conservation of agricultural land, investments in public transportation [17], limits to car-oriented projects [18], reduction in energy consumption and greenhouse gases emission in all sectors of production, water management [19], and improving water, soil, and air quality [20]. Of the multiple elements that impact the sustainability of urban development, this paper focuses on the quantitative parameter of population density on urbanized land, more specifically on the opportunity to increase the dwelling units in existing residential neighborhoods. The research presented in this paper applies urban-scale GIS analysis to make some hypotheses on densification at the community scale.

Density is neither a good nor a bad indicator: too much can result in overcrowding and eventual urban decay; too little can deprive the dwelling environment of the chance of social interaction and public transport efficiency [21]. The definition of density is, nevertheless, closely connected to sustainability because land is a limited resource; urbanization is "sustainable" when valuable agricultural land is saved, and the ecosystem is not jeopardized.

High density has both positive and negative effects. Many scholars have presented evidence that higher density in urban areas is associated with a variety of desirable outcomes, including increased use of public transportation, improved financial stability for local governments, walkable and healthy living environments, housing diversity and affordability, enhanced community character, and cultural vitality [9–13,22–28]. Other scholars, on the contrary, have noted how some high-density neighborhoods produce excessive concentration and congestion, difficult management, and inefficient facilities. The issue cannot be assessed in general terms because the output of high-density conditions depends on several context-related factors, such as urban planning strategies, technical capacity, building regulations, infrastructure, public utilities, regional services, and the economic conditions of the residents [29]. For example, the interplay between density and building types must be considered because the same FAR realized by different architectural forms and open spaces produces different urban environments. In China, the different ranges of FAR—usually classified as low, medium, and high in building regulations and master plans—often correspond to specific building types; for example, medium density is realized with multilayer condos [30–33].

Over the last two decades, the idea of urban densification has been explored worldwide in several large cities and also in medium or high-density cases [34–38]. These densification projects are promoted mainly for the following reasons:

Both reasons justify the actions of densification promoted in dense cities, such as Paris [39,40], Seoul [41–43], Rotterdam [44–46], or London [47–49]; policies and projects differ among cities because they adapt to local conditions, including the socio-cultural approach to dwelling.

Regarding China's concerns, often urban land expands faster than population growth. In recent years, this trend has impacted several Chinese cities, which have experienced declining densities, very likely due to the rising average income leading to the demand for more residential space per capita, and for the restrictions due to the hukou system (a household registration which prevents the change of residence in order to limit the number of permanent residents in cities) [50,51]. According to the World Bank, even though China underwent massive urban population growth, its urban population density in 2010 was lower than the average of the rest of the East Asia region [52].


**Chapter 2**

References

Renal Diseases

63:537-560

1990;XVI:585-589

4(6):1073-1082

647-654

4

[1] Luiano RL, Moekel GW. Update on

curriculum 2019. American Journal of Kidney Diseases. 2019;79:404-415

[2] Pollak VE, Pirani C, Schwarz FD. The

Laboratory and Clinical Medicine. 1964;

[4] Visconti L, Camaro V, Ricciardi CA, Laava V, Pelliano V, Laoquaniti A, et al. Renal biopsy - still a landmark for the nephrologist. World Journal of Nephrology. 2016;5(4):321-327

[5] Moutzouris DA, Herlitz L, Appel GB,

Markowitz GB, Freudenthal B, Radhkrshman R, et al. Renal biopsy in the very elderly. Clinical Journal of the American Society of Nephrology. 2009;

[6] Ramirez G, Saba SR. Primary glomerulonephritis in the elderly. In: Zawada ET Jr, Sica D, editors. Geriatric Nephrology and Urology. Littleton Massahussetts: PSG Publishing Company, Inc.; 1985. pp. 49-66

[7] Hae-Yoon G, Bomback AS,

Stokes MB, Santoriello D, Campenot ES, Batal I, et al. The spectrum of kidney biopsy findings in patients with morbid obesity. Kidney International. 2019;95:

[3] Kark RM. The development of percutaneous renal biopsy in man. American Journal of Kidney Diseases.

the native kidney biosy: Core

natural history of the renal manifestations of systemic lupus erythematosus. The Journal of


**7. Changes in gut barrier defense following PN**

The epithelial barrier of the intestine, which turns over rapidly due to highly proliferative pluripotent Lgr5+ stem cells in the intestinal crypts, serve as the first line of defense against the external environment of the gut. The epithelium turns over every 3–5 days [66]. One cell type that does not turn over rapidly are small intestinal Paneth cells, that turn over every 20–30 days. The selective barrier allows absorption of water, electrolytes, and some macromolecules via tight-junction proteins, including zonulins, occludins, and claudins between enterocytes. Epithelial permeability increases significantly during PN feeding [67, 68], in parallel with a loss of tight-junction proteins [69]. In addition to the physical barrier along the gut lining, subepithelial dendritic cells extend dendrites between epithelial cells, which are hypothesized to sample luminal antigens and augment barrier responses [70]. A smaller proportion (10%) of intestinal cells secretory cells, including enteroendocrine cells and mucous secreting goblet cells [47].

The epithelium's contribution to defense includes not only physical barrier formation, but release of antimicrobial molecules that influence microbial community composition and membership. Enterocytes comprise 90% of total epithelial cells in the gut and release β-defensins and RegIIIγ enzymes that limit microbial growth at the mucous barrier [71]. The far less abundant Paneth cells, found at the crypt bases, produce a large array of antimicrobial peptides and enzymes, including lysozyme, RegIIIγ, secretory PLA2, Angiogenin4, and α-defensins (murine cryptdins) [72]. These cationic antimicrobials localize to the negatively-charged mucous surface and work by targeting conserved aspects of microbial cell walls and membranes [67]. Studies demonstrate Paneth cell antimicrobial molecules reach 15–100 mg/mL within intestinal crypts, exceeding antimicrobial concentrations [73]. The colon exhibits two layers of mucous, an inner sterile layer and an outer more loosely colonized layer.

Paneth cell antimicrobial release is regulated by Th2 cytokines, including IL-4, IL-9, and IL-13, GLP-2, and insulin [37, 74, 75]. Stimulatory triggers include the ligands TLRs and NOD2 and parasympathetic cholinergic stimulation [76]. PN decreases antimicrobial production through lower levels of intestinal IL-4 and IL-13 [77, 78]. Enterocyte release of RegIIIγ is also lost during PN [79]. Exogenous administration of IL-25 stimulates production of IL-4 and IL-13 cytokines and levels of Lysozyme, sPLA2, and RegIIIγ compared with PN feeding alone [80]. Following PN, mucosal secretions contain less Paneth cell antimicrobial products, leading to decreased killing of bacteria in vitro [76]. Tissue explants from PN fed animals have been demonstrated to be more susceptible to enteroinvasive *E. coli* compared with controls [81].

The function of goblet cells is to produce the mucous barrier, composed of glycoprotein mucins. These proteins have numerous carbohydrate residues including *O*-glycosylation, *N*-actyl-galactosamine, galactose, and *N*-actyl-glucosamines [82]. Functionally, mucins mucous serves as a selective physical barrier that allows for movement of luminal digesta, absorption of nutrients, and limiting bacterial access to the gut wall. Mucins also help concentrate Paneth cell antimicrobial molecules and sIgA through charge interactions [83]. The importance of mucins is demonstrated by in mutant animals lacking the more abundant mucin, MUC2, which leads to inflammatory enteritis and increased risk of tumor formation [45].

PN decreases the release of MUC2, RELMβ, and trefoil factor 3 (TFF3) in neonatal piglets and adult mice [75, 84]. TFF3 promotes epithelial response to injury. Animals deficient in RELMβ display increased susceptibility to *Citrobacter rodentium* challenge. The release of goblet cell products are influenced by Th2 cytokines, including IL-4 and IL-13, which are decreased under PN feeding [85]. Exogenous administration of the Th2 stimulating cytokine, IL-25, elevates luminal MUC2 levels [77].

Although only 1% of the total intestinal epithelial cells are enteroendocrine cells (EECs), they collectively make up the largest endocrine organ, expressing almost 2 dozen peptide hormones that shape metabolism, immunity, and behavior [86, 87]. EECs respond to enteral nutrients as well as microbial ligands [88]. EECs are also intricately linked to innate immunity, indirectly activating and recruiting immune cells by producing the chemokines CXCL-1, CXCL-3, and the cytokine IL-32. EEC hormones can also influence epithelial cell function in the gut, including Paneth cell release of antimicrobial molecules. GLP-2 mutant animals are at increased susceptibility to gut infection than wild-type littermates [37].


*5.2. Research Models*

This paper examines the association between CG and the performance of Chinese SOEs by using regression analysis on panel data. The following two models are used, where *β*0 is the intercept, *β* is the regression coefficient, and *εi* is an error term.