SAR probit regression modeling using the fisher scoring approach: A case study of poverty levels on the island of Sumatra

The initial foundation of probit regression was introduced by Bliss in 1934, offering a nonlinear analytical framework to connect predictors with binary response variables. It is a nonlinear analytical instrument for evaluating the relationship between predictor variables (numerical, categorical, or mixed) and a binary categorical response variable. In this model, a value of 1 represents the existence of a specific attribute, whereas a value of 0 indicates the absence of that attribute. Empirically, categorical response data in this model often exhibit the influence of adjacent regions, known as spatial autocorrelation (Anselin, 1988). Spatial autocorrelation refers to a phenomenon in which observations at a specific location tend to depend on or are influenced by the values in the surrounding geographical areas. If this geographical relationship is ignored in a standard probit model specification, the resulting parameters will be both biased and inconsistent. This bias phenomenon refers to the estimator's inability to approach the true population parameter value owing to unaccommodated distortions in spatial patterns.
The main challenge in SAR Probit modeling lies in its likelihood function, which does not have a closed-form solution, thus requiring a numerical optimization approach. Standard approaches, such as the Newton-Raphson algorithm, iteratively use the Hessian matrix (second derivative), which often experiences convergence failure or matrix instability when dealing with high-dimensional spatial data. To overcome this limitation, this study applies the Fisher Scoring numerical iteration. This algorithm modifies the Newton-Raphson approach by substituting the Hessian matrix with the Fisher information matrix. Mathematically, this substitution guarantees that the matrix is always positive definite, resulting in more stable convergence and more consistent parameter estimates, particularly when handling regional poverty data with high spatial variations (McMillen, 1992).
Parameter estimation plays a crucial role in obtaining unbiased estimated values and determining the covariance matrix of the model. In the spatial probit regression framework, this parameter estimation can be implemented through the Maximum Likelihood Estimation (MLE) method integrated with the Fisher scoring numerical iteration scheme. This iterative procedure is essential because the equation function in the model does not have an analytical solution or is non-closed when estimated with MLE. The Fisher scoring technique was developed as a variant of the Newton-Raphson method and is specifically designed to find numerical solutions for complex maximum likelihood functions. Unlike the Newton-Raphson algorithm, which utilizes a second-derivative-based Hessian matrix and often faces constraints in achieving convergence, Fisher scoring replaces this component with the Fisher information matrix to ensure the stability of the estimation process (Schworer & Hovey, 2004).
Therefore, this study aims to analyze parameter estimation in the Spatial Autoregressive (SAR) probit model using the Fisher Scoring algorithm on the 2022 Sumatra poverty data. The main contributions of this study are twofold. First, methodologically, it demonstrates the efficiency and stability of the Fisher Scoring algorithm in handling complex spatial dependencies compared with the Recursive Importance Sampling (RIS) simulator. Second, empirically, it provides novel insights into the spatial patterns of poverty in Sumatra and offers targeted policy recommendations. By proving that poverty in a region is significantly influenced by neighboring areas, this study suggests that poverty alleviation strategies in Sumatra must be conducted through synchronized cross-regional policies rather than isolated local efforts.

Before delving into the mathematical formulations, it is essential to understand the conceptual transitions that build the Spatial Probit model. Standard probit regression is used to evaluate binary outcomes, such as a region having either 'high' or 'low' poverty based on various predictors. However, socio-economic phenomena such as poverty are rarely isolated; a region's economic status often overflows into neighboring areas. This geographical relationship was captured by spatial regression. By integrating these two concepts, spatial probit regression is formed. This advanced framework not only models the binary categorical response but also intuitively accommodates the spatial influence of adjacent regions, ensuring that the estimated parameters remain unbiased and consistent across the study area.

2.1. Probit Regression
Referring to the thought of (Greene, 2008), the binomial probit regression framework begins with the formation of a latent variable y_i^* formulated through the following linear relationship:

\[
y_i^{*} = \boldsymbol{\beta}^{\top}\mathbf{x}_i + \epsilon_i \tag{1}
\]

