The Acquisition of Tonal Hierarchies in Western Music During School Years: A Re-Analysis of 40 Years of Research

Aims and Objectives

We conducted a critical review of four decades of research on the development of IRTH sensitivity. Despite a sizeable volume of studies, the findings remain heterogeneous and partially contradictory, particularly regarding developmental trajectories and the influence of musical training. Building on Gordon’s (2012) theory of music learning, which posits that core aspects of tonal sensitivity are largely developed by age 9, we formulated four exploratory research questions. By examining these questions, this study aimed to (1) assess the magnitude of age-related changes in IRTH sensitivity, (2) examine whether sensitivity increases significantly beyond age 9, (3) investigate the role of musical training, and (4) explore the effects of different operationalizations of IRTH measurement.

Research Questions

Research Question 1 (RQ1). How substantial is the average age-related development of IRTH sensitivity?

Research Question 2 (RQ2). Does IRTH sensitivity continue to increase beyond the age of 9, or does it level off, as predicted by Gordon’s learning theory?

Research Question 3 (RQ3). Is IRTH sensitivity significantly higher in musically trained individuals than in those without musical training?

Research Question 4 (RQ4). Does the implicit measurement of IRTH sensitivity via response time yield lower estimates than more explicit operationalizations?

Procedure

We conducted a systematic and comprehensive literature search (What Works Clearinghouse, 2020) between January and April 2022 to identify eligible studies. The study selection followed a predefined process outlined in the review protocol (see Supplementary Material S1), which was not published before the review was conducted. In the final sample, we included only primary research studies of the preliminary literature corpus that used a hypothesis-testing approach, followed either a cross-sectional or longitudinal design, were published between 1982¹ and April 2022, investigated a sample of healthy elementary and secondary school students, and defined the measurement of IRTH sensitivity as a dependent variable based on either a listening task (probe tone paradigm, response time measures, and goodness-of-fit ratings in syntax-violation paradigms) or a creative production task (harmonization of a given melody, tonal composition, and improvisation to a given chord sequence). These operationalizations were further specified as described below.

First, the probe tone paradigm (Krumhansl & Shepard, 1979) involves presenting one or two of the 12 chromatic tones (the "probe") after establishing the tonal context (e.g., an ascending major scale), followed by a brief silence. Participants rate how well the probe tones fit the tonal context. This process continues until all 12 chromatic tones are rated. The perceived stability of each tone is measured based on participants’ goodness-of-fit ratings. The correlation of these ratings with the prototypical tonal hierarchy (see above) defines participants’ sensitivity to IRTH. Although no standardized protocols or psychometric criteria² have been established, the probe tone paradigm has been consistently replicated across various contexts, demonstrating its robustness (Morgan et al., 2019; Sauvé et al., 2021).

Second, some studies use goodness-of-fit ratings in a syntax-violation paradigm (e.g., Corrigall et al., 2022). Participants judge how well the final tone of short melodies fits within the preceding harmonic context with varying degrees of tonal congruence (e.g., a tonic note such as "C" in C major versus a non-diatonic note such as "C#"). A greater difference in ratings between congruent and incongruent endings indicates IRTH acquisition, with higher ratings for congruent endings reflecting higher IRTH sensitivity.

Third, Schellenberg et al. (2005) and Corrigall et al. (2022) focus on response times (as introduced by Janata & Reisberg, 1988). Rather than explicitly evaluating the tonal congruence of melody endings, participants rate other musical features, such as the timbre of the final tone, while hearing a priming sequence with either a tonally congruent (expected) or incongruent (unexpected) final tone. Faster and more accurate responses are expected when the target tone or chord is more tonally congruent with the preceding context, indicating stronger IRTH sensitivity.

Fourth, other methods involve creative production tasks, such as composing (Wilson & Wales, 1995), improvising while listening to a pre-recorded tonal chord sequence (Paananen, 2003), or harmonizing a melody (Paananen, 2009). The outcome variables are defined by the degree of tonal fit and timing of the tones or chords used by the participants. Wilson and Wales (1995, p. 102) report a substantial to almost perfect inter-rater agreement (cf., Landis & Koch, 1977) for expert assessments of compositions $\left(.69 \le \kappa \le .92\right)$.

