Estimation Curve of Mixed Spline Truncated and Fourier Series Estimator for Geographically Weighted Nonparametric Regression

Abstract

Geographically Weighted Regression (GWR) is the development of multiple linear regression models used in spatial data. The assumption of spatial heterogeneity results in each location having different characteristics and allows the relationships between the response variable and each predictor variable to be unknown, hence nonparametric regression becomes one of the alternatives that can be used. In addition, regression functions are not always the same between predictor variables. This study aims to use the Geographically Weighted Nonparametric Regression (GWNR) model with a mixed estimator of truncated spline and Fourier series. Both estimators are expected to overcome unknown data patterns in spatial data. The mixed GWNR model estimator is then determined using the Weighted Maximum Likelihood Estimator (WMLE) technique. The estimator’s characteristics are then determined. The results of the study found that the estimator of the mixed GWNR model is an estimator that is not biased and linear to the response variable y.

1. Introduction

Regression analysis is a statistical method used to determine the relationship between response variables and one or more predictor variables [1]. Regression is divided into three types: parametric regression, nonparametric regression, and semiparametric regression. Parametric regression is used when the shape of the regression curve is known, whether linear, quadratic, cubic or otherwise. Whereas, nonparametric regression is used when the shape of the regression curve is unknown. As for semiparametric regression, it is a combination of parametric and nonparametric regression.

Nonparametric regression is a method used to model regression curves of unknown shape [2]. This method is a more flexible approach because the data is expected to look for the estimation form of the regression curve itself without being influenced by the researcher’s subjectivity factor [3]. Some of the estimators used in nonparametric regression include spline estimators, Fourier series, kernels, and local polynomials. Each estimator has its characteristics to approach unknown regression functions. Research on nonparametric regression has been widely conducted by single estimators [4,5] and mixed estimators [6,7]. However, its application is still limited to nonspatial data. In fact, there are many problems related to spatial data. Spatial data is data that contains size and location information [8]. Methods used in spatial data analysis include Spatial Autocorrelation, Spatial Error Model, Geographically weighted regression, and others.

According to [8], Geographically Weighted Regression (GWR) is a statistical method that can analyze spatial heterogeneity. Spatial heterogeneity is one of the same predictor variables exerting unequal influences on different locations within a study site. The GWR model generates an estimator of model parameters that are local to each point or location where the data is observed. Research on GWR by [9,10] shows that the GWR model is better than the global model that can overcome spatial heterogeneity. It was further developed on a nonparametric GWR with a single estimator [11,12,13]. In addition, the relationship between response variables and some predictor variables can vary [14]. Therefore, a nonparametric GWR model with mixed truncated spline and Fourier series would be developed. The truncated spline estimator in GWR is expected to overcome the changing curve pattern at certain sub intervals [11]. In contrast, the Fourier series estimator is expected to model the repeating data pattern [6]. This study aims to create a Geographically Weighted Nonparametric Regression (GWNR) model with a mixed truncated spline and Fourier series estimator, to determine the parameter estimate of the GWNR model with the Weighted Maximum Likelihood Estimator method, and evaluate the properties the mixed GWNR model.

The following discussion in this paper is divided into three main topics. Section 2 discusses the GWNR Model and method estimation of WMLE. Section 3 presents the estimation parameter model GWNR; unbiased and linear estimator properties; and data application. Section 4 is the conclusions.

2. Materials and Methods

2.1. Geographically Weighted Nonparametric Regression (GWNR)

The GWNR model is a development of GWR in nonparametric regression. Provided paired data $\left(x_{1i},\dots ,x_{Pi},z_{1i},\dots ,z_{Qi},y_i\right)$ and assumed relationships between predictor variables $\left(x_{1i},\dots ,x_{Pi},z_{1i},\dots ,z_{Qi}\right)$ with response variables $\left(y_i\right)$ following a multivariable regression model [11] are as follows:

(1)$y_i=\mu \left(x_{1i},\dots ,x_{Pi},z_{1i},\dots ,z_{Qi}\right)+{\epsilon }_i,i=1,2,\dots ,n$

