There are many fundamentally different approaches of defining and classifying climate zones and biomes. To reduce ambiguity in our analyses, we defined climate zones and biomes based on the current literature. We also identified and defined disturbance regimes within these well-defined climate zones and biomes.

For the definition of broad climate zones, we followed the FAO classification used for global forest resources assessment and reporting (FAO, 2012). This classification identifies five major zones, namely tropical, subtropical, temperate, boreal and polar zones (FAO, 2012). The tropical zone is within the area bounded by latitudes of 23.5° north and 23.5° south. The subtropical zone is applied to the two belts between the tropics (± 23.5°) and 35° north and south of the equator. Temperate zones fall in the latitudes of 35–50° north or south, and have well-defined seasons with a distinct winter. The boreal zone falls between 50 and 60° north. When testing our first hypothesis, we grouped the data according to the above classification based on the geographic coordinates of the study sites in the databases.

The biome is an important construct for organizing knowledge about terrestrial ecosystems, for examining diversity-productivity relationships, modeling historical distributions and shifts following climate change. We followed the IUCN Global Ecosystem Typology (Keith et al., 2022), which recognizes seven terrestrial biomes: (1) tropical-subtropical forests (T1); (2) temperate-boreal forests and woodlands (T2); (3) shrublands and shrubby woodlands (T3); (4) savannas and grasslands (T4); (5) deserts and semi-deserts (T5); (6) polar-alpine biomes (T6); and (7) intensive land-use (T7) biome. Each of these biomes is characterized by different ecosystem functional groups described in detail in Keith et al. (2022). The present analysis was limited to T1, T2, T3, and T4 as these have trees as the main or co-dominating components. Deserts and semi-deserts and polar-alpine biomes were deemed outside the scope of this analysis. Although some samples may have come from T7, it was difficult to assign them with confidence in the existing databases.

A disturbance regime is defined as the combination of disturbance agents and disturbance attributes that characterize a particular landscape or region (Burton et al., 2020). Disturbance interactions are an important part of the disturbance regime (Burton et al., 2020; Sturtevant and Fortin, 2021). The main disturbance agents in forest ecosystems consist of abiotic (e.g., fire, drought, wind, snow, and ice) and biotic (e.g., insects and pathogens) agents (Fischer et al., 2013; Stephens et al., 2013; Sturtevant and Fortin, 2021). According to a recent systematic review of disturbance interactions studies (Sturtevant and Fortin, 2021), the most frequently investigated natural disturbance agent was fire accounting for 65% of studies, followed by wind (38%), insects (37%), water imbalance (drought or flooding; 15%), mammalian browsing/grazing (10%), and mass movement including erosion, debris flows, landslides and avalanche (7%). Accordingly, the focus of this analysis was on fire, wind and insects as disturbance agents.

Fire is an important natural disturbance factor in many boreal forests, in savannas and high mountain dry-land ecosystems (Fischer et al., 2013). Low-severity fire regimes are typified by frequent low-intensity fires where surface fuels are charred or consumed while damage to overstory canopy is minimal (Giunta et al., 2016). In contrast, high-severity fire regimes are characterized by transitions from surface fuels into the crowns of trees, consuming a majority of overstory vegetation (Giunta et al., 2016). On the other hand, mega-fires are defined as fires covering an area exceeding 10,000 ha arising from single or multiple related ignition events (Stephens et al., 2013; Linley et al., 2022). Megafires occur in a range of biomes, but were most frequently reported in temperate-boreal forests and woodlands biomes (Linley et al., 2022). Species that grow in areas where low-severity fires occur frequently have special adaptations like thick bark (e.g., Pinus ponderosa), the ability to re-sprout after fire or store seeds on the tree until a fire occurs (Fischer et al., 2013) or geoxylic life forms in African savannas and Brazilian Cerrado (Maurin et al., 2014). Many savanna trees in Africa and eucalypts in Australian are dependent on fire for regeneration (Tng et al., 2012). These kinds of adaptations are missing in tropical rainforests and temperate forests, which do not naturally experience frequent fires (Tng et al., 2012; Fischer et al., 2013).

