Modeling Junior High School Participation Rates in East Nusa Tenggara, Indonesia, Using K-Nearest Neighbor-Based Spatial Regression

Education is universally acknowledged as a cornerstone of human civilization and social progress. It not only equips individuals with cognitive and technical skills but also fosters moral reasoning, empathy, and critical awareness—elements that sustain a nation’s social fabric. Within the framework of the Sustainable Development Goals (SDG 4), quality education is recognized as both a human right and a transformative instrument for equity and empowerment [1, 2]. Across the globe, education reforms seek to reduce disparities and broaden participation among all social groups, emphasizing that knowledge accessibility should transcend geographic and economic boundaries [3-5]. However, despite these global commitments, educational participation remains unevenly distributed, particularly in developing countries, where spatial inequality continues to reflect socioeconomic disparities.

Indonesia, as a vast archipelagic country, faces a complex challenge in achieving equal access to education. The geographic fragmentation of islands, compounded by differences in infrastructure and economic development, often produces educational isolation in remote areas. This spatial divide is particularly evident in the eastern provinces, where geographical remoteness amplifies economic limitations. East Nusa Tenggara (NTT) stands as a vivid example of this paradox: a region rich in cultural and natural diversity but constrained by unequal educational access. According to Statistics Indonesia [6], participation at the junior high school level has improved—Gross Participation Rate (APS) rose from 86% in 2023 to 88.66% in 2024, and the Net Participation Rate (APM) increased modestly from 56% to 58.15%. However, nearly 130,000 school-age children remain outside formal education, reflecting deep-seated structural inequities that hinder inclusive development [7, 8].

This persistent challenge underscores the urgency to not only measure participation rates but to understand why and where inequalities persist. Education, in this sense, is not merely an institutional process but a spatial phenomenon—its distribution is shaped by topography, distance, and access. In the case of NTT, understanding how geography interacts with social and economic realities is essential to formulating responsive, evidence-based education policies that truly reach the periphery.

Education is not only a social institution but also a spatial process—its access, equity, and quality are deeply influenced by geography. The spatial distribution of schools, infrastructure, and human resources determines who benefits from education and who remains excluded. The school participation rate therefore reflects more than socioeconomic disparity; where physical distance, regional connectivity, and geographic isolation interact with social variables to shape educational outcomes. This understanding aligns with Tobler’s First Law of Geography, which posits that “everything is related to everything else, but near things are more related than distant things.” In this study, spatial theory is extended into the educational domain to explain how participation rates in one district may influence or be influenced by neighboring districts—an effect conventional non-spatial models cannot capture [9-14].

Spatial analysis is used in this research because educational phenomena, particularly in geographically fragmented regions like NTT, exhibit spatial dependence and spatial heterogeneity—conditions where observations in one area are statistically linked to nearby areas. Traditional regression models such as Ordinary Least Squares (OLS) assume independence among observations, making them unsuitable for data with spatial structure. When such assumptions are violated, OLS models often produce biased estimates, inflated significance levels, and misleading conclusions. In contrast, spatial regression explicitly incorporates the spatial relationships among regions through a spatial weight matrix, allowing for the identification of localized patterns, diffusion effects, and inter-district influences that are essential to understanding real-world education disparities.

Among various spatial weighting methods, this study adopts the k-nearest neighbor (KNN) approach because it is uniquely suited to the archipelagic and non-contiguous geography of NTT. Unlike Queen or Rook contiguity matrices, which define neighbors based solely on shared borders, the KNN method identifies neighboring regions based on geographical proximity rather than direct adjacency. This is crucial in contexts where districts are separated by sea but remain connected through social, economic, or educational interactions. The KNN model thus ensures that every district has a defined number of spatial relationships—preventing the exclusion of island districts that lack contiguous boundaries. This makes KNN particularly appropriate for spatial modeling in archipelagic provinces, where traditional contiguity-based models would underestimate spatial influence and distort policy interpretation [15-18].

