Poverty Modeling in North Sumatera Province Considering County Location Using Geographical Weighted Regression and LASSO

2.1 Multiple linear regression

Multiple linear regression is a model that explains the effect of one dependent variable (Y) with two or more independent variables $\left(X_1, X_2, \ldots, X_p\right)$. The multiple linear regression equation model can be expressed mathematically as follows [5]:

$Y_i=\beta_0+\beta_1 X_{i 1}+\cdots+\beta_p X_{i p}+\varepsilon_i ; i=1,2, \ldots, n$          (1)

where,

$Y_i$: the dependent variable value at the $i^{t h}$ observation,

$\beta_0$: intercept,

$\beta_k$: the $k^{t h}$ regression parameter,

$X_{i k}$: the $i^{t h}$ observable of the $k^{t h}$ independent variable, (k=1,2,3..., p)

$\varepsilon_i$: error at the $i^{t h}$ observation.

2.2 Spatial data

Spatial data refers to data related to location or space. It includes geographic information, such as geographic coordinates, borders, maps, satellite images, topographic data, and other information related to geographic dimensions.

2.3 Spatial heterogeneity test

Spatial heterogeneity is a condition where the global regression model is unable to explain between variables because of the diversity of characteristics between regions [6]. Spatial heterogeneity testing is done through Breusch-Pagan testing [6].

$B P=\left(\frac{1}{2}\right) f^T Z\left(Z^T Z\right)^{-1} Z^T f \sim \chi^2$          (2)

with:

$f_i=\frac{\varepsilon_i^2}{\sigma^2}-1$

where,

$\varepsilon_i$: error for the $i^{t h}$ observation,

$\sigma^2$: variance,

$Z$: matrix of size (n×1) which is the standardized value of X for each $i^{t h}$ observation.

2.4 Geographically Weighted Regression

Geographically Weighted Regression (GWR) is a development of global regression methods that takes into account location aspects and fulfills the assumption of spatial heterogeneity, resulting in different parameter estimates for each location. The GWR model depends on the weights used at each observation location based on their geographic location. The elements of the weight matrix are determined based on the regression point and the observation point. Large weights are obtained based on points that are close to the observation point. This weighting depends on the selection of the optimal bandwidth [7]. The GWR model is an extension of a regression model in which each parameter is estimated at each point [8]. The GWR model can be expressed as follows [9]:

$y_i=\beta_0\left(u_i, v_i\right)+\sum_{k=1}^p \beta_k\left(u_i, v_i\right) X_{i k}+\varepsilon_i$          (3)

where,

$y_i$: the value of the dependent variable at the $i^{\text {th }}$ observation location,

$X_{i k}$: the value of the $k^{t h}$ independent variable at the $i^{t h}$ observation location,

$u_i, v_i$: coordinates of point $i$ with $u_i$ is Longitude and $v_i$ is Latitude,

$\beta_k\left(u_i, v_i\right)$: $k^{t h}$ local parameter at observation location $\left(u_i, v_i\right)$,

In estimating GWR parameters, there are several steps, namely:

a. Calculating spatial weights

A spatial weighting matrix is essentially a matrix that shows how different regions relate to each other [10]. Spatial weights are contained in a diagonal matrix that shows the proximity between observation locations whose function is to estimate different parameters at each observation location [11].

$W\left(u_i, v_i\right)=\exp \left[-\frac{1}{2}\left(\frac{d_{i j}}{h}\right)^2\right]$          (4)

with [11]:

$d_{i j}=\sqrt{\left(u_i-u_j\right)^2+\left(v_i-v_j\right)^2}$          (5)

where, $d_{i j}$ is the euclidean distance measures the straight-line distance between two locations $\left(u_i, v_i\right)$ and $\left(u_j, v_j\right)$, while the smoothing parameter (bandwidth) controls the degree of smoothing applied in spatial analysis methods like kernel density estimation.

b. Optimum bandwidth selection

