3.1. Data selection

This study utilised data from the Organisation for Economic Co-operation and Development (OECD) PISA 2018 cycle to address the research questions. PISA assesses the academic knowledge and skills of 15-year-old students across various domains, including reading, mathematics, and science, while also collecting detailed background information through student, parent, teacher, and school-level questionnaires.

The 2018 dataset, which includes data from over 600,000 students across 79 countries and economies, is particularly valuable as it represents the first and, at the time of writing, only instance where intercultural communication skills, including FL learning, were included in the assessment (OECD, 2019a). This allowed for an in-depth analysis of these skills alongside student performance in mathematics.

To ensure alignment with the research questions and maintain data validity, a subsample of students was selected based on language-use criteria and data completeness. Our study focuses on language learners in mainstream programmes where FLs are taught as subjects, rather than through immersion. To minimise the likelihood of including proficient bilinguals – who, according to prior research, may exhibit different cognitive profiles (White & Greenfield, Reference White and Greenfield2017) – only students who reported using the same language at home and at school were included in the sample. For example, if a student’s home language (e.g., Spanish) differed from the language of assessment (e.g., English), their data were excluded. Additionally, cases with missing or invalid responses for key variables, as well as extreme outliers, were removed to meet model requirements and reduce bias in the results (see Supplementary Material Appendices S1 and S2 for details).

For the first research question – examining the relationship between time spent on school-based FL learning programmes and mathematics performance – we used data from the student, parent, and school questionnaires. After applying our filters, 34% of cases in the PISA dataset were removed due to invalid or incomplete data, 13% were excluded based on the language-mismatch filter, and a further 4% were excluded due to the outlier removal process. Consequently, the final sample for the first research question comprised 300,656 students from 73 countries and economies.

The second research question focused on the role of specific FL programme characteristics in predicting mathematics performance. This analysis incorporated all datasets in the first study, supplemented with data from teacher questionnaires. To ensure conceptual relevance, the dataset was filtered to include only language teachers, as the teacher-level variables incorporated in the models – such as years of experience, multicultural teaching practices, and cultural competence – are specific to FL instruction and would not apply to teachers of other subjects. Because only 19 countries administered the teacher questionnaire and the proportion of language teachers was relatively small, the final sample for the second question was substantially smaller, consisting of 53,459 students and 16,342 teachers, representing 31% and 23% of their respective total samples. Further details on the data selection and preparation processes are provided in the Supplementary Material (Appendices S1 and S2).

3.2. Variable description

3.2.1. Foreign language learning variables

Our primary predictor variable was FL learning time in minutes per week (FLMINS). FLMINS serves as a proxy for exposure to FL learning, as direct measures of proficiency were unavailable.

This variable was calculated by combining two PISA 2018 variables: FL periods per week (ST059Q04HA) and average minutes per class period (ST061Q01NA), both derived from students’ input in the student/parent questionnaire.

The calculation was performed as follows:

$$ FLMINS=\left( ST059Q04 HA\right)\times \left( ST061Q01 NA\right). $$

To ensure reliability, outliers were addressed using two approaches. First, FLMINS values were capped at the total minutes of instruction (TMINS), an OECD-derived index reflecting realistic instructional time. Second, within each school, the interquartile range (IQR) method was applied to exclude extreme values exceeding the upper bound (James et al., Reference James, Witten, Hastie and Tibshirani2023). This ensured that only FLMINS values disproportionately high relative to other values within the same school were removed. Tables 1 and 2 provide a summary of all the variables used, including their description, sources, and values.

Table 1. Description of variables: linear mixed model variables (Study 1)

^* Composite variables are calculated as described in the Methods section and Supplementary Material.

Table 2. Description of variables: random forest variables (Study 2)

Note: All variables from Study 1 are also included in Study 2.

^* Composite variables are calculated as described in the Methods section and Supplementary Material.

3.2.4. Teaching programme variables

The following variables were included to ensure that the effects of FL learning are assessed independently of the broader benefits of engagement with multicultural teaching and programmes.