The primary effects of wind include damage and breakage to the tops of crowns, branches, uprooting, and snapping of trees (Giunta et al., 2016). Wind thrown trees can also create suitable habitat for bark beetles or increase the fuel load (Giunta et al., 2016). Wind damage risk increases with increase in height, and trees are known to adjust their growth pattern to their local wind environment (Jackson et al., 2021).

Remarkable among insects as disturbance agents are the bark beetle outbreaks causing large-scale transformations of forest landscapes in the northern hemisphere (Giunta et al., 2016; Hlásny et al., 2021). The European spruce bark beetle (Ips typographus) is the primary outbreak species in Europe causing as much as 8% of all tree mortality in Europe (Hlásny et al., 2021). Its damage to Norway spruce (Picea abies), its primary host, has been historically very high in the temperate latitudes than in boreal forests (Hlásny et al., 2021). A similar trend occurs in western Canada and the USA, where tree mortality due to the mountain pine beetle (Dendroctonus ponderosae) exceeds 28 million ha (Hlásny et al., 2021). In the USA, a rise in bark beetle activity since the early 1990s has occurred across a range of forest types (Giunta et al., 2016). The Douglas-fir beetle (Dendroctonus pseudotsugae) is the primary insect pest of interior Douglas-fir (Pseudotsuga menziesii) forests (Giunta et al., 2016). Another example of insect outbreaks, is aspen (Populus tremuloides) defoliation by forest tent caterpillars (Malacosoma disstria) in North America. For the present analysis, we focussed on tree species affected by the bark beetles and forest tent caterpillar outbreaks.

For this analysis we chose allometric scaling between (1) stem height and stem diameter (H–D); (2) crown radius and stem diameter (CR–D) or stem height (CR–H); (3) crown diameter with stem diameter (CD–D), (4) crown depth with stem diameter (Cdep–D), (5) crown volume with stem diameter (Cvol-D); (6) crown depth with crown diameter (Cdep-CD); (7) aboveground biomass and stem diameter (A–D); and (8) belowground biomass and aboveground biomass (A-B).

The distributions of the H–D, CR–D, and CR–H allometry parameters revealed significant overlap between Gymnosperms and Angiosperms (Figure 1 and Supplementary Figure 1). When estimated at either the family or genus levels, the median values of the exponents of Gymnosperms and Angiosperms were also not significantly different from each other (Supplementary Table 1). Similarly, there was significant overlap between the distributions of exponents in the different climate zones (Figures 2A, B). The median values also did not significantly differ among climate zones (Supplementary Table 1). Except for the extreme values and outliers, the exponents predominantly fell within the range of theoretical values in Table 1. The distributions of exponents were more consistent with MST (β = 2/3) than the geometric (β = 1) or stress similarity (β = 1/2) hypothesis (Figures 1, 2). Values <1/2 or >1 were mostly artifacts of small sample size (N) estimation (Figure 8 and Supplementary Figure 2). All of the large outliers were associated with N < 10. When sample sizes exceeded 60, estimates of the allometry parameters [both β and ln(α)] fell in a narrower range than when N < 30 (Figures 8A, C, E and Supplementary Figure 2). The LOESS trend lines show that the exponents tended to stabilize around the MST exponent as sample sizes increased (Figures 8B, D, F). The minimum diameter in the sample also had significant effects on the H–D, CR–D, and CR–H exponents (Supplementary Figure 3). Across climatic zones, the H–D exponent was smaller when the minimum D was >5 cm compared to samples that include smaller trees D < 5 cm. The distributions of exponents of H–D,CR–D, and CR–H allometry calculated at the species level are shown here using Pinus and Quercus species as examples (Supplementary Figure 4). As in the family- and genus-level analysis above, the distributions of exponents were more consistent with MST than the geometric and stress similarity hypotheses.