In Geographically Weighted Regression (GWR), bandwidth acts as a weight to achieve a balance between how well the curve fits the data points and the overall smoothness of the curve in the model. The best bandwidth value is the one that minimizes the Cross-Validation (CV) value, which is a measure of the model's predictive accuracy.

$C V(h)=\sum_{i=1}^n\left(y_i-\hat{y}_{\neq i}(h)\right)^2$          (6)

c. GWR parameter estimation

Parameter estimation in the GWR model uses the Weighted Least Square (WLS) method, which is by giving different weights to each observation location [12]. Parameter estimates based on the weighting matrix will minimize the sum of squared residuals or Sum Square Error (SSE) so that it is obtained:

$\begin{aligned} \sum_{i=1}^n W\left(u_i, v_i\right) \varepsilon_i^2 & =\sum_{i=1}^n W\left(u_i, v_i\right)\left[y_i-\beta_0\left(u_i, v_i\right)\right] -\sum_{k=1}^p \beta_k\left(u_i, v_i\right) X_{i k}\end{aligned}$          (7)

so that the parameter estimates of the GWR model for each observation location are obtained as follows:

$\hat{\beta}\left(u_i, v_i\right)=\left[X^T W\left(u_i, v_i\right) X\right]^{-1} X^T W\left(u_i, v_i\right)$          (8)

d. Hypothesis testing of GWR model

   1) Goodness of Fit (GoF) test on GWR aims to determine the best model between GWR model and multiple linear regression model.

   2) t-Test is performed to identify which independent variables have a significant individual impact on each dependent variable.

2.5 Multicollinearity

Multicollinearity was initially identified by Ragnar Frisch, who observed a linear relationship among some or all of the independent variables in a regression model [13]. When multicollinearity symptoms appear in the model, it leads to an increase in the variance of the regression coefficients. This rise in variance results in the following effects [14]:

1. Testing regression parameters using the t test becomes invalid.

2. There is a contradiction between the results of simultaneous parameter hypothesis testing through the F test and the results of partial regression parameter testing through the t test.

According to the study conducted by Darnius and Tambunan [14], multicollinearity can be detected using the Variance Inflation Factor (VIF) value. The VIF value can be found using the formula:

$V I F_k=\frac{1}{\left(1-R_k^2\right)}, k=1,2, \ldots, p$          (9)

$R_k^2$ is the coefficient of determination obtained from the independent variable $X_k$ regressed with other independent variables. If the $VIF_k$ value is greater than 10, there is a multicollinearity problem.

2.6 Least Absolute Shrinkage and Selection Operator

Definition

Consider data $\left(x^i, y_i\right), \boldsymbol{i}=\mathbf{1}, \mathbf{2}, \ldots, \boldsymbol{n}$ with $\boldsymbol{x}^i=$ $\left(x_{i 1}, \ldots, x_{i p}\right)^T$ representing the independent variables and $y_i$ representing the dependent variables. As in the standard regression framework, we either assume that the observations are independent or that $y_i s$ are conditionally independent given $x_{i j} s$ [15]. We assume the $x_{i j}$ are standardized such that $\frac{\sum_i x_i}{n}=\mathbf{0}$ and $\frac{\Sigma_i x_{i k}^2}{n}=\mathbf{1}$. Let $\widehat{\boldsymbol{\beta}}=\left(\widehat{\boldsymbol{\beta}}_1, \ldots, \widehat{\boldsymbol{\beta}}_p\right)$, the LASSO estimate $\left(\widehat{\boldsymbol{\beta}}_0, \widehat{\boldsymbol{\beta}}_1\right)$ is defined by:

$\left(\hat{\beta}_0, \hat{\beta}_1\right)={argmin}\left\{\sum_{i=1}^n\left(y_i-\beta_0-\sum_k \beta_k x_{i k}\right)^2\right\}$         (10)

Subject to $\sum_j\left|\beta_j\right| \leq t$.