where, $y_i$ is a response variable, $\mu \left(x_{1i},\dots ,x_{Pi},z_{1i},\dots ,z_{Qi}\right)$ is a regression curve of unknown shape, P is the predictor variables with a truncated spline function, Q is the predictor variables approached with Fourier series functions, and n is a number of observations and is assumed to be additive. If the function $\mu \left(x_{1i},\dots ,x_{Pi},z_{1i},\dots ,z_{Qi}\right)$ is approached with truncated spline functions and Fourier series, then Equation (1) can be written:

$y_i=\mu \left(x_{1i},\dots ,x_{Pi},z_{1i},\dots ,z_{Qi}\right)+{\epsilon }_i$

$y_i={\displaystyle \sum _{p=1}^P{f}_p\left(x_{pi}\right)}+{\displaystyle \sum _{q=1}^Q{g}_q\left(z_{qi}\right)}+{\epsilon }_i,i=1,2,\dots ,n$

where,

$\begin{array}{ll}{\displaystyle \sum _{p=1}^P{f}_p}& {\displaystyle \left(x_{pi}\right)}=f_1\left(x_{1i}\right)+\dots +f_P\left(x_{Pi}\right)\\ & =X_1{\mathsf{\beta }}_1(u_i,v_i)+\dots +X_P{\mathsf{\beta }}_P(u_i,v_i)\\ & =\left[\begin{array}{ccc}{X}_1& \dots & X_P\end{array}\right]\left[\begin{array}{c}{\mathsf{\beta }}_1(u_i,v_i)\\ ⋮\\ {\mathsf{\beta }}_P(u_i,v_i)\end{array}\right]\\ & =X\mathsf{\beta }(u_i,v_i)\end{array}$

is a truncated spline component with P predictor variables and

$\begin{array}{ll}{\displaystyle \sum _{q=1}^Q{g}_q}& {\displaystyle \left(z_{qi}\right)}=g_1\left(z_{1i}\right)+\dots +g_Q\left(z_{Qi}\right)\\ & =Z_1{a}_1(u_i,v_i)+\dots +Z_Q{a}_Q(u_i,v_i)\\ & =\left[\begin{array}{ccc}{Z}_1& \cdots & Z_Q\end{array}\right]\left[\begin{array}{c}{a}_1(u_i,v_i)\\ ⋮\\ a_Q(u_i,v_i)\end{array}\right]\\ & =Za(u_i,v_i)\end{array}$

is a Fourier series component with Q other predictor variables. By matrix notation, it can be written with:

(2)$y=X\mathsf{\beta }(u_i,v_i)+Za(u_i,v_i)+\mathsf{\epsilon }$

where:

$y=\left[\begin{array}{c}{y}_1\\ ⋮\\ y_n\end{array}\right],\mathsf{\epsilon }=\left[\begin{array}{c}{\epsilon }_1\\ ⋮\\ {\epsilon }_n\end{array}\right]$

$X=\left[\begin{array}{ccccc}{X}_1& ⋮& \cdots & ⋮& X_P\end{array}\right]$

$Z=\left[\begin{array}{ccccc}{Z}_1& ⋮& \cdots & ⋮& Z_Q\end{array}\right]$

$\mathsf{\beta }(u_i,v_i)={\left[\begin{array}{ccccc}{\mathsf{\beta }}_1{}^T(u_i,v_i)& ⋮& \cdots & ⋮& {\mathsf{\beta }}_P{}^T(u_i,v_i)\end{array}\right]}^T$

$a(u_i,v_i)={\left[\begin{array}{ccccc}{a}_1{}^T(u_i,v_i)& ⋮& \cdots & ⋮& a_Q{}^T(u_i,v_i)\end{array}\right]}^T$

$(u_i,v_i)=(\mathrm{longitude},\mathrm{latitude}),i=1,2,\dots ,n$

2.2. Weighted Maximum Likelihood Estimator (WMLE)

Maximum likelihood estimation from the parameters $\mu$ and ${\sigma }^2$ with a distributed n-sized sample $y_i\sim N(\mu ,{\sigma }^2)$, $i=1,2,\dots ,n$ can be written with:

$f\left(y_1,y_2,\dots ,y_n\right)={\displaystyle \prod _{i=1}^n\frac{1}{\sqrt{2\pi {\sigma }^2}}}\exp {\left(-\frac{1}{2}\left(\frac{y_i-\mu }{\sigma }\right)\right)}^2$

so that the likelihood function for $y_i,i=1,2,\dots ,n$ is