In the next step, three independent and trained coders (including the first author) coded the included studies according to a previously developed protocol (see Supplemental Material S2). After calculating the initial inter-rater reliability, which showed almost perfect agreement $\left(\kappa = .92\right)$, coding discrepancies were resolved by consensus. The studies vary in whether and how they report participants' formal musical training; therefore, we calculated the percentage of those with training for each age group. Sufficient data are available in 11 studies, while other variables are coded as not applicable (“n/a”).

Data Analysis

Data were analyzed using two approaches. First, we conducted a Bayesian three-level meta-analysis to provide a quantitative summary of the systematic review, offer an overview, identify notable studies, and investigate differences in IRTH measurements. Second, we conducted a model comparison analysis of the cross-sectional data of a single study (Krumhansl & Keil, 1982) investigating the growth trajectory of IRTH acquisition.

In the meta-analytical approach, effect sizes were estimated from primary studies using the metafor package (Viechtbauer, 2010) in R (R Core Team, 2022) and, if needed, transformed into Cohen’s $d$ following Borenstein and Hedges (2019, pp. 214–234). The procedure for effect size calculations, raw data, and a markdown script for the Bayesian three-level meta-analysis are shown in Supplemental Materials S3, S4, and S5.

In the subsequent analysis step, the effect sizes were aggregated and weighted using a Bayesian three-level meta-analysis with the R package brms (Bürkner, 2017), following the procedure of Harrer et al. (2021). Unlike a one-level fixed effect model that only attributes variance to sampling error, the two-level random effects model separates variance into sampling errors³ and between-study heterogeneity (${\mathrm{\tau }}_1^2$). This model captures inherent differences between studies, such as measurement methods.

Our model further includes a third level of variance $\left({\mathrm{\tau }}_2^2\right),\mathrm{ }$reflecting the nested structure of effect sizes within studies and accounting for the variance within single studies, such as different age groups. This multilevel approach (for a formal model description, see Supplemental Material S6) better reflects the nested structure of our data and avoids the statistical issues of traditional univariate meta-analysis, which can misestimate heterogeneity and increase the risk of false positives (Harrer et al., 2021; Hedges, 2019). We compared models with differing levels of complexity by computing marginal likelihoods using the bridge sampling method implemented in the brms package.

We chose Bayesian meta-analysis for two key reasons. First, it allows for the explicit modeling of heterogeneity uncertainty (${\tau }^2$), which can yield more stable estimates in small-sample contexts—especially when applying weakly or moderately informative priors (Harrer et al., 2021). Notably, although Bayesian methods are not inherently robust to small samples or poor study quality, they provide a principled framework for incorporating prior knowledge and quantifying uncertainty. Second, prior information can improve the posterior estimation of key parameters, such as $\mathrm{\mu }$ (intercept) and $\tau$ (standard deviations; Röver, 2020).

To address the context-dependence of choosing the prior parameter settings (Gelman et al., 2015), we initially used two separate and weakly informative zero-centered priors for $\mathrm{\mu }$ and$\mathrm{ }\mathrm{\tau }$, as recommended by Williams et al. (2018) and Harrer et al. (2021, Setting Prior Distributions, para. 5):

Subsequently, we performed sensitivity analyses to estimate the impact of prior parameter choices on the posterior density distributions of the parameter estimates. Therefore, following Röver (2020) and Turner and Higgins (2019), we used two alternative weakly informative priors for both $\mathrm{\mu }$ and$\mathrm{ }\mathrm{\tau }$, while τ was modeled as a standard deviation parameter constrained to be positive, implying half-distributions for the respective priors:

We did not conduct a meta-regression using age as a moderator because of the limited number of studies, heterogeneity in age reporting, and insufficient information on age distribution in several samples.

We examined the magnitude of age-related development in IRTH sensitivity (RQ1) by fitting a Bayesian three-level meta-analytic model using the brms package in R. Small effects ($d \approx 0.2$) are often considered the minimum threshold for practical relevance in developmental research (e.g., Gignac & Szodorai, 2016), and the average difference in IRTH sensitivity between older and younger participants meets or exceeds the benchmark for a small effect. Insufficient information reported in the primary studies prevented us from conducting a meta-regression with musical training as a moderator (RQ2). Therefore, the impact of musical training on IRTH development remains an open research question.

