One of the fundamental properties of spatial data is the property of spatial dependence, which states that values that are close together in geographical space tend to be more similar than data values that are further apart. This property is known as Tobler’s law [23]. It can be expressed in a slightly different way when the observations are associated with an area (lattice data), such that an observation of a variable at location i carries some information about what is observed for the same variable in areas that are close to i [11]. The way to include the information of the neighbors is through a contiguity matrix, which states the neighbor structure of the sites where the variable of interest is observed.

We consider the final grades of the freshman students in each course as lattice data. To propose a multivariate model for the scores, some basic areal models are reviewed in the next section.

Let S ⊂ R 2 , and consider a finite partition of S given by the set { S 1 , … , S n }, i.e., S = ⋃ i = 1 n S i and S i ∩ S j = ∅ for all i ≠ j. The elements of the partition are referred to as areal units and are linked to a set of responses Y = ( Y 1 , … , Y n ) ⊤ .

The spatial variation in the response is modeled using covariates and a spatially structured random effect component Φ = Φ 1 , … , Φ n ⊤ , where each Φ i represents the spatial random effect for areal unit S i . A widely used approach to model these random effects is to define the value at S i conditionally on a linear combination of values from neighboring units.

Introduced by [4], CAR processes are commonly employed to capture this spatial dependence. A CAR model is specified via its full conditional distributions as: (2.1) Φ i ∣ Φ − ( i ) ∼ N ρ ∑ j ∈ ∂ i b i j Φ j , τ i − 1 , i = 1 , … , n , where j ∈ ∂ i denotes the indices of areal units that share a common boundary with S j , the b i j determine the strength of spatial interaction, ρ is a spatial autocorrelation parameter, τ i is the conditional precision.

By the Hammersley-Clifford theorem and Brook’s lemma, the set of full conditional distributions in (2.1) uniquely defines the joint distribution, Φ ∼ N n ( 0 , D τ ( I n − ρ B ) − 1 ) , where B = { b i j } i , j = 1 n ∈ R n × n satisfies b i i = 0, and D τ = diag { τ 1 , … , τ n }.

Typically, the neighborhood structure is encoded through a binary contiguity matrix W = { w i j } i , j = 1 n ∈ R n × n , where w i j = 1 if areal units S i and S j share a common boundary, and w i j = 0 otherwise. It is also conventional to set w i i = 0. The corresponding diagonal matrix of neighborhood counts is defined as D w = diag { w 1 + , … , w n + } ∈ R n × n where w i + = ∑ j ≠ i w i j represents the number of neighbors of S i .

These definitions are commonly used to specify the matrices B and D τ as follows: D τ = τ D w and B = D w − 1 W 1 , thereby providing the full specification of the CAR process.

[2] extended the univariate CAR model to the multivariate setting by considering Φ = ( Φ 1 , … , Φ p ) ∈ R n × p , where each Φ i = ( Φ 1 , i , … , Φ n , i ) ⊤ ∈ R n represents the spatial random effects for the i-th variable across all regions. Under suitable assumptions, the joint distribution of Φ can be expressed as vec ( Φ ) ∼ N n p ( 0 , Σ MCAR ) , where Σ MCAR = Λ − 1 ⊗ ( D w − ρ W ) − 1 . Following [12], Λ ∈ R p × p can be interpreted as a precision matrix that encodes the conditional association between variables within the same areal unit, while the matrix ( D w − ρ W ) − 1 captures the spatial association between different areal units for a given variable. The exact form of Λ − 1 is given in Web Appendix B.

This model is referred to as MCAR(ρ, Λ). If ρ is fixed at 1, the resulting model becomes an intrinsic multivariate CAR model, denoted IMCAR(Λ).

[26] laid the initial groundwork for extending the above model to a multivariate framework based on an adaptation of the linear model of coregionalization (LMC) for areal data. Subsequently, [5] introduced a unifying framework known as the M-model, which allows for distinct spatial autocorrelation parameters across the multivariate responses.

In the M-model, the multivariate spatial effect is represented as Φ = Θ M, where Θ = ( θ 1 , … , θ p ) ∈ R n × p is a matrix whose columns are mutually independent spatial processes θ i ∼ N n ( 0 , ( D w − ρ i W ) − 1 ), each following a CAR model with its own spatial autocorrelation parameter ρ i . The matrix M = { m i j } ∈ R p × p defines the loadings of the different underlying CAR spatial effects for each variable. Therefore, each spatial random vector effect can be written as: Φ j = ∑ i = 1 p θ i m i j , j = 1 , … , p .

The implementation of the IMCAR ( Λ ), MCAR ( ρ , Λ ), and M-models is commonly conducted under a Bayesian framework due to the explicit conditional structure [10]. The computational implementation can be constructed, e.g., using WinBUGS [2] or the R package INLAMSM [16], which allows the implementation of the models discussed in Section 2.2.