When data were analyzed after aggregating at different scales, the OLS estimate of the exponent was significantly smaller for angiosperms (CI: 0.56, 0.57) than the theoretical value of 2/3 but larger than 1/2. The exponent of gymnosperms (CI: 0.65, 0.67) was consistent with MST (Table 2). On the other hand, the RMA estimates for both angiosperms (CI: 0.75, 0.77) and gymnosperms (CI: 0.88, 0.89) were significantly larger than 1/2 or 2/3 but smaller than 1. Further analysis of the H–D and CR–D allometry did not reveal systematic variations in the exponent with taxonomic levels (Table 2). In the majority of cases the exponents were closer to 2/3 than 1/2 or 1. Unlike in Angiosperms, the exponents of the CR–H allometry in Gymnosperms were mostly biased downward (Table 2). The CR–H exponent was precisely biased downward in P. sylvestris in temperate zones (Table 2). To see whether different model fitting techniques will change the outcome, we compared non-linear regression, OLS, RRA, RMA, MA, and LMM analysis of the CR–H scaling in P. sylvestris in temperate zones. The non-linear regression exponent of 0.25 (CL: 0.24, 0.27), OLS exponent of 0.26 (CI: 0.25, 0.28), RRA exponent of 0.32 (CI: 0.30, 0.34), and LMM exponent of 0.45 (CI: 0.43, 0.47) were significantly smaller than the MA exponent of 0.61 (CI: 0.58, 0.65) and RMA exponent of 0.85 (CI: 0.84, 0.87) for the same species. The discrepancies between Model I (OLS and RRA) and Model II (RMA and RMA) estimates are indicative of large errors in measurement of H.

In-depth analyses of H–D, CR–D, and CR–H allometry in P. sylvestris revealed how study-level estimates of the allometry exponents can be biased by a combination of small sample sizes, sampling bias and attenuation bias (Table 3). Attenuation bias is more evident in CR–H and H–D than in CR–D allometry indicated by the large discrepancies between OLS and RMA estimates of the exponent (Table 3).

Comparisons of the allometry exponents of the H–D, CR–D, and CR–H scaling in genera with high and low disturbance regimes are summarized in Table 4. All genera in the high disturbance regime had significantly higher H–D exponents than those in low disturbance regimes. On the other hand, the CR–D and CR–H scaling exponents did not show consistent differences between trees growing in high and low disturbance regimes (Table 4). Almost all of the exponent were significantly higher than the stress similarity exponent and lower than the geometric similarity exponent, but closer to the MST exponent.

Comparison of the exponents of hyper-emergent trees species with medium and short-statured species did not reveal systematic variations (trends) with tree size (Table 5). The H–D and CR–D exponents of hyper-emergent species were in the same range as those of the medium to short statured species (Table 5). For example, the H–D exponents of S. sempervirens and S. giganteum (0.67) are statistically indistinguishable from P. tremuloides (0.64), V. seyal (0.65), and C. collinum (0.69) although they occur in different climate zones, biomes and disturbance regimes (Table 5). Similarly, the CR–D exponent of P. menziesii (0.57) is statistically indistinguishable from Picea abies (0.57) and C. mopane (0.63) despite the differences in the climate zones, biomes and disturbance regimes (Table 5). The H–D and CR–D exponents of almost all species in Table 5 did not conform to the geometric similarity hypothesis (β = 1). With the exception of E. regnans, the majority of the exponents fell within the range of values assumed in stress similarity and MST. In-depth analysis at the site level data indicated that curious case of E. regnans arose due to aggregation of data of varying quality (Supplementary Table 2). The probable causes of the bias become clear when the data from each Reference ID were analyzed separately (see Supplementary Table 2 for details). These problems could not be remedied by application of more sophisticated analyses such as robust regression (Supplementary Table 2).