To address RQ3, we estimated the impact of operationalization on effect size using a Bayesian three-level meta-regression in brms. The model included operationalization as a moderator, accounted for the nesting of effect sizes within studies, and incorporated their standard errors. Random effects were specified at the study and within-study levels.

We examined RQ4 using Krumhansl and Keil (1982), who report notable within-study heterogeneity, to compare learning-theory-informed models using data from. The dataset comprised four age groups assessed across five probe tone conditions. Based on age-specific difference scores, we fitted linear (linear, quadratic, cubic) and non-linear (sigmoidal, logistic, mixed, saturated) functions to model the development of tone judgment stability. Analyses were conducted using the R packages lme4 (Bates et al., 2015) and nlme (Pinheiro & Bates, 2024; see S7 for data and S8 for code).

Bayesian Three-Level Meta-Analysis

We initially estimated a four-level meta-analysis to aggregate effect sizes, as two articles (Matsunaga et al., 2020; Schellenberg et al., 2005) reported independent studies, suggesting potential within-article heterogeneity. However, although a comparison between the initial four-level model and a three-level meta-analysis model showed that the Bayes factor, BF₁₀ = 1.92, provided insufficient evidence to decisively favor either model (for a general discussion, see van Doorn et al., 2021), the three-level model was ultimately chosen because of its slightly better data fit and based on the principle of parsimony (Vandekerckhove et al., 2015). The model comparison did not support further simplifying the model to a random-effects structure⁴.

Before interpreting the results, we confirmed model convergence for our three-level structure by conducting posterior predictive checks and examining the $\widehat{R}$-values of the parameter estimates (Harrer et al., 2021). All $\widehat{R}$-values $\left(\widehat{R}\le 1.01\right)$ indicated successful convergence (Bürkner, 2017). Additionally, visual inspection confirmed that the posterior distributions aligned with the initial unimodal normal distribution, supporting our assumption of normality (see Appendix Figure A2-1). Sensitivity analyses using different weakly informed priors further support the robustness assumption of our results (Table 3).

We assessed publication bias using a Bayesian Egger test, which indicated a significant intercept (β₁ = 7.46, 95% CI [3.12, 11.67]), suggesting potential publication bias. However, previous studies have recommended conducting this approach on a corpus of at least 30 studies because asymmetry tests for standardized mean differences are prone to inflated Type I error rates (Pustejovsky & Rodgers, 2019; Renkewitz & Keiner, 2019). A visual inspection of the funnel plot (see Appendix Figure A2-2), showing the standard error of the reported effect sizes as a function of their magnitude, suggests the potential for slight asymmetry; however, this is primarily driven by Krumhansl and Keil (1982; Study ID 1). Thus, we found no significant evidence of general publication bias caused by selective reporting.

We further explored potential sources of bias by examining the model-based relationships between μ and τ₁ (Figure 2a), and μ and τ₂ (Figure 2b). Although visual inspection suggested only marginal associations, Pearson correlations revealed statistically significant but negligible effects in both cases: μ and τ₁, r(11998) = -.09, p < .001, 95% CI (-.11, -.07); μ and τ₂, r(11998) = .08, p < .001, 95% CI (.06, .10). These results indicated that any potential relationship between effect sizes and heterogeneity is minimal, thus providing little evidence of systematic publication bias.

Click to enlarge

Figure 2

Heatmap of the Joint Posterior Distribution of Model Parameters μ (x-axis) and Standard Deviations τ (y-axis)

Note. Figure 2a depicts the relationship between μ and τ₁ (between-study standard deviation). Figure 2b shows the relationship between μ and τ₂ (within-study standard deviation). The color gradient reflects the density scale, with higher values (lighter regions) indicating higher density and lower values (darker regions) indicating lower density. Red lines in each plot represent the maximum likelihood point estimates of both parameters.