For lattice data, units in close proximity to each other with similar values exhibit a spatial pattern indicative of positive spatial autocorrelation. Identifying groups of units in close proximity with high values is often of particular interest because they suggest a cluster of elevated risk, possibly originating from a common source [19, p. 11].

In this work, we have chosen to use multivariate CAR models due to their compelling advantages for modeling multivariate areal data. These models jointly capture spatial dependence both within and across variables, and their sparse precision matrices allow for computational efficiency, making them particularly suitable for large datasets. Compared to alternatives such as multivariate Gaussian processes, MCAR models are more scalable while maintaining modeling flexibility. Moreover, they integrate seamlessly into hierarchical Bayesian frameworks, facilitating robust parameter estimation and uncertainty quantification in complex spatial settings.

The dataset employed in this study encompasses a compilation of academic records for all students enrolled at Universidad Técnica Federico Santa María from 2018 to 2022. We restricted our analysis to three years, as 2018 marked the implementation of a revised teaching methodology in first-year mathematics courses, rendering data from earlier years not directly comparable.

Between 2018 and 2022, several disruptions prevented these years from being treated as a continuous sequence. In 2018, a student strike halted teaching for nearly three months. In 2019, Chile experienced a social uprising that forced universities to switch to online instruction from October onward, with the fall semester extending into February 2020. During 2020–2021, lectures were delivered online due to the pandemic, and 2022 marked the return to regular in-person operations.

With multiple entries, each record in the dataset represents an individual instance of a student’s academic performance in a specific subject, along with their personal and academic details.

The dataset was preprocessed to streamline its structure for a focused spatial and academic performance analysis. This included data cleaning and the exclusion of incomplete cases. A logarithmic transformation was applied to the USC variables to stabilize their variance, as their values ranged between 600 and 800. Additionally, since the social variables were scaled within the ( 0 , 1 ) interval, the transformation helped integrate all variables into a common analytical framework.

To examine differences in outcomes before, during, and after the pandemic, the first semesters of the years 2018, 2020, and 2022 were selected across the Valparaíso Region (V) and the Santiago Metropolitan Region (RM). The choice of these regions was driven by their having the largest student populations as can be seen in Table 1, and hosting the principal campuses of the university. In addition, other regions exhibit sparse enrollment data, often with several counties registering zero observations. It is important to highlight that the Valparaíso Region comprises 36 counties, whereas the Santiago Metropolitan Region consists of 52 counties.

In the Valparaíso Region, there were 689 students in 2018, 616 in 2020, and 567 in 2022. Similarly, in the Santiago Metropolitan Region, there were 881 students in 2018, 1010 in 2020, and 898 in 2022. A detailed breakdown of first-year student counts in these two regions is available in Tables 6 and 7 of the Supplementary Material.

Temporal effects were not incorporated into the analysis for two main reasons. First, the available data spans only a limited number of academic years (2018, 2020, and 2022), which is insufficient to reliably capture long-term temporal patterns. Second, the student population changes from year to year, with each cohort bringing different characteristics and backgrounds. This inter-cohort variability introduces additional complexity that makes it difficult to isolate temporal effects in a consistent and interpretable way. Therefore, we chose to focus the analysis on spatial patterns within each year, where the assumptions and comparisons are more stable and meaningful.

The grades (recorded on a 0–100 scale) obtained by first-year students in the Mathematics I (MAT021), Introduction to Physics (FIS100), Programming (IWI131), and Chemistry and Society (QUI010) courses were averaged by county across the Valparaíso Region and Santiago Metropolitan Region (refer to Tables 2 and 3 in the Supplementary Material).

Each record now directly links academic performance in specific subjects with spatial data. Additionally, in pursuit of fitting an appropriate model, various economic, social, and educational variables that might explain the grades obtained in each county are considered in an auxiliary dataset (refer to Tables 4 and 5 in the Supplementary Material). These are as follows:

The observed patterns (see Figure 1) in grades are similar in Santiago and Valparaíso, with the first term of 2020 yielding the highest scores in all subjects. This highlights that during the first semester of the pandemic, the average scores were significantly greater than those in the first semester before and after the pandemic.

This hypothesis has been confirmed through multiple tests (details not reported here). We also observe that on both campuses, the best performance is in Introduction to Physics, while Programming and Mathematics I show the lowest performance. The passing threshold for a course is set at 55 out of 100. In terms of score stability over the years, Chemistry and Society is the strongest course.