Here, y_i^* is a latent response variable influenced by a vector of predictor variables x_i and parameter coefficients β. The error component ϵ_i is assumed to follow a normal distribution with a mean of zero and a variance of σ^2. The probability density function (PDF) of the variable y_i^* is as follows (Greene, 2008):

\begin{equation}
f(y_i^{*}) = \frac{1}{\sqrt{2\pi\sigma^{2}}}
\exp\left(-\frac{1}{2\sigma^{2}}\left(y_i^{*}-\boldsymbol{\beta}^{\top}\mathbf{x}_{i}\right)^{2}\right)
\end{equation}

The determination of the response variable category y_i is based on a threshold value δ, where a value of 1 is assigned if the latent variable exceeds that threshold. Based on this assumption, the probability of an event (y_i=1) occurring can be calculated using the standard normal cumulative distribution function Φ as follows:

\[
\begin{aligned}
P(y_i = 0) &= P(y_i^{*} < \delta) \\
&= P(\boldsymbol{\beta}^{\top}\mathbf{x}_i + \epsilon_i < \delta) \\
&= P(\epsilon_i < \delta - \boldsymbol{\beta}^{\top}\mathbf{x}_i) \\
&= \Phi(\delta - \boldsymbol{\beta}^{\top}\mathbf{x}_i)
\end{aligned}
\tag{3}
\]

\[
\begin{aligned}
P(y_i = 1) &= P(y_i^{*} \geq \delta) \\
&= 1 - P(y_i^{*} < \delta) \\
&= 1 - P(\boldsymbol{\beta}^{\top}\mathbf{x}_i + \epsilon_i < \delta) \\
&= 1 - P(\epsilon_i < \delta - \boldsymbol{\beta}^{\top}\mathbf{x}_i) \\
&= 1 - \Phi(\delta - \boldsymbol{\beta}^{\top}\mathbf{x}_i)
\end{aligned}
\tag{4}
\]

2.2. Spatial Regression
The spatial regression approach for cross-sectional data was pioneered by Anselin (1988), one of which is through the Spatial Autoregressive (SAR) structure. For categorical response variables, the SAR model can be expressed as:

\[
y^{*} = \rho W y^{*} + X\beta + \epsilon
\tag{5}
\]

Equation (5) involves the parameter ρ, which indicates the spatial lag coefficient, as well as the weight matrix W, which represents the inter-regional interactions determined through the Queen Contiguity method. Queen Contiguity, or edge-and-corner contiguity, defines w_{ii^* } = 1 for entities that share an edge or a vertex with the region of interest, and w_(ii^* )=0 for all other regions. For region (3), the values obtained are w_32=1, w_34=1, w_35=1, while the remaining values are equal to zero. Figure 1 illustrates an example of five regions as they appear on a map.

Figure 1. Illustration of Queen Contiguity

The selection of this edge-and-corner contiguity weight (Queen) is based on the geographical characteristics of the administrative boundaries of the districts/cities on Sumatra Island, which form highly irregular polygons. Consequently, even the slightest intersection at the boundary corners is accounted for as a spatial neighbor that potentially exerts a poverty spillover effect. It is important to note that the diagonal elements of matrix W must be zero to prevent self-correlation.
In spatial analysis, dependence between locations generally appears in two main formats: spatial lag and spatial errors. The spatial lag phenomenon indicates that the conditions in a region are directly influenced by its neighboring regions. Spatial error refers to the existence of inter-regional linkages that are not only limited to the response variable but also include disturbances in predictor variables (Anselin, 1988). To validate the presence of spatial lag interaction effects in the model, researchers can apply the Lagrange Multiplier (LM) test. This procedure is carried out by testing the significance of the parameter ρ through the following hypotheses:  
H_0:ρ=0 (Assumption of no spatial lag dependence)  
H_1:ρ≠0 (Indication of spatial lag dependence)  
The test statistic used is as shown in Equation (6):

\[
LM = E^{-1} \left\{ (R_y)^2 Tr_{22} - 2R_y R_e Tr_{12} + (R_e)^2 (D + Tr_{11}) \right\} \sim \chi^2_{(m)}
\tag{6}
\]