Based on the observed data, hierarchical model, and prior choices, we produced a point estimation for the pooled medium-sized effect (Mdn_d = 0.57). As indicated by the credible intervals of the pooled effect size (Table 3), we could further conclude that the general so-called “true” effect size—expressed as Cohen’s d—reflects a medium to nearly large difference in IRTH sensitivity between younger and older participants. Moreover, the 95% credible interval for the effect lies entirely above zero (0.37, 0.77)⁵ along with decisive evidence against a point-null hypothesis (BF₁₀ = 11999). We further assess whether the age-related increase in IRTH sensitivity exceeds a negligible magnitude, that is, whether it is neither trivially different from zero nor within the range of a small effect $\left(H_0:0 \le \left|d\right| < 0.2, {BF}_{01} > 10\right)$. Therefore, we used Bayes factors to quantify the evidence that the result exceeded the small-effect threshold $\left(H_1: \left|d\right| \ge 0.2, {BF}_{10} > 10\right)$ by using a Bayesian model comparison approach. Our analyses revealed a Bayes Factor providing decisive support for the alternative assumption (BF₁₀ = 1332.33) indicating that the age-related increase in IRTH sensitivity is unlikely to be marginal. Instead, it should be best interpreted as a substantively meaningful skill-development effect with a high probability of at least medium-to-strong magnitude. We further illustrated the cumulative posterior evidence by plotting the empirical cumulative distribution functions for $\mathrm{\mu }$, ${\mathrm{\tau }}_1$ and ${\mathrm{\tau }}_2$ (see Appendix Figure A2-3) to facilitate the interpretation of the probability that these parameters exceed meaningful thresholds.

Furthermore, based on the hierarchical model and its estimated parameters, we revealed a model with parameter estimates as solution in which at least one effect size estimate from each study falls within the 95% credible interval of the pooled effect size (Figure 3). At first glance, Schwarzer et al. (1993) appear to be an exception; however, although the point estimates of effect sizes in their study are predicted to fall well outside the 95% credible interval of the pooled effect size, their respective 95% credible intervals still intersect with this range.

Click to enlarge

Figure 3

Forest Plot of the Bayesian Three-Level Model With Estimated Effect Sizes of Individual Studies and the Pooled Effect Size

Note. The densities represent their respective posterior distributions. Medians serve as point estimates for the effect sizes. Blue shading highlights estimated effect sizes and parts of the posterior distributions that fall within the 95% credible interval of the pooled effect size.

Although the vast majority of the effect sizes reported in investigated studies fall within the expected 50–80% credible intervals according to our hierarchical model (Figure 4), some observed effect sizes reported by Krumhansl and Keil (1982) deviate from this pattern. Specifically, the results of two studies (ID 11 and ID 12) exhibit exceptionally high effect sizes. Although our model considers these values possible, they are highly unlikely from a statistical standpoint. For example, according to our model, the probability that the estimate for the effect size for ID 12 in Krumhansl and Keil (1982), denoted as ${\mathrm{\theta }}_{\mathrm{1,12}}$ in the posterior samples, is greater than or equal to the observed effect size $({\text{ES}}_{\mathrm{1,12}}=3.02)$ is relatively low at 2.1%, $P\left({\mathrm{\theta }}_{\mathrm{1,12}}\ge {\text{ES}}_{\mathrm{1,12}}|y,M\right)=0.021$. In contrast, the observed effect size for ID 11 has an even lower probability of occurring, at less than 1% $(P\left({\theta }_{\mathrm{1,11}}\ge {\text{ES}}_{\mathrm{1,11}}|y,M\right)=0.005)$.

Click to enlarge

Figure 4

Model-Based Probability Predictions for the Occurrence of Each Effect Size

Note. Whereas the vast majority of effect sizes have a high probability of occurrence based on the model-based expectations, some values have extremely low statistical probability, such as ID 11 and ID 12.

The multilevel model revealed significant heterogeneity, divided into between- and within-study heterogeneity. Although the median between-study variance, ${\mathrm{\tau }}_1^2=0.04,\mathrm{ }95\mathrm{\%}\mathrm{ }CrI\mathrm{ }\left(0.00,\mathrm{ }0.16\right),$ accounted for only 18% of the total variance, the majority of the heterogeneity (82%) could be attributed to within-study variance, ${\mathrm{\tau }}_2^2=0.18,\mathrm{ }95\mathrm{\%}\mathrm{ }CrI\mathrm{ }\left(0.07,\mathrm{ }0.36\right)$.