(3)$L\left(\mu ,{\sigma }^2|y_i\right)={\left(2\pi \right)}^{-\frac{n}{2}}{\left({\sigma }^2\right)}^{-\frac{n}{2}}\exp \left(-\frac{1}{2{\sigma }^2}{\displaystyle \sum _{i=1}^n{\left(y_i-\mu \right)}^2}\right)$

Next, by multiplying by the weighting matrix $W=diag\left(w_{1i^{*}},w_{2i^{*}},\dots ,w_{n{i}^{*}}\right)$, the log likelihood function is obtained in terms of weighted function $w_{ii*}$ as [15]:

(4)$\ln L^{*}\left(\mu ,{\sigma }^2\left(u_{i*},v_{i*}\right)|y_i\right)={\displaystyle \sum _{i=1}^n{w}_{ii*}}\ln \left(\frac{1}{\sqrt{2\pi {\sigma }^2\left(u_{i*},v_{i*}\right)}}\exp \left(-\frac{1}{2}{\left(\frac{y_i-\mu }{\sigma \left(u_{i*},v_{i*}\right)}\right)}^2\right)\right),i=1,2,\dots ,n$

The role of weights in the GWNR model is very important because the weighting values will represent the location of observational data from one another. One method that can be used is the Gaussian Kernel function [8]. Equation (4) estimates the GWNR parameter with the WMLE method.

Furthermore, steps are given in estimating the parameters of a mixed GWNR model with the WMLE Method as follows:

• Defining a mixed GWNR model

• Assuming distribution $\mathsf{\epsilon }$

• Determining the distribution of y

• Forming a likelihood function

• Forming a weighted likelihood function

• Specifying the first partial derivative of the likelihood function against the mixed GWNR model parameter

• Getting an estimate of mixed GWNR model parameters.

3. Results

3.1. Parameter Estimation

Estimation of parameters on GWNR models with mixed estimators uses Weighted Maximum Likelihood Estimator (WMLE). The WMLE method is obtained by knowing the distribution of the response variable $y_i,i=1,2,\dots ,n$ in advance. Then, it is determined by the weighting matrix for each location to $i$, $i=1,2,\dots ,n$. The weighting used is a Fixed Gaussian Kernel function. Next, it is given the form of a distribution of the GWNR model that is presented on Lemma 1.

Lemma 1. 

Given the GWNR model in Equation (2), with ${\epsilon }_i,i=1,2,\dots ,n$ is normally distributed with mean equal to zero and variance ${\sigma }^2\left(u_i,v_i\right)$, hence $y_i,i=1,2,\dots ,n$ is normally distributed with mean

$\begin{array}{}{\beta }_0\left(u_i,v_i\right)+{\displaystyle \sum _{p=1}^P{\displaystyle \sum _{m=1}^M{\beta }_{mp}}}\left(u_i,v_i\right)x_{pi}^m+{\displaystyle \sum _{p=1}^P{\displaystyle \sum _{r=1}^R{\beta }_{(r+M)p}}}\left(u_i,v_i\right){\left(x_{pi}-t_{rp}\right)}_{+}^m\\ +{\displaystyle \sum _{q=1}^Q\left({\gamma }_q\left(u_i,v_i\right)z_{qi}+\frac{1}{2}{\theta }_{0q}\left(u_i,v_i\right)+{\displaystyle \sum _{h=1}^H{\theta }_{hq}\left(u_i,v_i\right)\cos \left(h{z}_{qi}\right)}\right)}\end{array}$

and variance, ${\sigma }^2\left(u_i,v_i\right)$.

where:

• P: number of spline components

• M: polynomial degree of spline

• R: number of knot points

• Q: number of Fourier components

• H: number of oscillation parameters.

Lemma 1 has been proven in Appendix A.

Theorem 1. 

If given a model on Equation (2) with ${\epsilon }_i,i=1,2,\dots ,n$ normally distributed with mean zero and variance ${\sigma }^2\left(u_i,v_i\right)$ and the weighted likelihood function given to (4), by the MLE method, an estimator is obtained $\widehat{\mathsf{\beta }}\left(u_i,v_i\right)$ and $\widehat{a}\left(u_i,v_i\right)$ as follows:

$\widehat{\mathsf{\beta }}\left(u_i,v_i\right)=A(t,h)y$

$\widehat{a}\left(u_i,v_i\right)=B(t,h)y$

where:

$\begin{array}{}A(t,h)=R{\left(X^{T}W\left(u_i,v_i\right)X\right)}^{-1}[X^{T}W\left(u_i,v_i\right)-X^{T}W\left(u_i,v_i\right)Z\\ {\left(Z^{T}W\left(u_i,v_i\right)Z\right)}^{-1}{Z}^{T}W\left(u_i,v_i\right)]\\ {\rm B}(t,h)=S{\left(Z^{T}W\left(u_i,v_i\right)Z\right)}^{-1}[Z^{T}W\left(u_i,v_i\right)-Z^{T}W\left(u_i,v_i\right)X\\ {\left(X^{T}W\left(u_i,v_i\right)X\right)}^{-1}{X}^{T}W\left(u_i,v_i\right)]\end{array}$

• t = knot point for spline component

• h = oscillation parameter component.

Proof of Theorem 1. 

Is given to Appendix B. □

3.2. Unbiased and Linear Estimator Properties

Lemma 2. 

If $\widehat{\mathsf{\beta }}\left(u_i,v_i\right)$is a truncated spline component parameter estimator of the GWNR model with a mixed estimator approach that follows Equation (4), so $\widehat{\mathsf{\beta }}\left(u_i,v_i\right)$is an unbiased estimator and belongs to the class of linear estimators in observation y.

Furthermore, it can be seen in Appendix C which is the proof of Lemma 2.

Lemma 3. 

If $\widehat{a}\left(u_i,v_i\right)$is a Fourier series component parameter estimator of the GWNR model with a mixed estimator approach that follows Equation (4), so $\widehat{a}\left(u_i,v_i\right)$is an unbiased estimator and belongs to the class of linear estimators in observation y.

Lemma 3 is the last proven lemma and is described in Appendix D.

Lemma 4. 

If  $\widehat{\mathsf{\beta }}\left(u_i,v_i\right)$  and  $\widehat{a}\left(u_i,v_i\right)$  are given by Theorem 1, hence the estimator for  $\widehat{f}$ ,  $\widehat{g}$  and  $\widehat{\mathsf{\mu }}$  is hence given by:

$\widehat{f}=X\widehat{\mathsf{\beta }}\left(u_i,v_i\right)$

$\widehat{g}=Z\widehat{a}\left(u_i,v_i\right)$

so that the following is obtained:

$\widehat{\mathsf{\mu }}=X\widehat{\mathsf{\beta }}\left(u_i,v_i\right)+Z\widehat{a}\left(u_i,v_i\right)=C(t,h)y$

Proof.  

Next, to determine the function estimator $\widehat{f}$, $\widehat{g}$ and $\widehat{\mathsf{\mu }}\left(u_i,v_i\right)$ are described as follows. Based on Theorem 1, it can be substituted $\widehat{\mathsf{\beta }}\left(u_i,v_i\right)$ and $\widehat{a}\left(u_i,v_i\right)$ so that it is obtained:

$\begin{array}{}\widehat{f}=X\widehat{\mathsf{\beta }}\left(u_i,v_i\right)\\ =XA(t,h)y\end{array}$

and

$\begin{array}{}\widehat{g}=Z\widehat{a}\left(u_i,v_i\right)\\ =ZB(t,h)y\end{array}$

As a result, obtained estimator $\widehat{\mathsf{\mu }}$ be

$\begin{array}{}\widehat{\mathsf{\mu }}=\widehat{f}+\widehat{g}\\ =X\widehat{\mathsf{\beta }}\left(u_i,v_i\right)+Z\widehat{a}\left(u_i,v_i\right)\\ =XA(t,h)y+ZB(t,h)y\\ =C(t,h)y\end{array}$

where $C(t,h)=XA(t,h)+ZB(t,h)$ is a hat matrix containing knot points t and the parameter of the h oscillation on the mixed GWNR model with the approach of truncated spline and Fourier series estimator. □

3.3. Data Application

The data used are secondary data from [16,17,18,19,20,21] with research variables, percentage of the poor population $\left(y\right)$, CPI $\left(x_1\right)$, TPT $\left(x_2\right)$, longitude and latitude coordinates $\left(u_i,v_i\right)$, and as many as 81 districts/cities on Sulawesi Island.

The following steps for applying the estimated parameters of the GWNR model to poverty data on Sulawesi Island in 2020 are as follows:

• Making a scatter plot between the variables $x_1$ and y, as well as $x_2$ and y

• Defining the initial model

• Selecting optimum knots and oscillation parameters

• Estimating parameters of global model with the OLS method based on the initial model formed

• Testing assumptions of spatial heterogeneity on residual values on global models

• Determining the weighting matrix

• Estimating parameters of the GWNR model with the WMLE method

• Choosing the best model based on MSE and R²

• Making conclusions

Based on the above step, the variable $x_1$ as a component of the Fourier series and the variable $x_2$ as the spline component are obtained based on the scatter plot shown. Furthermore, the test of spatial assumptions is obtained that the assumption of spatial heterogeneity is met, so the global model is less suitable for use because the residual properties are not homogeneous. One alternative model that can be used is the GWNR model. Use of this GWNR model is expected to overcome heteroskedasticity by generating a local model for each location. Here are some local models generated:

(5)${\widehat{y}}_{ken}=47.59-0.5x_1-0.17x_2+0.43\cos (x_1)+0.22(x_2-8.55)$

(6)${\widehat{y}}_{mks}=49.9-0.58x_1-0.07x_2+0.08\cos (x_1)+0.47(x_2-8.55)$

(7)${\widehat{y}}_{man}=44.19-0.46x_1-0.15x_2+0.11\cos (x_1)+0.19(x_2-8.55)$

(8)${\widehat{y}}_{pal}=40.79-0.47x_1+0.3x_2-0.48\cos (x_1)-0.27(x_2-8.55)$

(9)${\widehat{y}}_{gor}=44.12-0.3x_1-2.25x_2-0.05\cos (x_1)+2.38(x_2-8.55)$

(10)${\widehat{y}}_{maj}=39.88-0.43x_1+0.04x_2+0.02\cos (x_1)+0.22(x_2-8.55)$

where:

• y_ken = estimated poverty percentage for Kendari City

• y_mks = estimated poverty percentage for Makassar City

• y_man = estimated poverty percentage for Manado City

• y_pal = estimated poverty percentage for Palu City

• y_gor = estimated poverty percentage for Gorontalo City

• y_maj = estimated poverty percentage for Mamuju City.

GWNR mixed with oscillation parameters k =1 and linear spline t =1 resulted in MSE and R² values of 3.65 and 74.65 per cent, respectively. Based on several local models above, it shows that poverty in Sulawesi Island is influenced by HDI and TPT, where the increasing HDI will result in a decrease in the percentage of poverty. Conversely, an increase in TPT will increase the percentage of poverty.

4. Conclusions

Estimation of GWNR using the truncated spline and Fourier series was successfully formulated. It was found that:

• The GWNR model using a mixed estimator of truncated spline and Fourier series is

  $y=X\mathsf{\beta }(u_i,v_i)+Za(u_i,v_i)+\mathsf{\epsilon }$

  Where $f=X\mathsf{\beta }\left(u_i,v_i\right)$ is a truncated spline component, $g=Za\left(u_i,v_i\right)$ is a component of a Fourier series, and $\mathsf{\epsilon }$ is a residual component.

• Estimators of GWNR are $\widehat{\mathsf{\beta }}\left(u_i,v_i\right)=A(t,h)y$, $\widehat{a}\left(u_i,v_i\right)=B(t,h)y$, and $\widehat{\mathsf{\mu }}=C(t,h)y$. The estimator is an unbiased and linear estimator to observe the response variable.

Author Contributions

Conceptualization: L.L., I.N.B. and V.R.; methodology, L.L. and I.N.B.; writing-original draft preparation, L.L.; writing—review and editing, L.L., I.N.B. and V.R. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

References