Comparisons of tropical forests with savanna biomes in terms of CD–D, Cdep–D, Cdep–CD, and Cvol–D also revealed broadly overlapping distributions of the exponents (Figure 3). The distributions of exponents were closer to the MST exponents than the geometric similarity and stress similarity exponents (Figure 3). Comparison of the point estimates (and CIs) of the exponents of CD–D also did not reveal systematic variations with biomes or continents (Supplementary Table 3). Where such differences exist, differences in sample size and sampling bias seem to play a role. For example, the sample size for American forests was 113 times larger than for American savannas. The diameter range reflects sampling bias toward larger trees in American forests than in the savannas. Therefore, the comparisons among the forest and savanna biomes are confounded by the differences in the diameter ranges and sample sizes. Asian forests and Australian savannas had exponents biased downward partly due to aggregation of data of varying quality. Graphic comparison of CD–D allometry in African and Australian forests and savannas reveals further insights into the degree of sampling bias (Figure 4). Sampling bias toward large trees is more evident in data from Australian forests (Figure 4C) and African savannas (Figure 4B). In four out of the seven African countries the minimum D was 10 cm, while in the remaining three countries the minimum D was 0.2–6.3 cm. These differences are reflected in the two groups of data in Figure 4A. Measurement errors are also evident in Figures 4C, D. Close examination also revealed that data from one site in Asia (Mala_8) was weighing down the exponent estimated for the whole continent (Supplementary Table 3). The Mala_8 site had two populations with non-overlapping distributions, one of which was of very poor data quality (Supplementary Table 3).

Statistical tests of differences in Cdep–D, Cvol–D, and Cdep–CD allometry between tropical forest and savanna biomes are summarized in Supplementary Table 5. With the exception of Cvol–D allometry, tropical forests and savannas did not significantly differ in their exponents. In the case of Cvol–D, the OLS estimate of the exponent was biased downward in tropical forests than in savannas, but the RMA estimates are comparable (Supplementary Table 5).

The distributions of the A–D exponents in the different climate zones were significantly overlapping (Figure 5A). The median value of the exponents – 2.32 (CI: 2.25, 2.36) in boreal, 2.35 (CI: 2.34, 2.39) in temperate, 2.40 (CI: 2.34, 2.47) in subtropical and 2.42 (CI: 2.41, 2.48) in tropical zones – are statistically non-significant. The median values of the A–H exponents in the tropical biomes, where data were available, was 3.45 (CI: 3.18, 3.77). In all cases, the distributions of exponents were more consistent with the stress similarity and MST exponents than the geometric similarity exponent (Figure 5A). Tropical, subtropical and temperate regions had a large number of outliers of the A–D exponent. The total sample size available for analysis had significant effects on the A–D exponent (Figure 5C). There were large outliers when N < 30 than N > 60 trees (Figure 5C). The LOESS trend lines show that the exponent tended to approach the MST exponent as sample sizes increased (Figure 5C).

The 95% CI of the B–A exponents shows that in the majority of cases the exponent conforms with isometric scaling. The exponents also did not systematically vary with climate zones (Table 6). Except in a few cases, the RMA estimate of the B–A exponent in the different climate zones did not significantly differ from 1 (Table 6).

The analyses did not reveal significant trends in the allometry parameters of the H–D, CR–D, and CR–H scaling along the divergence times of Gymnosperm and Angiosperm families (Figure 6). Within Angiosperms the H–D allometry exponent of the more recent (∼29 MY old) family Bignoniaceae (β = 0.68) was not significantly different from those of the older (113–125 MY old) families Magnoliaceae (β = 0.68) and Winteraceae (β = 0.70). The H–D exponent of the more recent (∼24 MY old) Angiosperm family Staphyleaceae (β = 0.53) was also not significantly different from those of the older (∼236 MY old) Gymnosperm families Araucariaceae (β = 0.54), Podocarpaceae (β = 0.56), and Taxaceae (β = 0.44). The CR–D and CR–H allometries followed trends similar to H–D.