Here $\boldsymbol{t} \geq \mathbf{0}$ tuning parameter. For any value of $\boldsymbol{t}$, the solution for $\boldsymbol{\beta}_{\mathbf{0}}$ is $\widehat{\boldsymbol{\beta}}_{\mathbf{0}}=\overline{\boldsymbol{y}}$. Without loss of generality, we can assume $\overline{\boldsymbol{y}}=\mathbf{0}$ and thus omit $\boldsymbol{\beta}_{\mathbf{0}}$ [16].

The parameter $\boldsymbol{t} \geq \mathbf{0}$ governs the degree of shrinkage applied to the estimates. Let $\widehat{\boldsymbol{\beta}}_{\boldsymbol{k}}^{\mathbf{0}}$ denote the full least squares estimates, and let $\boldsymbol{t}_{\mathbf{0}}=\sum\left|\widehat{\boldsymbol{\beta}}_{\boldsymbol{k}}^{\mathbf{0}}\right|$. Values $\boldsymbol{t}<\boldsymbol{t}_{\mathbf{0}}$ will lead to shrinkage of the estimates towards zero, with some coefficients potentially being exactly zero.

LASSO is a technique used to address multicollinearity, a problem where independent variables are highly correlated. Introduced by Tibshirani in 1996 [17], LASSO works by reducing the regression coefficients of independent variables that are strongly correlated with the error term to zero or near zero [18]. LASSO estimation is obtained from the following equation:

$\boldsymbol{\beta}^{\text {lasso }}=\boldsymbol{a r g m i n}_{\boldsymbol{\beta}}\left\{\sum_{i=1}^n\left(y_i-\beta_0-\sum_{k=1}^p \beta_k x_{i k}\right)^2+\lambda \sum_{k=1}^p\left|\beta_k\right|\right\}$          (11)

provided that $\sum_{k=1}^p\left|\beta_k\right| \leq t$. The value of t is a tuning parameter that controls the shrinkage of the LASSO coefficient with $t \geq 0$. According to Djuraidah's study [19], LASSO coefficient estimation uses quadratic programming with inequality constraints.

2.7 Parameter estimation of GWR model with LASSO method

Estimation of GWR model parameters with the LASSO method is only the addition of a weighting matrix so that it is obtained:

$\hat{\beta}^{\text {lasso }}={argmin}\left\{\sum_{i=1}^n W\left(u_i, v_i\right)\left(y_i-\beta_{0\left(u_i \cdot v_i\right)}-\sum_{k=1}^p \beta_{k\left(u_i \cdot v_i\right)} x_{i k}+\lambda \sum_{k=1}^p\left|\beta_{k\left(u_i \cdot v_i\right)}\right|\right\}\right)^2$          (12)

Since the equation has an absolute value constraint on the regression coefficient, it causes a nonlinear pattern so that the solution process uses a quadratic program assisted by the LARS algorithm.

2.8 Least Angle Regression

LASSO has been solved using the LARS algorithm [20]. The LARS algorithm is generally as follows:

1) Assume all regression coefficients $\beta_k$ are zero so that $\varepsilon=0$.

2) Select the independent variable that has the highest correlation coefficient with $\varepsilon$.

3) Estimating the coefficient $\beta_k$ for $x_{i k}$ that has the highest correlation with $\varepsilon$.

4) Calculate the residuals $\varepsilon=y-\hat{y}$ with the independent variable $x_k$ included in the model.

5) Calculate the partial correlation between the remaining independent variables and the new $\varepsilon$.

6) Repeating steps 3 to 6 until all independent variables are included in the model and stopping if the correlation $\left(y, x_{i k}\right)=0$.

4.1 Data exploration

As this study involves spatial data, exploratory data analysis was conducted to gain initial insights. The data utilized for this study is the poverty rate data in North Sumatra for the year 2022. The results of this descriptive analysis are presented in Table 2.

Table 2. Data descriptive statistical analysis

4.2 Multiple linear regression with ordinary least square method

Below is a multiple linear regression model utilizing the OLS method to analyze the factors affecting the poverty rate in North Sumatra, as implemented in R Studio:

$\begin{aligned} Y=61.2976- & 2.8220 X_1+0.0683 X_2+1.4975 X_3 +2.9227 X_4+0.0591 X_5+0.0017 X_6\end{aligned}$          (13)

4.3 Spatial heterogeneity test

In Table 3, the p-value of 0.0136 is smaller than the value of α=0.05, so  is rejected. This means that there is spatial heterogeneity in the model, namely the difference in characteristics between observation locations. Therefore. it is necessary to do GWR modeling.

Table 3. Results of spatial heterogeneity testing with Breunsch Pagan

4.4 GWR modeling

Calculating Euclidean distance

The initial step in Geographically Weighted Regression (GWR) modeling involves establishing the geographical coordinates (longitude and latitude) for each regency/city in North Sumatra, representing the observation location of each region. Then calculate the Euclidean distance between  and  using Eq. (5).

The following is an example of calculating the Euclidean distance of several districts:

• Nias

$\begin{gathered}d_{11}=\sqrt{(1.033-1.033)^2+(97.766-97.766)^2}=\sqrt{0+0}=0\end{gathered}$          (14)

• Mandailing Natal

$\begin{gathered}d_{12}=\sqrt{(0.783-1.033)^2+(99.254-97.766)^2} =\sqrt{0.062+2.214}=1.509\end{gathered}$          (15)

The Euclidean distance of each regency/city in North Sumatra is presented in Table 4.

Table 4. Euclidean distance between regencies/cities in North Sumatra

Determining the bandwidth value

The next step is to determine the bandwidth value (h) obtained from the minimum CV value. The weighting function used is the Gaussian Kernel. Determination of bandwidth with the Cross Validation method using R-Studio software. For the minimum CV value and Bandwidth can be seen in Table 5.

Table 5. CV value and optimum bandwidth value

Furthermore, the GWR weighting function is calculated using Eq. (9). With the value of $h=0.6145937$ then obtained:

$\mathrm{W}\left(\mathrm{u}_{\mathrm{i}}, \mathrm{v}_{\mathrm{i}}\right)=\exp \left(-\frac{1}{2}\left(\frac{\mathrm{d}_{\mathrm{i}, \mathrm{i}+1}}{0.6571}\right)^2\right)$          (16)

Next, calculate the weight matrix based on the Fixed Gaussian Kernel weight function, namely:

$\begin{gathered}\mathrm{W}\left(\mathrm{u}_1, \mathrm{v}_1\right)=\exp \left(-\frac{1}{2}\left(\frac{\mathrm{d}_{11}}{\mathrm{~h}}\right)^2\right)=\exp \left(-\frac{1}{2}\left(\frac{0}{0.6571}\right)^2\right)=1\end{gathered}$          (17)

Substitute the weight values from each location into the matrix to obtain the overall W matrix as follows.

$W=\left[\begin{array}{ccccc}1 & 0.0715 & 0.0597 & \cdots & 0.9026 \\ 0.0715 & 1 & 0.5368 & \cdots & 0.0329 \\ 0.0597 & 0.5368 & 1 & \cdots & 0.0425 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0.9026 & 0.3297 & 0.0425 & \cdots & 1\end{array}\right]$          (18)

Table 6. Parameter estimation results of GWR model

GWR parameter estimation

The results of the parameter estimation are listed in Table 6.

GWR model hypothesis test

GWR model fit test. The $\boldsymbol{F}_{ {count }}$ value of GWR with $\boldsymbol{d} \boldsymbol{f}_{\boldsymbol{1}}$=26 and $\boldsymbol{d} \boldsymbol{f}_{\mathbf{2}}$=9.8272 is 4.876 with $\boldsymbol{F}_{ {table }}=\mathbf{2 . 3 5 5 2}, \boldsymbol{p}-$ value $=0.0065$. Because the $\boldsymbol{p}-$ value $<\boldsymbol{\alpha}=\mathbf{0 . 0 5}$ and the value of $\boldsymbol{F}_{\text {count }}>\boldsymbol{F}_{ {table }}$, the decision is made to reject $\boldsymbol{H}_{\mathbf{0}}$, this indicates a substantial difference between the global multiple linear regression model and the local Geographically Weighted Regression (GWR) model.