where:
\[
\begin{aligned}
R_y &= e^{\top} W_1 y \,/\, \sigma^2 \\
R_e &= e^{\top} W_2 e \,/\, \sigma^2 \\
M &= I - X(X^{\top}X)^{-1}X^{\top} \\
\mathrm{Tr}_{k k^{*}} &= \mathrm{tr}\left\{ W_k W_{k^{*}} + W_k^{\top} W_{k^{*}} \right\} \\
D &= \sigma^{-2}(W_1 X\beta)^{\top} M (W_1 X\beta) \\
E &= (D + T_{11})T_{22} - (T_{12})^2
\end{aligned}
\]

The null hypothesis will be rejected if the LM statistic exceeds the critical value on the χ^2 distribution with m degrees of freedom. This condition indicates the influence of spatial lag dependence on the model. Specifically for the SAR model structure, the value of m is set to one because it only uses the coefficient ρ as a single spatial parameter.

2.3. Spatial Probit Regression
The synergy between the probit regression framework and spatial analysis methodology yields a spatial probit model. The development of this model is rooted in the research of McMillen (1992), who proposed the use of the EM algorithm to produce consistent parameters under autocorrelated conditions. In the Spatial Autoregressive (SAR) scheme, the specification of the spatial probit model can generally be expressed as follows:

\[
y^{*} = \rho W y^{*} + X\beta + \epsilon
\tag{7}
\]

where
\[
\epsilon \sim N(0, \sigma^2 I_n)
\]

In this context, the vector y^* represents the latent response variable, while inter-regional interactions are captured through the spatial parameter ρ and the weight matrix W. For ease of estimation, this model can be transformed into a reduced form, as follows:

\[
y^{*} = (I - \rho W)^{-1} X\beta + (I - \rho W)^{-1} \epsilon
\tag{8}
\]

let
\[
(I - \rho W)^{-1} = \Delta
\quad \text{and} \quad
(I - \rho W)^{-1}\epsilon = u
\]
so that y^* forms the equation.

\[
\begin{aligned}
y^{*} &= \Delta X\beta + u, \\
u &\sim \mathrm{MVN}(0, \mathbf{T})
\end{aligned}
\tag{9}
\]

The variable y_i^* has a binary category defined as the variable y_i through the following equation:

\[
y_i =
\begin{cases}
0, & \text{jika } y_i^{*} < \delta \\
1, & \text{jika } y_i^{*} \geq \delta
\end{cases}
\tag{10}
\]

Then, the probability from Equation (10) can be written as follows:

\[
P(y_i = 0) = P(y_i^{*} < \delta)
= \Phi \left( \frac{\delta - [\Delta X\beta]_i}{\sqrt{\zeta_{ii}}} \right)
\tag{11}
\]

\[
P(y_i = 1) = P(y_i^{*} \geq \delta)
= 1 - \Phi \left( \frac{\delta - [\Delta X\beta]_i}{\sqrt{\zeta_{ii}}} \right)
\tag{12}
\]

From equations (11) and (12), the variable y_i will be obtained:

\[
y_i =
\begin{cases}
0, & \text{jika } u_i < \dfrac{\delta - [\Delta X\beta]_i}{\sqrt{\zeta_{ii}}} \\
1, & \text{jika } u_i \geq \dfrac{\delta - [\Gamma X\beta]_i}{\sqrt{\zeta_{ii}}}
\end{cases}
\tag{13}
\]

The variable Y_i in equation (13) is Bernoulli distributed with the following probability function:

\[
P(Y_i = y_i)=
\begin{cases}
p(x_i)^{y_i}\bigl(1-p(x_i)\bigr)^{1-y_i}, & y=0,1;\ 0<p(x_i)<1,\\
0, & \text{lainnya.}
\end{cases}
\tag{14}
\]

The mean and variance of Y_i are as follows:

\[
\begin{aligned}
E(Y_i) &= p(x_i), \\
\operatorname{Var}(Y_i) &= p(x_i)\bigl(1-p(x_i)\bigr)
\end{aligned}
\tag{15}
\]

where p(x_i ) is the probability of success of P(y_i=1). Since the variable u is assumed to have a multivariate normal distribution with a mean of 0 and a variance of T, the probabilities for P(y_i=1) and P(y_i=0) based on equations (11) and (12) are as follows:

\[
P(y_i = 0 \mid x_i) = \bigl(1 - p(x_i)\bigr)
= \Phi \left( \frac{\delta - [\Lambda X\beta]_i}{\sqrt{\zeta_{ii}}} \right)
\tag{16}
\]

\[
P(y_i = 1 \mid x_i) = p(x_i)
= 1 - \Phi \left( \frac{\delta - [\Lambda X\beta]_i}{\sqrt{\zeta_{ii}}} \right)
\tag{17}
\]

2.4. Fisher Scoring Parameter Estimation
In cases where the Maximum Likelihood Estimation (MLE) yields non-closed-form equations, an iterative procedure is required to find numerical solutions. One effective technique is the Fisher scoring algorithm, which modifies the Newton-Raphson method. According to Schworer and Hovey (2004), the iteration structure in this method is formulated to update parameter estimates until a certain convergence point is reached. The stability of this algorithm is guaranteed by the use of the Fisher information matrix in its updating process, which is mathematically expressed as:

\[
\boldsymbol{\beta}^{*(t+1)}=\boldsymbol{\beta}^{*(t)}+\left(\widetilde{\mathcal{I}}\!\left(\boldsymbol{\beta}^{*(t)}\right)\right)^{-1} S\!\left(\boldsymbol{\beta}^{*(t)}\right)
\]

where β^* is the vector of model parameter estimators, namely
\[
\beta^{*} = [\beta,\rho]^{T}
\]
when applied to the spatial probit regression model. t is the iteration level.
\[
\widetilde{\mathcal{I}}(\beta^{*})
\]
is the Fisher information matrix of β^*, namely:
\[
\widetilde{\mathcal{I}}(\beta,\rho)
=
- E \left[
\frac{\partial^{2}\bigl(l(\beta,\rho)\bigr)}
{\partial(\beta,\rho)\,\partial(\beta,\rho)^{T}}
\right]
\]

Meanwhile, S(β^* ) is the gradient vector to be estimated, having the following form:

\[
S(\beta,\rho)=\frac{\partial\bigl(l(\beta,\rho)\bigr)}{\partial(\beta,\rho)}
\]

The Fisher scoring iteration process stops if the convergence condition is met, namely

\[
\left\| \beta^{*(t+1)} - \beta^{*(t)} \right\| < \Delta
\]

where Δ is a verll positive real value. If the criteria are met, the iteration process stops at the T-th iteration, so the parameter estimator is
\[
\widehat{\beta}^{*} = \beta^{*(T)}
\]
After obtaining the value of the parameter estimator, parameter testing on the spatial probit model can be continued.

2.5. Simultaneous Parameter Testing
The significance of all predictor variables in the spatial probit model was evaluated using the Maximum Likelihood Ratio Test (MLRT). The operational basis of this procedure is to contrast the capability of the model that has been equipped with predictors against a model without independent variables (under the null hypothesis). Statistical decisions are made by comparing the G_SP^2 value against the χ^2 distribution at the appropriate degrees of freedom. The hypothesis framework for simultaneously evaluating the parameters is formulated as follows:
\[
H_0:\ \beta_1=\beta_2=\cdots=\beta_p=0
\]

\[
H_1:\ \text{at least one }\beta_j\neq 0,\ \text{for }j=1,2,\ldots,p
\]
Suppose L(Ω ̂_SP ) is the likelihood function under the population with
\[
\widehat{\Omega}_{SP} = \{\rho,\beta\},\, L(\widehat{\omega}_{SP})
\]
is the likelihood function under H_0 with
\[
\widehat{\omega}_{SP} = \{\rho_{\omega 0}, \beta_0\}
\quad \text{and} \quad
G_{SP}^{2}
\]
is the test statistic used for simultaneous parameter significance testing, also called the likelihood ratio test. Then G_SP^2 can be written in equations (18) and (19).

\[
V=\left[\frac{L(\widehat{\omega}_{SP})}{L(\widehat{\Omega}_{SP})}\right]
\tag{18}
\]

\[
G_{SP}^{2}=-2\ln\!\left[V_{SP}\right]
=2\left[\ln L(\widehat{\Omega}_{SP})-\ln L(\widehat{\omega}_{SP})\right]
\tag{19}
\]