The analyses also revealed significant negative correlations between ln(α) and β in the allometries analyzed at the different scales (Figure 7 and Supplementary Figure 4). As a result, when ln(α) was biased upward, β was biased downward. We noted significant positive correlations between the CR-D and CR–H exponents (r = 0.647; P < 0.0001). Similarly, we observed a positively correlation (r = 0.543; P < 0.0001) between H–D and A–D exponents (Supplementary Figure 5) indicating tight co-variations among D, H, CR, and biomass.

Our analyses have uncovered a tight co-variation among stem height, diameter, crown dimensions and tree biomass, consistent with allometric theory (e.g., Niklas and Spatz, 2004; West et al., 2009). The various analyses did not reveal evidence of systematic variations in the H–D, CR–D, and CR–H allometry with taxonomic level, evolutionary time, climate zones or biomes. These findings are consistent with the allometry constraint hypothesis, which posits that the evolutionary divergence of traits is restricted by integrated growth regulation, and thus β remains constant at macro-evolutionary timescales (Voje et al., 2013). The results are also in agreement with those of Blanchard et al. (2016), who found no significant variation in CA–D allometry across five biogeographic areas. According to Blanchard et al. (2016), the stability in CA–D allometry suggests that some universal constraints are sufficiently pervasive to restrict the variation of the exponent in a narrow range.

Unlike the H–D and CR–D scaling, the CR–H allometry was more variable probably because both H and CR are associated with greater measurement error than D (Ducey, 2012; Calders et al., 2015; Blanchard et al., 2016). Contrary to recent reports (e.g., Shenkin et al., 2020; Panzou et al., 2021), our analyses did not reveal significant differences in crown allometry between savanna and forest biomes. Our analysis also casts doubt over the notion that the H–D allometry depends on the context in which a tree grows and hence the exponent varies widely [references cited in Osunkoya et al. (2007), Watta and Kirschbaum (2011); Fortier et al. (2015), Hulshof et al. (2015), and Motallebi and Kangur (2016)]. Some of the past studies did not have adequate sample sizes and/or did not adequately control for differences in sample size and size-frequency distributions resulting in conflation of statistical artifacts with violations of allometric relationships. For example, Osunkoya et al. (2007), used small sample sizes (10–39 trees) and found wide variability in the exponent of the H–D scaling in 22 species. Watta and Kirschbaum (2011) reported that the exponents of the H–D relationship ranged between 0.73 and 1.43 across 84 plots, and this clearly violates the assumptions underlying allometric relationships. However, close examination of their report reveals differences in both the diameter range (D < 10 cm vs. D > 10 cm) and sample size between the two groups. Similarly, Fortier et al. (2015) reported differences in H–D exponents between moderate and high fertility sites. However, the diameter range of trees from the moderate fertility site (D = 10–27 cm) was smaller than the high fertility site (D = 20–38 cm) besides the very small sample sizes used (12) at each site. Even if sample sizes were large, the difference in the diameter range alone could have resulted in the observed differences in exponents. Here, we do not claim that the allometry exponent is static, but its variability is conflated with sampling bias and low statistical power due to small sample sizes. Unlike previous studies, our use of parameter distributions reveals the reality better than point estimates and their 95% CIs which are sensitive to differences in sample sizes and sampling bias. Our inferences are also based on disaggregated data and hence the agreement with allometry theory.

In previous studies, inferences have been solely based on point estimates and their 95% CI derived from aggregated data. For example, in Panzou et al. (2021) complex LMMs were implemented using aggregated data, which are prone to violations of LMM assumptions. We suspect this is the source of the counterintuitive values of the exponent reported in Panzou et al. (2021) for the American forests, American savannas, Asian forests, Australian forests and Australian savanna (Supplementary Table 3). The problem with aggregated data is that they are prone to contamination with poor quality data as is the case with E. regnans (Supplementary Table 2) and Asian forests (Supplementary Table 3). Data aggregation is known to leave analyses vulnerable to aggregation bias or ecological fallacy, i.e., drawing of false inferences about individual behavior on the basis of aggregate level statistics (Pollet et al., 2015; Salkeld and Antolin, 2020).

