The binary spatial autoregressive (BSAR) model is a class of spatial regression that explicitly accounts for dependence in a binary outcome by allowing each observation to be influenced by neighboring observations. When applied to networks, those spatial relationships are naturally represented as edges between nodes, so BSAR models extend directly to network data. One example of these spatial econometric models is credit scoring, where financial outcomes such as default or non-default can exhibit dependencies across related entities. Capturing these dependencies can improve predictive performance and provide a more realistic representation of risk propagation in financial systems. Although estimation for a single network feature is well understood, scaling these models to real-world datasets is challenging: feature selection, i.e. identifying which features should contribute to the spatial component and which should enter the standard covariate structure, must be performed in tandem with estimation. However, a large number of covariates leads to substantial computational burden. In this work, we employ a one-step procedure that simultaneously discovers a useful partition of features and estimates the BSAR parameters using evolutionary algorithms as a global optimizer. The proposed approach enables the construction of implicit network structures from real-world data and supports the development of more accurate and scalable binary BSAR models.