The decision to reject the null hypothesis (H_0 ) is taken if the statistical value of G_SP^2 is greater than the critical value of χ_df^2. The degrees of freedom (df) magnitude is obtained from the difference between the number of parameters in the population and the total parameters under the H_0 condition. The comparison value of χ_df^2 can be seen through the Chi-Square distribution table reference.

2.6. Partial Testing
After the overall model significance was proven, the next step was to dissect the contribution of each individual predictor. This procedure was performed using a partial test based on the Wald test statistic. The Z_hitung statistic is computed based on the comparison between the estimated parameter value and its standard error, which is then compared with the critical value in the standard normal distribution table. The hypotheses for conducting the partial tests were as follows:
Η_0:β_j=0, for j=1,2,…,p.
Η_1:β_j≠0.
The test statistic used in the partial model parameter testing is as in Equation (20) below:

\[
Z_{\mathit{hitung}}=\frac{\widehat{\beta}_j}{\mathrm{SE}(\widehat{\beta}_j)}
\tag{20}
\]

\[
\text{where } \widehat{\beta}_j \text{ is the parameter estimator of } \beta_j
\text{ and } \mathrm{SE}(\widehat{\beta}_j)
\text{ is } \sqrt{\widehat{\mathrm{var}}(\widehat{\beta}_j)},
\text{ with } \widehat{\mathrm{var}}(\widehat{\beta}_j)
\text{ being the } (j+1)\text{-th diagonal of }
\widehat{\mathrm{Cov}}(\widehat{\boldsymbol{\beta}}^{*})
\cong n^{-1}\widetilde{\mathcal{I}}^{-1}(\widehat{\boldsymbol{\beta}}^{*}).
\]
\[
\text{The decision is to reject } H_0
\text{ if } Z_{\mathit{hitung}} < -Z_{\alpha/2}
\text{ or } Z_{\mathit{hitung}} > Z_{\alpha/2},
\text{ where the value of } Z_{\alpha/2}
\text{ can be obtained from the standard normal table.}
\]

3.1. Research Data and Variables
This study utilizes secondary data from Badan Pusat Statistik (2022), covering 131 observation units of regencies/cities on Sumatra Island. The analyzed response variable was the binary poverty status with a threshold value of 9.54%. The factors tested for their influence included GRDP growth, unemployment rate, per capita expenditure, and expected years of schooling (see Table 1).

Table 1. Research Variables

3.2. Research Stages  
The research process was divided into several interconnected stages. The complete sequence for each stage is presented in the flowchart below (see Figure 2):

Figure 2. Research Phase Flowchart

4.1. Overview of Poverty Rates in Regencies/Cities in Sumatra Island
A visualization of the distribution of poverty rates at the regency/city level on Sumatra Island for 2022 is presented in Figure 3. In the map, the low poverty rate classification is represented by the yellow zone, while areas with high poverty rates are marked in brown. Based on visual observation, there is a pattern where geographically adjacent regions show similarities in the degree of poverty, which explicitly indicates the presence of spatial dependency in poverty data in the Sumatra region.

Figure 3. Poverty Rates by Regency/City on the Island of Sumatra in 2022

4.2. Pre-Modeling Analysis
The pre-estimation stage in the SAR Probit model requires verification of classic assumptions to ensure that the resulting model has high validity. One crucial aspect is ensuring the absence of multicollinearity, which is a phenomenon in which predictor variables correlate significantly with each other. As a preventive measure in modeling the poverty rate, early detection is performed by evaluating the Variance Inflation Factors (VIF) values. The following are the details of the VIF parameters for all the independent variables analyzed in this study.

Table 2. Predictor Variable Multicollinearity

Source: RStudio Software Processing Results