Methodologically, the KNN spatial regression model reflects a relational approach to education research: it does not view each district as an isolated unit, but as part of a network of spatial interactions. This approach captures the notion that improvements in one district—such as better infrastructure investment—may create spillover benefits for neighboring districts. Hence, the KNN-based spatial framework does more than describe statistical correlation; it reveals the geography of interdependence that defines educational opportunity in multi-island regions. In this sense, spatial analysis is not merely a technical choice but a philosophical one: it recognizes that place and proximity are central to understanding inequality and to designing education policies that are equitable, targeted, and contextually relevant.

Despite the growing body of literature on education inequality in Indonesia, few studies have systematically examined the spatial dependence of school participation rates, particularly within island provinces like NTT where geography itself becomes an educational barrier. Most existing works focus narrowly on socioeconomic determinants—poverty, income, or employment—without acknowledging the geographic interactions that shape how one district’s development can influence another’s. This oversight limits understanding of the spatial diffusion and clustering of educational participation, resulting in policy interventions that often generalize across contexts and overlook regional diversity.

In addressing this gap, this study advances a novel perspective by integrating spatial econometric modeling into the analysis of educational participation, using the KNN weighting approach. This framework allows for a more nuanced exploration of education inequality in non-contiguous island regions, where districts may be physically separated but remain socially and economically interconnected.

The originality of this study does not derive from the routine application of KNN weighting to educational data—an approach that has appeared in various regional analyses. Its distinctiveness lies in questioning the spatial assumptions that underpin such applications when transferred to territorially fragmented systems.

Prevailing spatial education studies implicitly operate within geographically contiguous settings, where distance is treated as a stable and uniform metric of interaction. Such assumptions obscure the structural discontinuities that define archipelagic provinces like East Nusa Tenggara. Here, proximity is not merely geometric adjacency; it is conditioned by maritime separation, uneven transport infrastructures, and historically entrenched peripheralization. Spatial relations are therefore asymmetrical and mediated, rather than continuous and homogeneous.

Against this backdrop, the study does not treat KNN weighting as a technical default. Instead, it interrogates how neighborhood specification behaves under non-contiguous geography and socio-spatial fragmentation. By calibrating the k-parameter through the joint evaluation of information criteria and spatial autocorrelation diagnostics, the analysis reframes weighting matrix construction as an empirical question rather than a procedural step.

Substantively, junior high school participation is examined as a relational outcome shaped by inter-district spillover mechanisms. The findings indicate that educational disparities in peripheral regions are not confined within administrative borders but propagate through spatial interdependence. This challenges policy frameworks that conceptualize inequality as territorially isolated.

Through this analytical repositioning, the study demonstrates that spatial weighting choices carry substantive implications in discontinuous territories—an issue that remains insufficiently scrutinized in contemporary spatial education modeling.

Specifically, the study pursues three objectives:

(1) To identify spatial patterns and dependence in junior high school participation across NTT’s districts;

(2) To evaluate the extent to which socioeconomic variables—poverty, employment, and teacher availability—affect participation levels; and

(3) To determine the most suitable spatial model for representing the inter-district relationships that underpin regional education disparities.

The study’s contributions are twofold. Empirically, it enriches understanding of spatial inequality in education through the lens of regional connectivity, revealing how geographic structures shape learning access. Conceptually, it bridges spatial analysis and educational policy, providing a foundation for evidence-based, geographically targeted interventions. By highlighting the interplay of place and policy, this research contributes to the broader goal of achieving equitable, inclusive, and contextually grounded education for all—particularly in geographically fragmented regions such as NTT.

3.1 Multiple regression linear

The multiple linear regression model is a statistical model used to examine the relationship between more than one independent variable and one dependent variable. The equation of the multiple linear regression model can be defined as follows [42].

$Y_i=\beta_0+\beta_1 X_{i 1}+\beta_2 X_{i 2}+\ldots+\beta_{p-1} X_{i, p-1}+\varepsilon_i$      (1)

where,

$Y_i$: the value of the dependent variable in observation i

$\beta_0, \beta_1, \ldots, \beta_{p-1}$: Model parameter regresi linear berganda