GWR model parameter significance test. Parameter significance test is determined by testing the parameters partially. The $\boldsymbol{t}_{ {table }}=2.0595$ was obtained. If the $\left|\boldsymbol{t}_{ {count }}\right|$ value of each district/city is greater than 2.0595. This indicates that the independent variable affects the dependent variable.

Table 7. Summary of t test on each independent variable

Based on the results of the t-test as shown in Table 7, the independent variables demonstrate varying significance across different regions: Variable $X_1$ is significant in Gunungsitoli, indicating that the factors represented by $X_1$ have a noticeable impact on the outcomes in this area. This may be due to unique local conditions or specific policies that amplify the influence of this variable. Variable $X_2$ shows significance in North Nias, which could be attributed to demographic or economic characteristics in the region that align with the aspects measured by $X_2$. Variable $X_3$ is significant in several regions, including Nias, South Tapanuli, South Nias, North Nias, West Nias, Sibolga, Padang Sidempuan, and Gunungsitoli. This widespread influence suggests that $X_3$ may be related to broader regional factors affecting many areas, such as development policies or widespread environmental issues. Variable $X_5$ shows significance in Nias, North Nias, South Nias, West Nias, and Gunungsitoli, indicating that $X_5$ has specific relevance in these areas, possibly linked to social structures or relevant infrastructure impacting this variable. These results highlight the importance of contextual analysis when assessing the impact of independent variables, as their influence can vary significantly based on geographic location and the specific characteristics of each region.

4.5 Detection of local multicollinearity

The next step is to assess the presence of multicollinearity in the GWR model, which can be determined by examining the local VIF (Variance Inflation Factor) values. From the table, it is evident that some VIF values exceed 10. For instance, the HDI (Human Development Index) variable in Nias Regency has a VIF value of 133.31897.

This indicates that the standard error of the estimated coefficient for the HDI variable will be significantly higher 11. $5463(\sqrt{133.3189})$ times than it would be if there were no correlation with other independent variables. A summary of VIF values in the research data is presented in Table 8.

Table 8. Summary of local multicollinearity

4.6 GWR modeling with LASSO

Since the GWR model exhibits local multicollinearity, it can be addressed using the LASSO method. The following steps will be presented. GWR Modeling with LASSO in Nias Regency from the LARS algorithm obtained Table 9.

Table 9. Stages in the Nias Regency LARS algorithm

The selection of the best LASSO model is done by the 5-fold Cross Validation method with fraction mode, which is by calculating the cross-validation value for each step with one variable entered into the model.

Based on Figure 1, it is shown that there are 100 steps generated by mode fraction. Therefore, the resulting $s$ value is also 100 different values. Where for the best minimum $s$ value is when the cross-validation value is minimum. In this model, the smallest cross validation value is at step 31, namely with a cross validation value of 0.0322 with an $s$ value of 0.303. Next, the $s$ value obtained from the fraction mode will be compared with the $s$ value of the LARS algorithm.

1.png

Figure 1. CV value using mode fraction Nias Regency

The coefficient value at each stage will be compared with the sum of the OLS method coefficients and the $s$ value obtained in the LARS algorithm will be adjusted to the s value in the fraction method and the $s$ value is the shrinkage value.

Because the value of s in the fraction mode is 0.303 and in the LARS method the value of $s$ that is close to 0.303 is 0.288, namely at the $\mathbf{6}^{\text {th }}$ stage in Table 10.

Table 10. Comparison of LASSO coefficient values with OLS coefficients Nias Regency

Note: The final model is obtained:

$Y=0.0962 X_2+0.1101 X_3+0.1125 X_4-0.7019 X_5-1.1092 X_6$