The data presented in Table 2 confirm that all Variance Inflation Factor (VIF) values are below the threshold of 10. These findings provide empirical evidence that there are no multicollinearity issues among the independent variables in this case study. With this assumption fulfilled, all predictor variables were eligible for integration into the subsequent statistical analysis stages.  
The next step in this research procedure was to evaluate the presence of autocorrelation in the dependent variable through spatial dependency testing. The identification of inter-regional interactions is implemented using the Lagrange Multiplier (LM-Lag) test. The hypothesis formulation used to detect this spatial effect specifically refers to the theoretical framework.
Η_0∶ρ=0 (no spatial lag dependency)
Η_1∶ρ≠0 (spatial lag dependency exists)
Referring to the statistical calculation results using the formula with a 5% error rate, an LM lag value of 19.51115 was obtained. This figure is greater than the Chi-square table value (3.8415) for one degree of freedom, so it is decided to reject H_0. This indicates strong evidence of inter-regional dependence in poverty patterns on Sumatra Island in 2022. Thus, the analysis process can proceed to the core modeling stage.

4.3. Poverty Rate Modeling with SAR Probit Method  
The application of the spatial probit model in this study is intended to identify various variables that significantly contribute to poverty. The analysis process began with parameter estimation using the Fisher Scoring algorithm. This technique was chosen to overcome the complexity of non-closed-form functions in the Maximum Likelihood procedure. The parameter estimation results are as follows (see Table 3):

Table 3. Parameter Estimation Results Using Fisher Scoring Iteration

Source: RStudio Software Processing Results

From a spatial planning policy perspective, the presence of a significant spatial dependence effect indicates that poverty in Sumatra is closely linked. Policy interventions, such as the development of educational infrastructure, job creation, or economic stimulus, should not be implemented in isolation in each regency. Instead, local governments in Sumatra must adopt a coordinated, cross-regional approach, as economic improvements or setbacks in one district significantly impact the surrounding areas.
Based on the parameter estimation results using the Fisher Scoring approach, the SAR Probit regression model is as follows:

\[
[\Lambda X\beta]_i
=
3,384
+
2,535\left(\sum_{i=1,\; i^{*}\neq i}^{n} w_{ij}y_i^{*}\right)
-0,159X_{1i}
+0,114X_{2i}
-0,1677X_{3i}
-0,205X_{4i}
\]

The analysis process then continues by verifying the parameters in the SAR Probit model through simultaneous and partial testing approaches. In the spatial probit regression model, simultaneous testing is conducted to evaluate the collective effect of the coefficient vector β. This joint significance test utilizes the likelihood ratio test method based on the following hypothesis:
Η_0:β_1=β_2=β_3=β_4=0,
Η_1: at least one β_j≠0, for j=1,2,3,4.
Based on the collective parameter significance test results, a G_SP^2 statistical value of 44.30691 is obtained. This value is proven to be greater than the critical Chi-Square table value at a significance level of α=0.05 with 4 degrees of freedom, which is 9.487729, so the statistical decision taken is to reject the null hypothesis (H_0 ). This finding indicates that in the spatial probit regression model with a 95% confidence level, at least one independent variable has a tangible effect on the model as a whole. Based on these results, the analysis proceeded to the individual parameter testing stage (partial test) to identify which predictor variables specifically contributed significantly to the response variable.
The decision to reject H_0 in the simultaneous test requires a more in-depth analysis through partial hypothesis testing of the spatial probit model parameters. The focus of this test was to dissect the level of significance of the contribution of each predictor variable to the response variable separately. In its implementation, the Wald test statistic is used as the main instrument to verify the individual influence of each variable in the model.
Η_0:β_j=0,
Η_1: β_j≠0, for j=1,2,3,4.
The details of the individual parameter evaluation results with a significance level of α=0.05 are presented in Table 4. Based on these data, it is identified that the GRDP growth rate variable (X_1 ), open unemployment rate (X_2 ), average per capita expenditure (X_3 ), and expected years of schooling (X_4 ) have a partially significant effect on the poverty classification at the regency/city level in Sumatra Island. This finding confirms that these four economic and social indicators are important predictors of mapping poverty in the region.

Table 4. Partial Parameter Testing of Spatial Probit Model

Source: RStudio Software Processing Results