In Figure 2, the maps displaying the final grades in each course for the Valparaíso and Santiago regions are presented. In cases where data are unavailable for a specific county, the respective areal unit is colored orange, and the missing value is imputed using the average of the other counties to avoid bias. Visually, it is evident that in 2020, the grades are higher than those in both the prepandemic and postpandemic periods. Various explanations have been proposed for this phenomenon, making it challenging to investigate from a scientific perspective. For instance, the limited experience of professors with online teaching methods may have led to the use of inadequate evaluation methods in this context. Additionally, the unintended collaboration among student groups during midterms and exams, which involves sharing crucial information, has been identified as a significant factor contributing to unprecedented approval rates. Despite the challenges in providing a rigorous explanation for this behavior, the data underscore an existing pattern that calls for further analysis.

According to Moran’s index (Anselin, 1995 [1]), significant autocorrelation exists for most of the scoring variables examined in the study (refer to Web Appendix A). Consequently, it is worthwhile to develop spatial models that account for spatial interactions both between and within variables.

To incorporate the covariates described in Section 3.1, we model the average student performance per county and year using a Bayesian hierarchical spatial model. Let Y i t denote the average score for county i = 1 , … , n at year t = 1 , … , p. In our study, the number of years is p = 3 and the number of counties is n = 36 for the Valparaiso Region and n = 52 for the Santiago Metropolitan Region.

We assume a Gaussian likelihood log ( Y i t ) ∼ N ( μ i t , τ − 1 ), with mean μ i t = α t + x i t ⊤ β t + ϕ i t , where α t is a year-specific intercept, x i t ∈ R q is the vector of q = 4 covariates (PI, MDI, USC, RPS), β t ∈ R q is the vector of regression coefficients that vary by year, and ϕ i t denotes spatially structured random effects modeling residual spatial variation.

The spatial random effects ϕ i t for all counties and years are modeled jointly as vec ( Φ ) ∼ N n p ( 0 , Σ ), where Σ is specified according to the spatial models described in Section 2.2

Inference is performed using INLA [18]. Further details on priors and posterior inference are provided in the Supplementary Material (Web Appendix C). In particular, Figure 1 and Table 1 in the Supplementary Material present forest plots with 95% credible intervals for each course, method, year, and area. These plots provide an overall perspective on the dispersion associated with the variables considered.

As a sensitivity analysis, we refitted the MCAR model under a likelihood scaled by the county–year sample size n i t . The fixed-effect posterior means differ by at most 0.06 on the log scale, and credible-interval widths differ by less than 2 %, indicating that the conclusions are robust to sample-size weighting (Web Appendix F of the Supplementary Material).

The full posterior summaries for all fitted models are reported in the Supplementary Material (Tables 8 and 9). According to the Deviance Information Criterion (DIC) [21] and the Widely Applicable Information Criterion (WAIC) [24, 25], where lower values indicate better fit (see also Tables 10 and 11), the MCAR model provides the best overall fit, supporting the estimation of the spatial autocorrelation parameter rather than fixing it. Consequently, Figure 3 shows the fitted maps for the four subjects and the three years using the MCAR specification only.

To assess model adequacy, we inspected standardized predictive residuals from the MCAR specification. Figure 4 presents QQ-plots based on residuals computed on the log-response scale for the four subjects in both regions. The plots indicate an overall reasonable agreement with normality, especially in the central part of the distribution, although some noticeable deviations are present in the tails.

Across the four covariates, only one association has a credible interval that excludes zero: in Santiago, the municipal development index is positively associated with grades in MAT021 in 2018, with β ˆ = 0.77 and a 95 % credible interval of [ 0.15 , 1.40 ]. On the original grade scale, this corresponds to roughly an 8 % higher county-average grade per 0.1 increase in the index. The remaining patterns are weaker. The admission score is positively associated with grades in IWI131 in Santiago across the three years, with point estimates near 0.8 but credible intervals that always include zero. The poverty index changes sign across courses and years, with large but uncertain magnitudes; the rate of private schools shows no stable pattern. Most estimates therefore carry substantial posterior uncertainty, and the associations above should be read as suggestive rather than conclusive.

The change in the sign of the PI coefficient between the lockdown semester (2020) and the post-pandemic semester (2022), particularly for IWI131 in Valparaíso, is consistent with the temporary buffering role of pandemic-related social policies documented across Latin America [22, 6]. With the data at hand we cannot test this mechanism formally, and we present this only as one possible reading of the observed pattern.

These findings suggest the value of admissions policies that adopt a more holistic approach, taking into account the disparities that students endure because of their socioeconomic background. Universities could explore initiatives to implement useful indicators based on socioeconomic factors for freshman students to mitigate the disadvantages faced by these students, aligning with broader educational equity objectives. These new policies could complement existing initiatives driven by the government of Chile, such as the PACE1¹

https://acceso.mineduc.cl/admision-universidades-2023/portal-pace/