$X_{i 1}, X_{i 2}, \ldots, X_{i, p-1}$: the value of the independent variable to $\mathrm{p}-1$ in observation i

$\varepsilon_i$: error in observation i

3.2 Spatial regression model

Spatial data inherently embodies locational influence, where proximity among observations gives rise to spatial dependence, while regional diversity produces spatial heterogeneity. These twin properties underpin the logic of spatial regression, which extends conventional econometric thinking by explicitly modeling spatial interaction and variation. As established in the canon of spatial econometrics [43], spatial effects are not statistical noise but integral structures that shape social and educational outcomes. The evolution of this discipline—driven by innovations in spatial data science and open-source computation—has enabled more transparent and reproducible analyses of complex regional pattern [44]. Recent developments, including the formalization of endogenous spatial regimes and advances in model specification search, refine how spatial heterogeneity and dependence are simultaneously identified and estimated [45-48]. Such methodological precision strengthens the analytical foundation of spatial studies in education, allowing researchers to uncover not only where inequalities occur, but how spatial structure itself sustains or mitigates those disparities. The location effect consists of two types, namely spatial dependence and spatial heterogeneity. The general model of spatial regression can be written as follows [49].

$\begin{gathered}y=\rho W y+X \beta+u \\ u=\lambda W u+\varepsilon, \varepsilon \sim N\left(0, \sigma_{\varepsilon}^2 I_n\right)\end{gathered}$    (2)

With:

$y$: Vector Response Variable Size $n \times 1$

$\beta$: vector regression parameter size $p \times 1$

$\rho$: parameter spatial coefficient lag dependent variable

$X$: a matrix of predictor variables measured $n \times(p+1)$

$W$: weighting matrix with size $n \times n$

$\lambda$: the spatial coefficient parameter on the error

$u$: an error vector that has a spatial effect with a size $n \times 1$

$\varepsilon$: error vector with size $n \times 1$

$n$: number of observations or Location

The general equation of the spatial regression model in Eq. (2) can be formed by other models as follows:

• If and called Spatial Autoregressive Model (SAR) with equations $\rho \neq 0 \lambda=0$

$y=\rho W y+X \beta+\varepsilon$      (3)

• If and called a SEM with the equation $\rho=0 \lambda \neq 0$

$y=X \beta+u$   (4)

$u=\lambda W u+\varepsilon$  (5)

• If and called SARMA with the equation $\rho \neq 0 \lambda=0$

$y=\rho W y+X \beta+u$     (6)

$u=\lambda W u+\varepsilon$    (7)

3.3 Spatial weighting matrix

A spatial weighting matrix is a square matrix of size N × N that expresses the effect of spatial dependency between locations. The configuration of spatial units can be represented by a matrix W, where the proximity (distance) between locations is used to construct the matrix.

$W=\left[\begin{array}{ccc}W_{11} & \cdots & W_{N 1} \\ \vdots & \ddots & \vdots \\ W_{1 N} & \cdots & W_{N N}\end{array}\right]$    (8)

Each element of the W matrix can be defined as follows:

$w_{i j}=\left\{\begin{array}{c}1, j i k a j \in N(i) \\ 0, \text { lainnya }\end{array}\right.$

$N(i)$ is a neighborhood group from location j. The matrix concept defined above is based on only two closest neighbors. The most common weighting matrix used when two locations are not neighbors to each other is the concept of distance. The distance in question consists of geographical, economic, and social distance.

3.4 Nearest neighbor k-approach

According to Tobler's law, all things are interrelated, but the closer they are, the stronger the relationship. Therefore, a spatial weighting matrix can be formed with a nearest neighbor approach as an alternative to form a weighting matrix for locations that are not neighboring each other. Previous research conducted by Jaya et al. [50] stated that the formation of the W matrix with the KNN approach has the goal of determining the most optimal number of neighbors with the function of choosing the Moran's I value that will provide the most significant spatial effect and minimizing the value of Akaike's Information Criterion (AIC). In the W matrix of the KNN spatial weighter, each row i has a K conjunct of column j with element 1 and other columns of value 0. In the KNN spatial weighter, the K-value of the nearest location is determined by the researcher. According to Jaya and Andriyana the KNN spatial weighting formula is defined as follows:

$W=\left\{\begin{array}{c}1, \text { central point } \mathrm{k} \text { is close to central point } \mathrm{i} \\ 0, \quad \text { others }\end{array}\right.$

Spatial dependency testing or spatial autocorrelation using Moran's I statistics is as follows:

$I=\frac{n \sum_{i=1}^n \sum_{j=1}^n W_{i j}\left(Z_i-\bar{Z}\right)\left(Z_j-\bar{Z}\right)}{\left(\sum_{i=1}^n \sum_{j=1}^n W_{i j}\right) \sum_{i=1}^n\left(Z_j-\bar{Z}\right)^2}$      (9)

With:

$i: 1,2, \ldots, n$

$j: 1,2, \ldots, n ; i \neq j$

$Z_i$: The value of the variable at location to i

$Z_j$: Variable value at location to j

$\bar{Z}$: the average of the value of the variable

$W_{i j}$: the weights used to compare locations to i and j

$n$: Sample size

$I$: I from Coefasia Global Moran

The value of the Moran index ranging from -1 to 1 indicates that there is a large autocorrelation between residuals at one location and another. Autocorrelation does not occur if the value of the Moran index is equal to zero. Hypothesis testing of Moran index parameters as follows:

$H_0: I=0$ (No spatial autocorrelation)

$H_1: I \neq 0$ (There is spatial autocorrelation)

$Z(I)$ is the statistical value of the Moran index, is the expected value of the Moran index, is the value of the variance of the Moran index. This test rejects when $E(I) \operatorname{Var}(I) H_0|Z(I)|>Z_{\alpha / 2}$. Meanwhile, to see the relationship between standardized observation values and standardized neighbor averages, you can use Moran Scatterplot. This can be used to identify spatial equilibrium or influences where the type of spatial relationship can be seen in Figure 1.

Figure 1. Plot scatter Moran

There are four quadrants that show the four types of spatial relationships between an area and other adjacent areas as follows.

The Moran scatter plot groups observations into four patterns based on how values relate to their neighbors, showing areas where high values tend to be near high values, low near high, low near low, and high near low.

Furthermore, the global spatial autocorrelation or Moran index is not able to accommodate the need to identify the relationship between observation sites and other observation sites. Therefore, it is necessary to provide information related to the indication of spatial relationships in each region with the Local Indicator of Spatial Autocorrelation (LISA) which can be written for each region i as follows:

$L_i=\frac{\left(x_i-\bar{x}\right)}{m_2} \sum_{j=1}^n w_{i j}\left(x_j-\bar{x}\right)$     (10)

With:

$m_2=\sum_{j=1}^n \frac{\left(x_j-\bar{x}\right)^2}{n}$       (11)

$L_i$ is LISA at Location i, n is the number of Observation Locations, is the observation value at location i, is the observation value at location j, is the average of the observation value, is the weighting matrix element between location i and j. Hypothesis testing of the parameters can be carried out as follows: $x_i x_i \bar{x} w_{i j} L_i$

$H_0: L_i=0$ (There is no spatial autocorrelation at the ith location)

$H_1: L_i \neq 0$ (There is spatial autocorrelation at the ith

Test statistics:

$Z\left(L_i\right)=\frac{L_i-E\left(L_i\right)}{\sqrt{\operatorname{Var}\left(L_i\right)}}$    (12)

$Z\left(L_i\right)$ is the statistical value of the LISA test of the i th location, is the expected value of the LISA of the ith location, is the value of the variance of the LISA of the ith location. This test rejects when $E\left(L_i\right) \operatorname{Var}\left(L_i\right) H_0\left|Z\left(L_i\right)\right|>Z \alpha / 2$.

The selection of k = 2 in the KNN weighting matrix was not arbitrary but grounded in both theoretical reasoning and empirical validation. A sensitivity analysis was conducted by testing k values ranging from 1 to 5, and the resulting models were compared using the AIC and Moran’s I statistic to evaluate model fit and spatial dependence, respectively. The analysis revealed that k = 2 produced the lowest AIC value, indicating the best balance between model complexity and explanatory power, and the most stable spatial autocorrelation structure, as reflected by consistently significant and positive Moran’s I value. This outcome aligns with established spatial econometric practices that recommend optimizing k to minimize AIC while maintaining statistically meaningful spatial interaction [50]. Thus, k = 2 represents the optimal trade-off between capturing local spatial dependence and preventing excessive smoothing that could obscure meaningful regional variations.

From a contextual standpoint, the geographical characteristics of NTT further justify the adoption of the KNN weighting scheme. As an archipelagic province composed of multiple non-contiguous islands, many districts are physically separated by sea boundaries, making conventional contiguity-based matrices (such as Queen or Rook) less appropriate, since they rely on shared borders to define spatial neighbors. The KNN approach, by contrast, defines spatial proximity based on Euclidean distance rather than border adjacency, ensuring that each district maintains a valid set of neighboring relationships even when geographically isolated. This methodological choice allows for a more realistic representation of social and economic interconnections that extend across island boundaries. Consequently, the KNN weighting matrix is not only statistically robust but also geographically consistent with the empirical nature of the study area, enhancing both the validity and interpretability of the spatial regression results.

3.5 Research location

The dependent variable in this study, Gross Enrollment Rate (GER) for junior high school, measures the total number of students enrolled—regardless of age—relative to the official school-age population (13–15 years). This indicator provides an overview of the inclusiveness and accessibility of the education system and is widely used by the Badan Pusat Statistik (BPS) and UNESCO Institute for Statistics as a global benchmark for assessing education equity [7, 51]. A higher GER reflects broader access to education but may also indicate over-age or under-age enrollment, which is common in developing regions with uneven education access [52].

Among the independent variables, poverty level (X₁) represents the percentage of the population living below the poverty line. Poverty remains a decisive socioeconomic determinant of educational participation, as financial hardship restricts families’ capacity to invest in schooling and often forces children to withdraw prematurely from formal education [53]. This relationship is strongly supported in the spatial and development literature, which underscores the interdependence between economic deprivation and educational outcomes. For instance, poverty has been shown to exhibit clear spatial patterns, where regions of entrenched deprivation reinforce educational disadvantage through localized feedback mechanisms [54]. Likewise, education is recognized as a critical driver of poverty reduction and socioeconomic sustainability, indicating a bidirectional relationship between these two dimensions [55]. Broader spatial analyses further reveal that poverty alleviation and social development are geographically contingent processes influenced by structural disparities across regions [56]. Consequently, incorporating poverty as an explanatory variable within a spatial regression framework is theoretically warranted, as it captures the embedded spatial dynamics linking economic inequality to educational participation across territorial contexts.

Conversely, GDP per capita (X₂) captures local economic capacity and household welfare, reflecting how regional income disparities influence investment in education (OECD, 2020). In addition, the Special Allocation Fund (DAK) (X₃) represents fiscal transfers from the national government aimed at supporting local infrastructure and education service improvements. Previous studies highlight that fiscal decentralization, such as DAK, plays a vital role in addressing education inequality across regions [57, 58].

This research was conducted in the province of East Nusa Tenggara, Indonesia. The data used in this study is secondary data obtained from the BPS-Statistics Indonesia, related to the participation rate of junior high schools along with the variables that are suspected to affect it in 22 districts/cities in NTT.

The variables used along with the operational definition can be seen in Table 1.

Labor and institutional indicators further explain the socioeconomic dynamics that shape education outcomes. The Labor Force Participation Rate (LFPR) (X₄) measures the share of the working-age population participating in the labor market. A high LFPR may indicate economic pressure that reduces children’s school attendance, especially in poorer households [19]. Lastly, the number of teachers (X₅) captures the availability of teaching personnel relative to the population, reflecting the human resource capacity of the education system. Teacher availability is consistently linked to improved access, learning quality, and school participation, particularly in remote and low-come regions [21, 59].

Table 1. Variables and operational definitions