Multicultural curriculum (MCCUR). This index was derived from five PISA school questionnaire items (SC167Q01HA–SC167Q05HA) asking whether the formal curriculum for 15-year-olds includes components such as intercultural communication, cultural knowledge, openness to intercultural experiences, and FL instruction. Each item was coded 1 = Yes, 0 = No, and equally weighted. The score represents the proportion of these elements present at the school level, yielding a range from 0 to 100:

$$ \mathrm{MCCUR}=\frac{\mathrm{Number}\ \mathrm{of}``\mathrm{Yes}"\mathrm{responses}}{5}\times 100. $$

This approach follows PISA’s aggregation conventions (OECD, 2020) and reflects the extent to which multicultural content is formally integrated into the curriculum. Moreover, all components are derived from question SC167 and designed to be interpreted together, ensuring consistency in meaning. The inclusion of these elements aligns with research emphasising the role of cultural education in supporting language learning and cognitive development (Dlaska, Reference Dlaska2000).

Multicultural teaching practices (MCTEACH). This composite index captures the extent to which schools implement practices that actively promote cross-cultural interaction and understanding, such as hosting teachers from other countries, teaching the history of international cultural groups, offering exchange programmes, celebrating cultural festivities, encouraging intercultural communication, and promoting FL skills (Cok, Reference Cok2021; Savignon & Sysoyev, Reference Savignon and Sysoyev2005).

The index was calculated from seven PISA items (SC159Q01HA, SC165Q02HA, SC165Q06HA–SC165Q09HA, TC207Q05HA). Each item was coded 1 = Yes, 0 = No, and equally weighted. The composite score represents the proportion of these practices implemented at the school level, yielding a range from 0 to 100:

$$ \mathrm{MCTEACH}=\frac{\mathrm{Number}\ \mathrm{of}``\mathrm{Yes}"\mathrm{responses}}{7}\times 100. $$

Internal consistency was high (Cronbach’s α = 0.923), supporting the reliability of this measure.

Teacher’s self-efficacy in multicultural environments (GCSELF). It is an OECD-calculated index reflecting teachers’ confidence in addressing multicultural challenges, based on self-assessment across various dimensions such as raising awareness of cultural differences and managing multicultural classrooms effectively.

4.1. Method: Linear mixed modeling

To investigate the association between FL learning time (FLMINS) and students’ mathematics performance, a linear mixed-effects model (LMM) was fitted using the lme4 package in R (Bates et al., Reference Bates, Mächler, Bolker and Walker2015). This model accounts for the nested structure of the data, where students are grouped within countries, and it controls for sources of variability at both levels.

The model was specified as follows:

$$ {\displaystyle \begin{array}{l} Math\ Score\sim FLMINS+ MMINS+ GENDER+ GRADE\\ {}\hskip9em +\hskip2px ESCS+ MISCED+ FISCED+ SCHLOC\\ {}\hskip9em + SCHLTYPE+\hskip2px ADMISSION\\ {}\hskip9em +\left(1|\mathrm{CNT}\mathrm{STUID}\right)+\left(1|\mathrm{CNT}\right).\end{array}} $$

This specification estimates the association between FLMINS and mathematical performance (Math Score), while controlling for time spent on mathematics (MMINS), demographic factors (gender, grade), socioeconomic and cultural background (ESCS), parental education level (MISCED, FISCED), and school characteristics (location, type, and admission process). Random intercepts at the student (CNTSTUID) and country (CNT) levels account for individual differences and country-level variability, respectively. This structure addresses unobserved heterogeneity within countries and among individual students.

An initial screening of FLMINS revealed the presence of outliers, even after applying data cleaning filters described earlier. To address this, FLMINS was winsorised at 1,500 minutes per week, the maximum value for mathematics instruction time (MMINS). This threshold aligns with educational practice, where FL instruction typically receives fewer periods compared to core subjects like mathematics (OECD, 2023a). Winsorisation, a statistical technique that limits extreme values in a dataset by replacing outliers with a specified value, was employed to preserve meaningful data points – such as those from schools prioritising FL study – while preventing FLMINS to exceed realistic bounds and having disproportionate influence on analyses (Gelman & Hill, Reference Gelman and Hill2006; Tabachnick & Fidell, Reference Tabachnick and Fidell2018; Wilcox, Reference Wilcox2017).

In addition, we tested for a potential nonlinear relationship between FLMINS and Math Score using a generalised additive model (GAM) with the mgcv package for R (Wood, Reference Wood2017). The analysis indicated a linear relationship between the variables of interest, supporting the use of an LMM (for details, see Supplementary Material Appendix S2). Given the size and complexity of the data, the use of a GAM, instead of an LMM, would have introduced substantial computational challenges without significantly improving model performance.

Importantly, school size (SCHSIZE) and class size (CLSIZE) were initially considered as predictors but were ultimately excluded, as their inclusion would have removed all countries for which the OECD did not report or compute values for these variables. Because complete data were required to estimate the LMM, including SCHSIZE and CLSIZE would have substantially reduced the sample without offering commensurate analytical benefit. To include the effect of class numbers without significantly reducing the sample, student–teacher ratio (STRATIO) was used instead. Other variables were initially considered but later removed, such as student age, which was replaced by GRADE, and highest parental income, which was replaced by parental education level (MISCED, FISCED).

4.2. Study 1 results

Given the potential associations among predictors, Pearson’s correlations (Olivoto & Lúcio, Reference Olivoto and Lúcio2020) and variance inflation factors (VIFs) (Lüdecke et al., Reference Lüdecke, Ben-Shachar, Patil, Waggoner and Makowski2021) were calculated to assess multicollinearity. Although some variables exhibited moderate-to-strong correlations (r = .70 for ESCS/MISCED, r = .67 for ESCS/FISCED), no significant collinearity concerns were identified (VIFs: 2.43 for ESCS, 1.78 for FISCED, and 1.89 for MISCED) (for details, see Supplementary Material Appendix S2).

The linear model, which included both fixed and random effects, achieved a conditional R ² of .865, while the marginal R ² was .162. This indicates that the inclusion of random intercepts for students and countries considerably improved the model fit, with 86.5% of the variance in mathematics scores explained by the full model compared to 16.2% achieved by fixed effects alone. Table 3 presents a summary of the results obtained.

Table 3. Linear mixed-effects model results: fixed effects estimates and random effects intercepts

All fixed effects were statistically significant (p < .001), which was expected given the large sample size, suggesting a robust relationship between the predictors and mathematics scores. FLMINS showed a positive association with maths scores, estimated at 4.8 PISA points higher per additional hour of FL instruction (or b = .08×60 minutes). For context, across the middle 50% (110 to 240 min/week), this effect corresponds to a positive variation of 8.8–19.2 points in predicted achievement.

Note that some of the variables included in the analysis were categorical – specifically GENDER, GRADE, MISCED, FISCED, SCHLOC, SCHLTYPE, and ADMISSION – and were treated accordingly. However, instead of presenting the individual effects for each category level, we only report their combined effect.

Among the demographic variables, SES (ESCS) and GRADE had the strongest effects. A one-standard-deviation increase in ESCS corresponded to a 29.33-point difference in mathematical scores – for the middle 50% of the sample (−0.86 to +0.66), this represents a 44.5-point spread. Similarly, each additional grade level was associated with an estimated difference of 28.08 points (e.g., 56.16 points between grades 8 and 10).

Other school-level factors showed smaller but significant associations. Schools located in larger cities (SCHLOC = 5) reported scores 25.24 points higher than those in rural areas (SCHLOC = 1), averaging a 6.31-point difference in scores per level. Likewise, students attending private independent schools (SCHLTYPE = 1) showed scores 10.40 points higher than those attending public schools (SCHLTYPE = 3) and about 5.20 points higher than those in government-dependent private schools. Similarly, academic selectivity (ADMISSION) accounted for a 4.84-point difference between schools that “always” versus “never” use placement tests (b = 2.42).

The remaining variables showed lesser effects on mathematics scores. Results indicate that the time spent learning mathematics (MMINS) had a limited positive association, corresponding to about 1.8 points per hour, a significantly smaller effect than the time spent learning FL. Gender differences were similarly modest, with boys scoring 10.43 points higher than girls on average. Surprisingly, higher maternal education levels (MISCED) was associated with an average difference of −2.16 point, while each increase in paternal education level (FISCED) corresponded to an estimated difference of −0.83 point in mathematics scores.

To explore differences among countries, these were grouped by economic status according to the latest United Nations World Economic Situation and Prospects report (United Nations, 2025). This approach was selected since national economic growth has been shown to mediate the effect of SES – the strongest predictor in our main model – on academic achievement (Kim et al., Reference Kim, Cho and Kim2019). Research also indicates there is a reciprocal relationship between economic growth and academic performance, further emphasising the rationale for this stratification (Hanushek & Woessmann, Reference Hanushek and Woessmann2021; Valero, Reference Valero2021; Woessmann, Reference Woessmann2016).

This model revealed a consistent positive association between increased FL instruction time and mathematics performance across all economic categories (Figure 1). However, the magnitude of this effect varied notably by economic status. Economies in transition, countries such as Azerbaijan, Belarus, and Georgia, exhibited the steepest slope, with mathematics scores on average 11.4 points higher per additional hour of FL learning per week. In contrast, developing countries (e.g., Argentina, Morocco, and Vietnam) recorded the lowest baseline scores and the smallest overall estimated difference, with an average of 4.2 points per hour of FL instruction. Developed countries (e.g., Canada, Japan, and Norway) also displayed a positive trend, with an average difference of 7.2 points per hour of FL learning (for details, see Supplementary Material Appendix S2).

Figure 1. Linear mixed effects model results: Effect of FLMINS on mathematics scores according to countries’ economic status.

4.3. Robustness tests

To assess the stability of the observed association between foreign-language learning and mathematics achievement, we conducted additional analyses using alternative subsamples and model specifications.

4.4. Study 1 discussion

Study 1 examined the relationship between foreign language instruction time (FLMINS) and mathematics achievement using a linear mixed model (LMM). The main model revealed a positive association: Mathematics scores were approximately 4.8 points higher per additional weekly hour of foreign-language (FL) instruction. This finding reinforces prior research on the academic benefits of FL learning (Nucette et al., Reference Nucette, Hamamura, Leitao and Biedermann2024; Woll & Wei, Reference Woll and Wei2019), as well as previous research on the indirect link between FL learning and improved numeracy performance, mediated by EF (Al-khresheh & Karmi, Reference Al-khresheh and Karmi2024; Shoghi & Ghonsooly, Reference Shoghi and Ghonsooly2018).

Although the model explained a significant portion of the variance in mathematical scores (R ² = 86.5%), the analysis highlighted that most of the variability is attributable to country and individual differences, particularly socioeconomic level and grade. School-level characteristics (e.g., school location, type), while influential, played a comparatively smaller role in the variation.

Robustness checks further strengthened these findings against methodological challenges. The association persisted, though at a reduced magnitude (3 PISA points per FL-hour/week), when the analysis was restricted to multilingual education systems, indicating the effect is not confined to language learning contexts. Expanding the sample to include language-mismatch cases confirmed the overall stability of the association and indicated that it is not driven by language-background selection or test-language disadvantage. Crucially, decomposing exposure into within- and between-school components (Mundlak adjustment) revealed that the effect holds true among peers in the same school (b = 0.06) and across schools based on their FLMINS allocation (b = 0.09). These patterns suggest that the observed relationship is not only attributable to fixed differences between schools or non-random placement of students. While they do not establish causality, they are consistent with our theoretical interpretation that FL learning may support mathematics achievement indirectly through cognitive processes.

Stratification by country-level economic status demonstrated that FLMINS was consistently associated with higher mathematical performance across all economic categories. The strongest effects were observed for economies in transition, an outcome that could be attributable to specific pedagogical mechanisms, similar educational reforms, or shared attitudes toward FL learning. Another possible explanation is cultural in nature, with intrinsic motivation and resilience playing a role (Huisman et al., Reference Huisman, Smolentseva and Froumin2018). As neither of these mechanisms was the primary focus of the present study, additional research is needed to explore these hypotheses. In contrast, developing and developed economies exhibited smaller gains. In developing economies, this may reflect structural constraints such as limited access to quality education or inadequate teacher training, whereas in developed economies, well-established teaching practices and curricula may lead to smaller effects. These findings align with previous literature highlighting the influence of national economic contexts on educational outcomes (Hanushek & Woessmann, Reference Hanushek, Woessmann, Rosén, Hansen and Wolff2017).

While the positive association between FLMINS and mathematics achievement was robust across specifications, it should be interpreted as correlational rather than causal, given the cross-sectional nature of the data. We acknowledge that unobserved school- and country-level factors – such as specific teaching practices, student–teacher rapport, curriculum design, and educational policies – likely contribute to the observed relationship and require further exploration.

To address the limitations of traditional approaches, a machine learning (ML) framework was adopted for the second research question. ML techniques offer enhanced capacity to manage complex datasets, tolerate noise, and identify influential predictors with high accuracy – making them particularly suitable for analysing large-scale educational datasets such as PISA. This approach also enables cross-validation of findings from Study 1, contributing to the robustness and generalisability of the results.

5.1. Method: Machine learning approach (random forest)

To evaluate the impact of FL programme characteristics on mathematical achievement, an ML approach was adopted, specifically the random forest (RF) method, implemented via the R package randomForest (Breiman, Reference Breiman2001; Liaw & Wiener, Reference Liaw and Wiener2002).

The primary rationale for selecting this method was its high predictive accuracy. ML methods, particularly tree-based models, are designed to construct statistical models that generalise effectively to unseen data, whereas traditional models often perform optimally on the data to which they were fitted. In other words, ML models tend to produce predictions that closely align with actual outcomes (Hastie et al., Reference Hastie, Tibshirani and Friedman2009). Due to their predictive power, such algorithms have been widely applied across various fields, including psychology and education, to uncover patterns that might otherwise remain undetected (Balabied & Eid, Reference Balabied and Eid2023; Liew et al., Reference Liew, Hamamura and Uchida2025; Nachouki et al., Reference Nachouki, Mohamed, Mehdi and Abou Naaj2023; Yarkoni & Westfall, Reference Yarkoni and Westfall2017).

RF algorithms, along with other ML approaches, have been increasingly applied to identify predictors of mathematics and science literacy, using large-scale international student assessments, such as PISA and Trends in International Mathematics and Science Study (TIMSS) (Arroyo Resino et al., Reference Arroyo Resino, Constante-Amores, Gil Madrona and Carrillo López2024; Gil-Madrona et al., Reference Gil-Madrona, Guerrero-Muguerza, Infantes-Paniagua and Martínez-López2025; Song & Cutumisu, Reference Song and Cutumisu2025). Their growing use likely reflects their suitability for large and complex datasets. These algorithms also demonstrate a greater resilience to missing data and are less influenced by outliers, ensuring robust and reliable results (Hastie et al., Reference Hastie, Tibshirani and Friedman2009; Liaw & Wiener, Reference Liaw and Wiener2002). RF models also provide an important interpretative advantage by quantifying the relative importance of predictor variables, offering insights into their practical significance on the outcome (James et al., Reference James, Witten, Hastie and Tibshirani2023; Varian, Reference Varian2014; Yarkoni & Westfall, Reference Yarkoni and Westfall2017). These features make them superior to more traditional statistical models, including LMMs, which have known issues to detect complex interactions and may overestimate the significance of certain variables in large samples (Lin et al., Reference Lin, Lucas and Shmueli2013; Thiese et al., Reference Thiese, Ronna and Ott2016).

Building on these strengths, the present study applies an RF model to examine a broader range of predictors than those considered in Study 1 by incorporating factors related to FL teaching programmes, such as multicultural curriculum and teaching practices. The primary objective was to identify the variables that most robustly and accurately predict the outcome variable – Math Score. Accordingly, the model included established demographic predictors of academic achievement, alongside variables specific to FL programmes (see Table 2).

The predictors can be categorised as follows:

• • Student-level predictors: Gender, grade level, SES (ESCS), parental involvement (PARINV), and parents’ highest educational levels (MISCED, FISCED).

• • Foreign language programme characteristics: Time allocated to FL study (FLMINS), teacher experience in years (TEACHEXP), teacher self-efficacy in multicultural environments (GCSELF), and indicators of multicultural teaching practices and curriculum (MCTEACH, MCCUR).

• • School-level predictors: Student–teacher ratio (STRATIO), school size (SCHSIZE), school ownership (SCHTYPE), school location (SCHLOC), and mathematics instructional time (MMINS) as a control for time-on-task effects.

The dataset was randomly split into training (80%) and testing (20%) subsets. This approach ensured that model fit indices reflected the model’s ability to predict outcomes on unseen data (Varian, Reference Varian2014), thereby providing external validation of model performance (hold-out testing). Internal model performance was evaluated using out-of-bag (OOB) validation, which provides a cross-validation estimate based on bootstrap resampling during training. Performance was quantified using root-mean-square error (RMSE), mean absolute error (MAE), and the coefficient of determination (R ²).

5.2. Study 2 results

The RF model performed well across both internal and external validation methods. The OOB cross-validation yielded an R ² of 0.962, an RMSE of 17.29, and an MAE of 9.43. Performance on the independent hold-out test set was nearly identical (R ² = 0.966, RMSE = 17.16, MAE = 9.41), indicating minimal overfitting and strong generalisability within the PISA dataset. Given the range of mathematics scores (136–799.6) and the mean of 485.1 in the data, the RMSE and MAE correspond to deviations of approximately 3.5% and 1.9% of the mean score, respectively.

These results indicate that the model explained approximately 96% of the variance in the OECD-generated PVs for mathematics scores, with a low prediction error, as illustrated in Figure 2.

Figure 2. Random forest model results: Predicted vs actual values.

Predictor importance was assessed using the percentage increase in mean-squared error (%IncMSE), a metric that quantifies the reduction in the model’s predictive performance when the values of a given variable are permuted while keeping all other variables unchanged. A higher %IncMSE value indicates greater predictor importance of a particular variable, as altering its values results in an error of higher magnitude (Breiman, Reference Breiman2001; Liaw & Wiener, Reference Liaw and Wiener2002).

Using this approach, Economic, Social, and Cultural Status (ESCS) emerged as the most influential predictor (%IncMSE = 364.9%). This means that randomly changing ESCS values, while keeping all other variables constant, introduces an error 3.6 times greater than that of the original model. The second most important predictor, closely following ESCS, was FLMINS (363.3%), further confirming its strong association with mathematics performance.

Other important predictors included the implementation of multicultural curricula (MCCUR = 277.9%) and maths instruction time (MMINS = 274.9%), while SCHSIZE (259.4%) and STRATIO (228.0%) demonstrated moderate importance. Notably, teachers’ multicultural self-efficacy (GCSELF) had a %IncMSE of 0%, indicating that it had no predictive power for mathematics scores in this model. Figure 3 presents a ranked summary of predictor importance.

Figure 3. Random forest model results: Variable importance based on the percentage increase in mean-square error (%IncMSE – mean decrease in accuracy).

5.3. Robustness tests

5.4. Study 2 discussion

Study 2 employed an RF model to examine associations between FL programme characteristics and mathematics achievement. The model explained a large share of variance (R ² = 0.97, RMSE = 3.5%, MAE = 1.9% of mean score), with OOB performance closely matching test-set results, indicating strong internal consistency and generalisability. However, the high level of model fit requires careful interpretation. PISA PVs are conditioned on background variables (such as SES, grade, and gender) that are also included as predictors – creating an outcome-predictor dependency known to inflate fit indices (Zieger et al., Reference Zieger, Jerrim, Anders and Shure2022). The large sample size and the presence of several highly predictive variables may have further contributed to this effect.

A robustness check using only PV3 produced a modestly lower R ², consistent with the expected influence of pooling results across PVs, which smooths random error and can slightly increase fit. Crucially, when key PV-conditioning variables (ESCS, GRADE, and GENDER) were removed from the model, FLMINS emerged as the most important predictor across both importance metrics, and the ranking of the remaining variables remained highly stable. This pattern indicates that although absolute accuracy metrics are influenced by PV construction, the relative importance structure – particularly the prominence of FLMINS – reflects a genuine and robust signal rather than a modeling artifact.

Taken together, these findings suggest that the RF results should be interpreted primarily as exploratory, highlighting variable importance and nonlinear relationships rather than providing definitive estimates of predictive accuracy. Nonetheless, the stability of the importance rankings offers meaningful insights into the role of FL programme features within the broader ecosystem of predictors of mathematics achievement.

Consistent with Study 1, time dedicated to FL instruction (FLMINS) emerged as a highly influential predictor – second only to SES (ESCS) in the full model – a pattern that held across specifications. The RF model’s capacity to capture nonlinearities and interactions further enriches this interpretation, offering a more nuanced understanding of how FL-related factors operate alongside other student- and school-level characteristics (Bates et al., Reference Bates, Mächler, Bolker and Walker2015; Breiman, Reference Breiman2001).

Among the FL programme features, multicultural curriculum integration (MCCUR) demonstrated the strongest association with mathematics scores. This effect is likely mediated by the cognitive benefits of robust FL learning and cultural exposure, which is strongly linked to enhanced EFs – including working memory, response inhibition, and cognitive flexibility – all critical predictors of mathematical reasoning (Cragg et al., Reference Cragg, Keeble, Richardson, Roome and Gilmore2017; ten Braak et al., Reference ten Braak, Lenes, Purpura, Schmitt and Størksen2022). Additionally, multicultural education has been associated with enhanced critical thinking and problem-solving skills (Aslan & Aybek, Reference Aslan and Aybek2020; Qondias et al., Reference Qondias, Lasmawan, Dantes and Arnyana2022), which also contribute to successful mathematical performance.

While this study offers credible evidence that features of FL programmes, particularly instructional time (FLMINS) and multicultural curriculum (MCCUR), are meaningfully associated with mathematical achievement, these findings do not establish causality. As in Study 1, additional research is required to disentangle the mechanisms underlying these associations.

6.1. Dual pathways and contextual moderation

The robustness of this positive association provides theoretical nuance, highlighting overlapping student-level and school-level components. The within–between decomposition showed a significant within-school effect: even when comparing peers in the same school, students with greater FL exposure tend to score higher in mathematics – consistent with our cognitive transfer theory. The stronger between-school effect suggests that unobserved institutional features may covary with FL time (e.g., curriculum, instructional quality, selection/admissions, timetabling).

These findings also engage with potential trade-offs in highly multilingual systems (Dentella et al., Reference Dentella, Masullo and Leivada2024). In a robustness check restricted to Luxembourg, Switzerland, Canada, Finland, Spain, and Singapore, the association persisted but was modest (3 points per additional FL hour/week). This pattern suggests that the link is sensitive to system-level features (curricular alignment, language-of-instruction policy), and learner-level factors (FL acquisition, usage), which may interact or influence cognitive development (Schiltz et al., Reference Schiltz, Lachelin, Hilger and Marinova2024). Together, these considerations provide a framework through which our findings can coexist with prior reports of reduced maths achievement in contexts of high language burden, rather than contradict them.

6.2. Programme features and equity

Consistent with prior research, SES (ESCS) emerged as the strongest determinant of mathematics achievement (b = 29.33), highlighting persistent educational disparities linked to economic resources (Kim et al., Reference Kim, Cho and Kim2019; OECD, 2023b; Rodríguez-Hernández et al., Reference Rodríguez-Hernández, Cascallar and Kyndt2020). However, our findings suggest that FL study may serve as a potential avenue for promoting equity. Indeed, the magnitude of these associations is considerable: the predicted score difference associated with two hours of weekly FL instruction (9.6 points) is roughly equivalent to the observed gender gap, while a six-hour instructional volume (28.8 points) parallels the impact of a full SD increase in SES. Moreover, the positive association between FLMINS and maths achievement was observed across all economic status groupings, showing the strongest effect size in economies in transition. This suggests a compensatory role for FL programmes, where the strongest positive associations are observed in settings with fewer existing structural advantages.

The ML approach (Study 2) complemented this finding by confirming FLMINS as a strong predictive factor, comparable in importance to ESCS in predicting mathematical performance. Furthermore, multicultural curriculum integration (MCCUR) emerged as the most significant feature among the evaluated FL programme components, supporting the view that FL benefits may accrue not just from time on task, but from instruction delivered authentically, promoting higher-order thinking and culturally responsive learning environments.

6.3. Other determinants of achievement

Beyond the primary determinants, the models identified several conventional and anomalous predictors of mathematics achievement. Consistent with expectations, grade level exhibited a strong effect (b = 28.08 per year of schooling), as students at higher grade levels have been exposed to more advanced topics and have had more opportunities to practice, apply concepts, and master fundamental skills.

Regarding institutional factors, school location (SCHLOC) and school type (SCHLTYPE) demonstrated moderate effects, indicating that students in urban settings and private schools achieve higher scores – a finding often linked to differential access to specialised resources and infrastructure. As expected, admission type (ADMISSION) also showed a significant positive effect, showing that selective or merit-based enrolment policies contribute to higher mean achievement. Institutional size factors, such as school size (SCHSIZE) and student-teacher ratio (STRATIO), were found to be of moderate predictive importance in Study 2, confirming that while these structural attributes may shape the learning environment, their impact on achievement is less pronounced than curricular content (Mohammadpour & Abdul Ghafar, Reference Mohammadpour and Abdul Ghafar2014).

Our results also showed that boys outperformed girls by an average of 10.43 points, consistent with the documented gender gaps in mathematics performance (Anaya et al., Reference Anaya, Stafford and Zamarro2022), which are likely influenced by contextual factors unrelated to biology, such as SES and cultural expectations (Breda et al., Reference Breda, Jouini and Napp2018; Johnson et al., Reference Johnson, Burgoyne, Mix, Young and Levine2022).

An interesting anomaly emerged regarding parental education (MISCED, FISCED). Contrary to typical findings (Davis-Kean et al., Reference Davis-Kean, Tighe and Waters2021), both maternal and paternal education showed a small negative association with scores. This finding is likely explained by the high correlation between parental education and the already-included ESCS variable (r ≈ 0.70), suggesting that SES (ESCS) accounts for the primary variance. The residual, small negative effect might reflect subtle mechanisms, such as parental time constraints associated with highly demanding careers, or could simply be the result of multicollinearity within the structural model.

Finally, instructional time allocated to mathematics (MMINS) and teacher quality variables exhibited nuanced roles. The effect of increasing MMINS was nearly three times smaller than the association observed for FLMINS (1.8 versus 4.8 PISA points per hr/week), suggesting that gains from instructional time alone are incremental and likely secondary to other factors, such as pedagogical quality or student engagement. Correspondingly, teacher self-efficacy in multicultural environments (GCSELF), a predictor in Study 2, showed no predictive power for mathematics achievement. Contrary to previous studies (Abacioglu et al., Reference Abacioglu, Volman and Fischer2020; Nuenay et al., Reference Nuenay, Cariga, Bualan and Banes2024), this non-finding suggests a potential intention-action gap, where teacher confidence in managing diverse classrooms primarily impacts noncognitive outcomes (e.g., classroom climate) rather than translating directly into measurable gains in mathematics proficiency (Feng et al., Reference Feng, Zhang, Yang, Lin and Maulana2024).

Collectively, these studies offer significant theoretical insights and suggest promising directions for policy. Our findings provide empirical support for the hypothesis that FL learning is a robust correlate of mathematical achievement. This relationship, established through a dual-pathway framework of cognitive transfer and institutional quality signaling, aligns with recent syntheses highlighting the potential of FL programmes to foster mathematical performance through relatively short instructional periods (Nucette et al., Reference Nucette, Hamamura, Leitao and Biedermann2024).

Given the cross-sectional nature of the PISA data, these results are hypothesis generating for policy. They suggest that FL education should be viewed as a complementary component of a well-rounded curriculum, rather than a competitor for STEM hours. Notably, the 2.6-to-1 efficiency ratio – where FL instruction exhibits a significantly stronger association with maths scores than additional mathematics time – suggests that at the margin, prioritising cognitive breadth may be more effective for developing mathematical literacy than increasing subject-specific volume.

Furthermore, the consistent predictive role of FLMINS indicates it may serve as a lever for reducing academic disparities. Strategic FL investment offers a hypothesised “cognitive lift” capable of narrowing achievement gaps associated with gender and SES. Specifically, policy should prioritise integrating multicultural curricula to foster the abstract reasoning essential for numeracy and ensuring equitable access to high-quality FL programmes in under-resourced schools. However, it is crucial that any implementation of such initiatives consider the system-level constraints and potential trade-offs noted in multilingual contexts.