In a nutshell, the data does not support our hypothesis that the allometry exponents systematically vary with taxonomic levels, divergence time or climate zones. The data does not also support our hypothesis that trees adapted to different disturbance regimes have different allometry exponents. Although disturbance regimes are said to be critical in driving local to regional-scale variance in tree allometry relationships (Wenyan et al., 2022), the evidence seems to be weak. For example, Wenyan et al. (2022) found crown allometry exponents that conform to theoretical predictions in gap forest sites (created by cutting trees), and greater deviations from theoretical predictions in unmanaged forest sites. Indeed, their finding is contrary to the common notion that forests that have undergone disturbance deviate substantially from theoretical predictions.

Our initial hypothesis that hyper-emergent and short-statured tree species follow different allometric trajectories was also not supported. According to Osunkoya et al. (2007), H–D scaling in understory species follows a geometric similarity, while mid-canopy and emergent species follow stress–elasticity similarity. Close examination of their reports revealed that some of the heterogeneity in exponents is an artifacts of small sample sizes and inferences based on point estimates.

Despite the wide variability in sample size and sampling conditions, the A–D exponents did not significantly vary with climate zones. They were also closer to the theoretical values of 5/2 and 8/3 for A–D allometry and 3–4 for A–H allometry. This is consistent with Zianis and Mencuccini (2004) who noted that there is a general convergence of the scaling exponents despite the multitude of factors affecting tree growth in different sites. Allometry theory (Enquist et al., 1999; Enquist and Niklas, 2002; Niklas and Enquist, 2002) and empirical results (e.g., Kuyah et al., 2013; Sun et al., 2022) also predict the same exponent for the allometric relationships between the different biomass components (i.e., stem, branch, foliage and root) and D. For example, in our re-analysis of data from Kuyah et al. (2013), the exponent was 2.50 for stem biomass to D allometry, 2.68 for branch biomass to D allometry, 1.87 for foliage biomass to D allometry and 2.45 for root biomass to D allometry. Similarly, in Sun et al. (2022) the exponent was 2.176 for stem biomass to D, 2.334 for branch biomass to D, 1.578 for foliage biomass to D and 2.063 for root biomass to D allometry although biomass allocation to stem, branch, foliage, and root biomass significantly differed. The exponent was also estimated at 2.196 for total woody biomass to D allometry, 2.235 for aboveground woody biomass to D allometry (Sun et al., 2022). The closeness of these figures indicate that the different biomass components follow the same scaling rules and allocation patterns.

Our earlier review (Sileshi, 2014) and the present analysis reveal that the exponents reported in the literature are biased downward partly due to the small sample sizes used in individual studies. Duncanson et al. (2015) showed that the use of small sample sizes in allometric equations can result in positive bias of ∼70% in the average site-level biomass estimates. In the past, the limiting factor has been the time and resources needed for destructive sampling of a large number of trees (Duncanson et al., 2015). With the increasing availability of portable ground-based LiDAR (Calders et al., 2015), data for large sample sizes can be acquired quickly, including samples of very large trees for which destructive sampling would be logistically impractical (Duncanson et al., 2015).