However, differences in age-group compositions across studies prevented us from answering RQ2 using this analytical approach, which is discussed below. Furthermore, RQ3 could only be addressed in an exploratory manner, as the limited number of studies and inconsistent reporting of relevant sample characteristics hindered formal statistical testing. Table 2 summarizes the proportions of musically trained and untrained individuals across the included samples.

We addressed RQ4 by examining whether differences in IRTH measurement type could explain the variability in effect sizes. The central goal of a meta-analysis is to account for heterogeneity in the data (Borenstein, 2019). Therefore, we followed Corrigall et al. (2022) and introduced a categorical moderator representing six measurement types: probe tone, response time, goodness-of-fit rating, composition, improvisation, and harmonization. The analysis revealed largely overlapping 95% credible intervals and only small differences in posterior means, suggesting that the measurement type had no systematic effect (Figure 5). However, comparing model fit using R² values (Hayes, 2022) showed a slight improvement in the moderated model, R² = .51, 95% CrI (.30, .72), compared to the model without a moderator, R² = .48, 95% CrI (.27, .70). Although the Bayes factor of BF₁₀ = 38.80 indicated a strong statistical preference for the moderated model (van Doorn et al., 2021), the practical significance of this improvement remained uncertain, given the modest differences in parameter estimates. Therefore, we directly addressed RQ4 by comparing the response time measures with all other operationalizations. These comparisons yielded weak-to-moderate evidence against a systematic difference for goodness-of-fit ratings (BF₁₀ = 0.48), harmonization (BF₁₀ = 1.00), and improvisation (BF₁₀ = 0.84) in lower estimates for response time measures compared to probe tone ratings ${(BF}_{10 }= 12.99)$. Overall, the results indicated smaller effect sizes for response time measurements in probe tone ratings only, with no consistent reduction in effect sizes for response time measurements.

Click to enlarge

Figure 5

Posterior Distributions of the Intercept and Moderator Levels of the Bayesian Three-Level Moderated Meta-Analysis

Note. Effect sizes are standardized mean differences (d). Thin error bars are 50% CrI of d. Thick error bars indicate 95% CrI of d.

Taken together, these results suggest that response time measurements do not consistently yield lower IRTH sensitivity estimates compared to more explicit operationalizations. Therefore, the current evidence does not support RQ4, except for a notable difference in probe tone ratings.

The observed within- and between-study heterogeneity was likely caused by specific artifacts, such as sample and characteristics. However, we cannot pinpoint these factors more precisely because of the limited information in the primary studies. However, certain possibilities were excluded based on the model comparison analysis presented below.

Thus far, we have reported evidence of a medium-to-strong age-related increase in IRTH sensitivity, although the high variance somewhat obscures this trend. However, these analyses do not clarify whether age-specific learning gains are present, particularly (1) whether acquisition continues beyond age 9 and (2) whether there are sensitive periods during school years that would be reflected in a non-linear learning trajectory. Thus, we conducted further analyses to address these uncertainties.

Model Comparison Analysis

Based on the AIC, BIC, and RMSEA model comparison criteria, the non-linear mixed model best fits Krumhansl and Keil’s (1982) data (see Appendix Table A1-1), suggesting no linear or sigmoidal increase in IRTH sensitivity between ages 6 and 20. Figure 6 shows the age-related development of various IRTH sensitivity outcomes according to the non-linear mixed model. Based on the maximum increase in IRTH sensitivity at age 20 (k = 0.08), we concluded that IRTH continues to develop into adulthood and is not fully developed by age 9.

Click to enlarge

Figure 6

Development of IRTH Modelled as a Non-Linear Function of Age

Note. Standardized mean differences (d) as indicator of tonal recognition stability based on reanalyzed cross-sectional data of Krumhansl and Keil (1982, pp. 247–248) as measured with the two sample tones probe tone paradigm (see section procedure for explanation).

Furthermore, the distinct trajectories shown in Figure 6 suggest that outcomes may strongly depend on specific task characteristics, even within a single-measurement paradigm. Although the differentiation between pairs of diatonic and non-diatonic tones, as well as the tonic triad versus other diatonic tones, shows a steep increase, more subtle IRTHs, such as the preference for the tonic over other tonic triad tones, show little to no increase, even in adults.