The value of $s$ is 0.288 with an R-squared value of 0.9403. This indicates that in Nias Regency, 94.03% of the variation in the poverty rate is explained by the independent variables in the model, while the remaining 5.97% is influenced by factors not included in the model.

GWR Modeling with LASSO in Medan City. The outcomes of the LARS algorithm are shown in Table 11.

Table 11. Stages in LARS algorithm for Medan City

The selection of the best LASSO model is done with the n-fold Cross Validation process with mode fraction, namely by calculating the cross validation value for each step with one variable entered into the model.

Based on Figure 2, it is shown that there are 100 steps generated by mode fraction. Therefore, the resulting s value is also 100 different values. Where for the best minimum s value is when the cross validation value is minimum. In this model, the smallest cross validation value is at step 14, namely with a cross validation value of 0.0215 with an s value of 0.131. Next, the s value obtained from the fraction mode will be compared with the s value of the LARS algorithm.

2.png

Figure 2. CV value using mode fraction Medan City

The coefficient value at each stage will be compared with the sum of the OLS method coefficients and the $s$ value obtained in the LARS algorithm will be adjusted to the $s$ value in the fraction method and the $s$ value is the shrinkage value.

Because the value of $s$ in the fraction mode is 0.131 and in the LARS method the value of s that is close to 0.131 is 0.1376, namely at the $6^{t h}$ stage presented in Table 12.

Table 12. Comparison of LASSO coefficient value with OLS coefficient of Medan City

Note: The final model is obtained:

$Y=-0.2255 X_2-0.3152 X_3-0.1936 X_4+0.1598 X_5+0.1206 X_6$

The value of $s$ is 0.137 with an R-squared value of 0.5143. This indicates that in Medan, 51.43% of the variation in the poverty rate is explained by the independent variables in the model, while the remaining 48.57% is affected by other factors outside the model.

Spatial patterns in parameter coefficient measures regarding geographic heterogeneity in the determinants of poverty focus on understanding how the variables that influence poverty differ across locations or regions. LASSO helps reduce model complexity and automatically selects the most significant variables. This is particularly relevant in the context of Geographically Weighted Regression (GWR), where local multicollinearity can influence the results. LASSO's advantage in identifying and measuring differences in poverty factors across regions is that it takes into account spatial influences in parameter estimation, thereby providing a better understanding of how poverty determinant variables behave at the local level.

Additionally, LASSO automatically simplifies the model by shrinking the coefficients of less significant variables to zero, increasing interpretability. From a statistical perspective, in Nias Regency, the LASSO model produces an R-squared of 0.9403, which means that 94.03% of the variation in poverty levels can be explained by independent variables. In contrast, in Medan City, the R-squared was recorded at 0.5143, indicating 51.43% of the variation was explained by the model. AIC is also used to evaluate model quality, where lower values indicate a good balance between goodness of fit and complexity. With its computational efficiency, LASSO not only fixes technical problems such as multicollinearity but also provides more precise insights into the determinants of poverty.

4.7 Limitations and future research

In this analysis, there are several data limitations and model assumptions that can affect the accuracy of the results. First, data limitations include temporal coverage limited to one year (2022), which may not reflect temporal variations in poverty trends. Therefore, the use of data from multiple years can provide more in-depth longitudinal analysis. In addition, the variables currently selected only include life expectancy, open unemployment, life expectancy, average years of schooling, and per capita income. Exploring additional socio-economic and environmental variables can provide a more comprehensive understanding. In terms of model assumptions, the spatial independence assumption may not fully hold across regions, so integrating spatial autocorrelation measures can improve model accuracy.

Although LASSO helps overcome multicollinearity, interactions between variables may not have been fully detected, so adding interaction terms or non-linear models can provide deeper insights. From a methodological perspective, experiments with different kernel functions other than Gaussian can still help to better capture spatial dependencies. Validating the model with out-of-sample data is also important to test the model's predictive power on unseen data. Further research should consider geographic differences in factors influencing poverty to guide more targeted policy interventions, as well as use causal inference techniques to understand the relationships between variables and poverty outcomes.