The power function has been reported to perform equally or even outperform other models in many situations (e.g., Sileshi, 2014; Sun et al., 2022; Wenyan et al., 2022). However, the debates over perceived variability of β has been discouraging practitioners from using simple power law models. As a result, empirical models of statistically dubious quality continue to proliferate the biomass estimation literature (Sileshi, 2014). Here we have shown that the variability in β is mostly an artifact of small sample sizes and sampling bias toward large trees. In forest inventories, it is a common practice to measure trees above a certain stem diameter, e.g., D > 5 cm or even >10 cm. Biomass estimation models developed using trees with D > 10 cm tended to overestimate the mean diameter (X¯) and ln(α) resulting in underestimation of β (see Equation 4 below). Therefore, we strongly recommend destructive sampling of trees including smaller stem diameter (e.g., 2.5 cm) during the development of biomass estimation models. We also recommend sampling roughly equal number trees from each diameter classes to get a truly representative sample of the size-frequency distribution of the target population [see Kuyah et al. (2013)]. When these statistical problems are remedied, the power law model could provide a more convenient tool for predicting forest biomass and carbon stocks at different scales. A model based on a tested and established theory is more likely to be robust to new information than those purely based on observed patterns and correlations, which may prove unstable when new information emerges.

Our analysis indicates that belowground biomass scales with aboveground biomass isometrically (β = 1) regardless of the context. This is consistent with allometric theory (Cheng and Niklas, 2007). This makes allometric models a more powerful tool than the traditional use of root-to-shoot ratios to predict belowground biomass. Root to shoot ratios often vary across biomes, vegetation types, taxonomic groups (Qi et al., 2019), growth stages and tree sizes (Peichl and Arain, 2007; Kuyah et al., 2013; Mašková and Herben, 2018). For example, a global analysis by Qi et al. (2019) found significantly higher root-to-shoot ratios in angiosperms than in gymnosperms. Similarly, Peichl and Arain (2007) and Mašková and Herben (2018) found decrease in root to shoot ratios with time and substrate nutrients. In contrast, a single allometric equation could predict total belowground biomass from aboveground biomass across the entire age-sequence (Peichl and Arain, 2007; Robinson and Peterkin, 2019).

Since there were not many datasets on changes in tree biomass allocation with taxonomic levels, disturbance regimes or stages of forest stand development, we were unable to investigate the B–A allometry in more detail. Biomass allocation also remains poorly documented in geoxyles, whose biomass is disproportionately found belowground. Geoxyles are plants with short-lived reproductive aerial branches and woody underground structures (xylopodia) (Maurin et al., 2014; Meller et al., 2022), and these are common in regions with frequent fires such as African savannas and the cerrado in Brazil (Maurin et al., 2014; Meller et al., 2022). Geoxyles may provide the next frontier of research in biomass allocation in biomes that experience frequent disturbance by fire.

In the literature, counterintuitive values of β have been reported and in some cases such values have been used to challenge predictions of allometry theories. For example, the MST has been challenged by a number of forest ecologist [see Shenkin et al. (2020) for details] for not providing coherent explanations for the variability in β. In the following sections we show how (1) spurious values of β can arise as statistical artifacts; (2) β varies with ln(α), and (3) the accuracy with which β is estimated depends on the sample size (N), measurement errors, the representativeness of the sample available for analysis and the regression technique used.

The various analyses in this study and earlier studies (e.g., Zianis and Mencuccini, 2004; Ducey, 2012; Sileshi, 2014; Zhang et al., 2016) have demonstrated an inverse linear relationship between the exponent and the intercept. This points to some kind of trade-off between ln(α) and β. One possible explanation is the principle of optimality in biological design (Popescu, 1998). Since natural selection leads to an economy of design, the parameter trade-off may be regarded as a manifestation of an attempt to achieve optimality. The observed co-variation between ln(α) and β is consistent with the power-law relationship in the logarithmic domain:

where X¯ is the mean of log(X) and Y¯ is the mean of log(Y). From equation 3 it follows that

This relationship implies that bias in ln(α) can result in biased estimates of β. For example, data in Figures 2, 3, 5 revealed that β values that are biased in one direction are associated with ln (α) values biased in the opposite direction. From Equation 4 it follows that an upward bias in ln(α) or X¯ will result in a downward bias in β and vice versa. Biases in X¯ typically arise from small sample sizes, sampling biases and measurement errors. For example, sampling bias toward large trees will result in an upward bias in X¯ and ln(α), therefore a commensurate downward bias in in β.