Table 4 shows that, in essence, the negative impact of GRDP growth rates and average per capita expenditure indicates that regions with higher economic growth and purchasing power tend to have a lower probability of experiencing high poverty rates, as these factors reflect a higher standard of living and the region’s economic capacity. The positive impact of the open unemployment rate implies that a lack of job opportunities directly increases the risk of poverty owing to the loss of a stable source of income. Meanwhile, an increase in the expected years of education significantly reduces poverty, highlighting the critical role of education in enhancing human capital and productivity.
After the model structure is successfully formed, the next crucial stage involves calculating the predicted probability. The estimated probability for observation y_i=1 in modeling the poverty rate across regencies/cities in Sumatra Island is formulated as follows:

\[
\widehat{P}(y_i = 1 \mid X)
=
\widehat{P}(y_i^{*} \geq 0 \mid X)
=
1-\Phi\left(
\frac{-[\Lambda X\widehat{\beta}]_i}{\sqrt{\zeta_{ii}}}
\right)
=
1-\Phi\left(
\frac{-\left[(I-2{,}5349W)^{-1}X\widehat{\beta}\right]_i}{\sqrt{\zeta_{ii}}}
\right)
\]

As an illustration of the model's application, the calculation of the predicted probability value for Medan City, which occupies observation index i=49, is carried out. The magnitude of the probability y_i=1 in modeling the poverty rate in the region is determined through the following mathematical procedure.

\[
[\Lambda X\widehat{\beta}]_{49}
=
3,38
+
2,54\left(
\sum_{i=1,\; i^{*}\neq i}^{131} w_{49i^{*}}\,y_{49}^{*}
\right)
-0,16X_{1;49}
+0,11X_{2;49}
-0,17X_{3;49}
-0,22X_{4;49}
\]

Therefore, if substituted into Equation (17), the model becomes

\[
P(y_{49}=1 \mid X_{49})
=
1-\Phi\left(
\frac{-[\Lambda X\widehat{\beta}]_{49}}{\sqrt{\zeta_{49,49}}}
\right)
=
0,096536984
\]

The estimated predicted probability result for Medan City shows a value of 0.096536984, which implies that the probability of the region being in the high poverty category is only 9.65%, while the probability of remaining in the low poverty category is 90.35%. This finding indicates that the poverty rate in Medan City is influenced by internal indicators such as the economic growth rate, open unemployment rate, per capita income, and expected years of schooling. The existence of spatial linkages with surrounding regions, especially Deli Serdang Regency, which is also predicted to have a low poverty rate, makes a tangible contribution to Medan City's economic status. This reinforces the strong geographical dependence between the two bordering regions in the context of poverty distribution.  
A visualization of the distribution of the predicted poverty results estimated using the SAR Probit model is presented in Figure 4. The findings from this modeling show that there are 66 regencies/cities included in the high poverty rate classification, while the other 65 regions are predicted to be in the low category.

Figure 4. Map of Actual and Predicted Poverty Levels from The Spatial Probit Model

A more detailed explanation of the regencies/cities between the actual and prediction data is shown in Table 5.

Table 5. Results of Actual and Prediction Data Grouping by Regency/City

4.4. Evaluation and Model Comparison
After conducting spatial probit modeling using the Fisher scoring approach, a comparison was made with the spatial probit model using the RIS simulator approach. The models were compared based on their accuracy, sensitivity, and specificity. The comparison is also made by examining which SSE value is smaller for each model.  
The results of comparing the classification capabilities of the SAR Probit model estimated using the Fisher Scoring approach and the Recursive Importance Sampling (RIS) simulator approach are presented in the following table.

Table 6. Comparison Results of Fisher Scoring Iteration Approach with RIS Simulator

Source: RStudio Software Processing Results