1. Draper, N.R.; Smith, H. Applied Regression Analysis; 3rd ed. John Wiley & Sons Inc.: Hoboken, NJ, USA, 2014; pp. 1-716. [DOI: https://dx.doi.org/10.1002/9781118625590]

2. Budiantara, I.N. The combination of spline and kernel estimator for nonparametric regression and its properties. Appl. Math. Sci.; 2015; 9, pp. 6083-6094. [DOI: https://dx.doi.org/10.12988/ams.2015.58517]

3. Cheng, M.; Paige, R.L.; Sun, S.; Yan, K. Variance reduction for kernel estimators in clustered/longitudinal data analysis. J. Stat. Plan. Inference; 2010; 140, pp. 1389-1397. [DOI: https://dx.doi.org/10.1016/j.jspi.2009.09.026]

4. Octavanny, M.A.D.; Budiantara, I.N.; Ratnasari, V. Pemodelan faktor-faktor yang memengaruhi provinsi jawa timur menggunakan pendekatan regresi semiparametrik spline. J. Sains Seni ITS; 2017; 6, pp. 1-7.

5. Hidayat, M.F.; Achmad, R.F.A.; Solimun,. Estimation of truncated spline function in non-parametric path analysis based on weighted least square (WLS). IOP Conf. Ser. Mater. Sci. Eng.; 2019; 546, pp. 5-10. [DOI: https://dx.doi.org/10.1088/1757-899X/546/5/052027]

6. Sudiarsa, W.; Budiantara, I.N.; Suhartono, S.; Purnami, S.W. Combined estimator fourier series and spline truncated in multivariable nonparametric regression. Appl. Math. Sci.; 2015; 9, pp. 4997-5010. [DOI: https://dx.doi.org/10.12988/ams.2015.55394]

7. Mariati, N.P.A.M.; Budiantara, I.N.; Ratnasari, V. Combination estimation of smoothing spline and fourier series in nonparametric regression. J. Math.; 2020; 2020, 4712531. [DOI: https://dx.doi.org/10.1155/2020/4712531]

8. Fotheringham, A.S.; Brundson, C.; Charlton, M. Geographically Weighted Regression, the Analysis of Spatially Varying Relationships; John Wiley & Sons Ltd.: Chichester, UK, 2002.

9. Tang, J.; Gao, F.; Liu, F.; Zhang, W.; Qi, Y. Understanding spatio-temporal characteristics of urban travel demand based on the combination of GWR and GLM. Sustainabily; 2019; 11, 5525. [DOI: https://dx.doi.org/10.3390/su11195525]

10. Dziauddin, M.F. Estimating land value uplift around light rail transit stations in Greater Kuala Lumpur: An empirical study based on geographically weighted regression (GWR). Res. Transp. Econ.; 2019; 74, pp. 10-20. [DOI: https://dx.doi.org/10.1016/j.retrec.2019.01.003]

11. Sifriyani,; Kartiko, S.H.; Budiantara, I.N.; Gunardi,. Development of nonparametric geographically weighted regression using truncated spline approach. Songklanakarin J. Sci. Technol.; 2018; 40, pp. 909-920.

12. Fitri, N.; Sifriyani, S.; Yuniarti, D. Nonparametric geographically weighted regression dengan pendekatan spline truncated. Pros. Semin. Nas. Mat. Dan Stat.; 2019; pp. 98-105.

13. Sifriyani, S. Simultaneous hypothesis testing of multivariable nonparametric spline regression in the GWR model. Int. J. Stat. Probab.; 2019; 8, 32. [DOI: https://dx.doi.org/10.5539/ijsp.v8n4p32]

14. Nurcahayani, H.; Budiantara, I.N.; Zain, I. The curve estimation of combined truncated spline and fourier series estimator for multiresponse nonparametric regression. Mathematics; 2021; 9, 1141. [DOI: https://dx.doi.org/10.3390/math9101141]

15. Akbarov, A.; Wu, S. Waranty claim forecasting based on weighted maximum likelihood estimator. Qual. Reliab. Eng. Int.; 2012; 28, pp. 663-669. [DOI: https://dx.doi.org/10.1002/qre.1399]

16. BPS. Provinsi Sulawesi Tenggara dalam Angka 2021; UD. Resky Bersama: Kendari, Indonesia, 2021.

17. BPS. Provinsi Sulawesi Selatan dalam Angka 2021; BPS Press: Makassar, Indonesia, 2021.

18. BPS. Provinsi Sulawesi Utara dalam Angka 2021; Perum Percetakan NRI: Manado, Indonesia, 2021.

19. BPS. Provinsi Sulawesi Tengah dalam Angka 2021; UD. Rio: Palu, Indonesia, 2021.

20. BPS. Provinsi Sulawesi Barat dalam Angka 2021; Erlangga: Mamuju, Indonesia, 2021.

21. BPS. Provinsi Gorontalo dalam Angka 2021; CV. Rifaldi: Gorontalo, Indonesia, 2021.