From Figures 5, 8 it is evident that β estimates tend to be biased either downward or upward when sample sizes are small. Our findings reinforce earlier reports that allometry parameters are highly sensitive to sample size (Duncanson et al., 2015). When N < 30, the 95% CI also tend to be very wide, and cover different theoretical values of β (see Table 3). This makes it difficult to distinguish between the different theoretical predictions. Our review of the published exponents indicates that site-specific biomass models were particularly based on small sample sizes. For example, over 75% of the A–D scaling exponents in the literature were estimated using N < 60, of which 43% were based on N < 30 trees. Therefore, it is not surprising that the median value of the A–D exponents were smaller than the theoretical value of 8/3 (Figure 5). To achieve sufficient statistical power, we recommend the use of a minimum of 66 sample trees when estimating parameters of the power law model following the rule of thumb N > 50+8P proposed by Green (1991).

The accuracy with which β is estimated in OLS regression depends on the accuracy with which the X and Y variables were measured. Measurement errors cause a pervasive problem known as attenuation bias or regression dilution (Maroco, 2007; Hutcheon et al., 2010). OLS regression assumes that the X-variable is measured without error. In reality, tree dimensions such as D and H are measured with a great deal of error (see Supplementary methods). In the presence of measurement errors in X, estimators of the correlation coefficient (r) do not converge to their true population values (ρ) even if the sample size is infinitely large. Accordingly, β will be biased downward because it is a function of the correlation coefficient (r) between log(X) and log(Y), their variances (V_Y and V_X) and the covariance (i.e., COV_X,Y):

where σ_D and σ_Y are the standard deviation of log(X) and log(Y), respectively. Loken and Gelman (2017) showed that with large N, adding measurement error will almost always reduce the observed correlation (r) between X and Y. As such, β becomes more precisely biased toward zero as sample sizes increase (Maroco, 2007). From equation 5 it also follows that the sign of β depends entirely on the sign of r or COV_X,Y since σ_X, σ_Y and V_X are always positive. If the measurement errors in X are very large, the sign of r tends to be negative. This is probably why counterintuitive values of the exponent reported in the literature emerge.

The accuracy with which of β is estimated also depends on the representativeness of the sample of the underlying population. P. sylvestris (Table 3), Australian savannas (Figure 4) and E. regnans (see Supplementary Table 2) provide vivid examples of how sampling bias result in biased estimates of allometry parameters.

The various analyses (Tables 3, 6 and Supplementary Table 3) show that different regression techniques can lead to entirely different conclusions about the size of β for the same dataset. These differences arise due to the differences in the way in which measurement errors and outliers are handled within the different regression techniques (see Supplementary methods). At the different scales of analysis, we have shown that the RMA estimates of β are consistently larger than the OLS estimates, while LMM estimates are usually smaller than OLS estimates (Tables 3, 6 and Supplementary Tables 2, 3, 5). Elsewhere RMA was also reported to overestimate the true regression slope (Kilmer and Rodríguez, 2017). This is due to the mathematical relationship between the OLS and the RMA estimators of β. For r ≠ 0, β_RMA is the ratio between the β_OLS and r:

β_RMA is also the ratio of SD_Y to SD_X (Sokal and Rohlf, 1995):

This implies that β_RMA will be biased upward if r is biased toward 0. Therefore, RMA can give a false impression of isometry when in reality the actual relationship is allometric.

Although LMMs provide a convenient framework to account for factors associated with site quality, climate, or stand history, in our analyses it yielded β estimators that were biased downward (Supplementary Table 5). We have also noted this problem in earlier studies (e.g., Panzou et al., 2021). It must be noted that LMMs are sensitive to violations of various assumptions (Schielzeth et al., 2020) and imbalances in sample size between categories. Allometry parameters in LMM can be biased due to these violations, and therefore empirical estimates should not be taken at face value.