The comparison results in Table 6 show that the Fisher Scoring approach is superior in terms of total accuracy and sensitivity. The much lower SSE value (15.87) indicates that this model has better estimation precision than the RIS Simulator in modeling poverty in Sumatra. Therefore, it can be concluded that the spatial probit model with the Fisher scoring approach is more effective and accurate in modeling poverty rates, providing predictions closer to the true values and producing better statistical results.
The pattern revealed by the SAR probit model should not be read only as a statistical dependency but as evidence that poverty in Sumatra is territorially structured. The coexistence of high-poverty clusters and low-poverty clusters implies that neighboring labor markets, commodity chains, and public-service environments may transmit deprivation across administrative borders. Similar interpretations are common in spatial growth and convergence research, where adjacent regions often move together because shocks and opportunities diffuse geographically (Rey & Montouri, 1999; Fingleton & López-Bazo, 2006). In practical terms, the findings suggest that district-level poverty programs will be less effective when they are implemented without coordination with surrounding jurisdictions.
This result is also consistent with the broader logic of spatial discrete-choice modeling. When poverty status is coded as a binary outcome, the probability that a district enters the high-poverty category is shaped not only by its own socioeconomic conditions but also by the configuration of neighboring areas. Earlier work on neighborhood influence and spatial probit estimation reached a similar conclusion: local choices and outcomes can become interdependent once interaction effects are admitted into the model (Case, 1992; LeSage, 2000). Therefore, the significant spatial parameter in this study should be understood as substantive spillover rather than a mere nuisance term.
From a methodological standpoint, the superior performance of the Fisher Scoring approach over the RIS simulator is important because it shows that estimation stability can materially affect substantive inference. A method that converges more reliably and yields smaller prediction error is especially valuable when researchers work with regional data characterized by uneven distributions and nontrivial spatial dependence. This interpretation aligns with the wider spatial econometric literature, which emphasizes that implementation details, optimization routines, and numerical stability shape model quality as much as theoretical specification does (Elhorst, 2010; Bivand & Piras, 2015). In this study, the higher sensitivity of Fisher Scoring is particularly useful because it improves the identification of districts that are genuinely vulnerable to high poverty.
At the same time, the current specification should be read as a parsimonious starting point rather than the final word on spatial poverty modeling in Sumatra. Recent spatial-econometric work has shown that alternative formulations, including models that more explicitly trace local spillovers or simplify likelihood computation, can provide complementary insights depending on the theoretical objective and data structure (LeSage & Pace, 2007; Halleck Vega & Elhorst, 2015). Future research could therefore compare SAR probit estimates with SLX-style or MESS-based formulations, especially when the aim is to separate direct local effects from broader neighborhood transmission mechanisms.
The policy implication is that poverty reduction on Sumatra Island should combine place-based targeting with intergovernmental coordination. Improvements in education, labor absorption, and household expenditure in one district may generate indirect benefits for adjacent districts, while persistent weakness in one area can slow progress in nearby locations. This reading is consistent with foundational discussions of indirect effects and spatial spillovers in regional econometrics, as well as later work on the large-sample behavior of alternative spatial specifications (LeSage & Pace, 2009; Debarsy et al., 2015). Accordingly, future studies should move beyond binary classification accuracy alone and estimate marginal effects, spillover magnitudes, and regional policy scenarios.

Based on the series of analyses and discussions, this study draws several conclusions. First, parameter estimation in the spatial probit regression model was conducted using the Maximum Likelihood Estimator (MLE) method. Because the resulting equation is non-closed form and does not have an analytical solution, the Fisher Scoring numerical iteration procedure was applied to obtain accurate parameter estimates. Second, the significance of the spatial probit model was tested using the Maximum Likelihood Ratio Test (MLRT), while the estimation process under both population conditions and the spatial model framework continued to rely on the Fisher Scoring algorithm to address the complexity of the equation. The hypothesis testing results on poverty data in Sumatra in 2022 demonstrate that per capita expenditure, unemployment rate, economic growth as measured by GRDP, and expected years of schooling significantly influence the classification of poverty rates, both simultaneously and partially. Third, the spatial probit model results show a high level of predictive accuracy, with 110 regencies or cities correctly classified in accordance with the observed data, while 21 regions were classified differently. Overall, the spatial probit model with the Fisher Scoring approach achieved an accuracy rate of 83.97 percent, indicating reliable model performance. Finally, comparative evaluation reveals that Fisher Scoring-based spatial probit regression performs better than the RIS Simulator approach in terms of accuracy, sensitivity, and efficiency, as reflected in its lower Sum of Squared Errors (SSE). However, the RIS Simulator approach produced a higher specificity. Taken together, these findings indicate that the Fisher Scoring approach is more effective for generating precise poverty predictions, whereas the RIS Simulator offers a particular advantage in identifying the negative category in poverty data.