Royal Society Publishing

Species traits and the form of individual species–energy relationships

Karl L Evans, Sarah F Jackson, Jeremy J.D Greenwood, Kevin J Gaston

Abstract

Environmental energy availability explains much of the spatial variation in species richness at regional scales. While numerous mechanisms that may drive such total species–energy relationships have been identified, knowledge of their relative contributions is scant. Here, we adopt a novel approach to identify these drivers that exploits the composite nature of species richness, i.e. its summation from individual species distributions. We construct individual species–energy relationships (ISERs) for each species in the British breeding avifauna using both solar (temperature) and productive energy metrics (normalized difference vegetation index) as measures of environmental energy availability. We use the slopes of these relationships and the resultant change in deviance, relative to a null model, as measures of their strength and use them as response variables in multiple regressions that use ecological traits as predictors. The commonest species exhibit the strongest ISERs, which is counter to the prediction derived from the more individuals hypothesis. There is no evidence that predatory species have stronger ISERs, which is incompatible with the suggestion that high levels of energy availability increase the length of the food chain allowing larger numbers of predators to exist. We find some evidence that species with narrow niche breadths have stronger ISERs, thus providing one of the few pieces of supportive evidence that high-energy availability promotes species richness by increasing the occurrence of specialist species that use a narrow range of resources.

Keywords:

1. Introduction

Much of the marked spatial variation in species richness can be explained by the availability of environmental energy (Hawkins et al. 2003; Gaston & Spicer 2004; Pimm & Brown 2004). Such total species–energy relationships have been documented in different geological time periods, and across a wide variety of taxa and habitats (Roy et al. 1998; Crame 2001; Hawkins et al. 2003; Pimm & Brown 2004). Much effort has been invested in describing the form of these relationships, and how they vary with spatial scale (for reviews see Waide et al. 1999; Mittelbach et al. 2001, 2003; Whittaker & Heegaard 2003). Positive relationships tend to dominate at broad spatial scales and a number of potential causal mechanisms of these patterns have been identified, but their relative importance remains poorly understood (Currie et al. 2004; Evans et al. 2005a).

The more individuals hypothesis (MIH) suggests that in areas with high levels of environmental energy, as measured by plant productivity, consumers can maintain larger populations that reduce their extinction risk, thus elevating species richness (Wright 1983). Extinction risk is related to population size by a negative decelerating function (Lande 1993), thus rare species should respond more strongly to increasing energy availability than common ones (Evans et al. 2005b). Second, the niche breadth hypothesis suggests that, in areas with high plant productivity, resources will be sufficiently abundant to promote the occurrence of species that specialize on a few resource types and elevate species richness directly and also possibly by reducing competitive exclusion (Abrams 1995; Vázquez & Stevens 2004; Evans et al. 2005a). Third, the niche position hypothesis suggests that high-energy areas are the only ones that contain relatively scarce resources in sufficient abundance to maintain viable populations of the niche position specialists (sensu Shugart & Patten 1972) that use them (Abrams 1995; Evans et al. 2005a). Thus, both the niche breadth and niche position hypotheses predict that specialist species will respond more strongly to energy availability than less specialized ones. Fourth, the trophic level hypothesis assumes that the transfer of energy between trophic levels is inefficient, so the number of levels may be regulated by the amount of energy at the base of the food chain (Oksanen et al. 1981; Fretwell 1987). Longer food chains in highly productive areas should thus enable greater numbers of predatory species to occur, and predators would therefore respond more strongly to energy than herbivores or omnivores. While there is strong evidence that food chain length is only limited by environmental energy availability at much lower levels of environmental energy availability than occur in many regions (Post 2002), we include this hypothesis in the interests of thoroughness. Fifth, the thermoregulatory load hypothesis entails higher temperatures reducing the amount of energy that endotherms require for thermogenesis (Turner et al. 1988). Smaller individuals, other things being equal, are more susceptible to heat loss than larger ones, thus the thermoregulatory load hypothesis predicts that species with smaller body sizes will exhibit stronger responses to temperature (Lennon et al. 2000).

Predictions can be drawn from these mechanisms concerning which types of species should exhibit the strongest species–energy relationships. These predictions can be tested by assessing how the richness of different species groups, such as common and rare species, respond to energy availability (Jetz & Rahbek 2002; Ruggiero & Kitzberger 2004; Vázquez & Gaston 2004; Evans et al. 2005b,c). An alternative and complementary approach, that as far as we are aware has not been previously used, is based on the principle that species richness results from the summation of occurrence data for each species in the focal assemblage. Here, we assess how the relative strength of individual species responses to energy availability, in terms of their presence/absence, is related to their ecological characteristics. Such an approach has great potential to reveal the underlying mechanisms driving spatial variation in biodiversity (Marquet et al. 2004).

We are not aware of any other study that has systematically investigated the nature and determinants of the strength of individual species–energy relationships (ISERs). In this manuscript, we illustrate this technique using the British breeding avifauna as a case study.

2. Material and methods

(a) Avian distribution and energy data

We used the summer breeding distribution of birds in Great Britain recorded in April–July 1988–1991 (Gibbons et al. 1993). These data are species presence/absence (combined over the four years) at a resolution of 10×10 km quadrats on a continuous grid. Fieldwork was coordinated by a network of regional organizers and undertaken by experienced volunteer ornithologists. Data are based on timed visits, of two hours' duration, to at least eight tetrads (2×2 km quadrats) within each 10 km quadrat and supplemented with additional records to ensure that the species lists were as complete as possible. We excluded marine species and many vagrants, but retained self-sustaining introductions and those species that bred in each season between 1988 and 1991, giving a total of 189 species. Some filtering was performed on the distributional data; 10 km quadrats that contained less than 50% land were excluded, leaving 2262 quadrats.

There are two main types of environmental energy that are relevant in the present context (Evans et al. 2005a). Solar energy is the amount of solar radiation reaching the Earth's surface and can be measured by surrogates, such as temperature; productive energy is the amount of resources available for consumers, which can be measured by plant productivity or its surrogates, such as the normalized difference vegetation index (NDVI; Boelman et al. 2003; Kerr & Ostrovsky 2003). We used both types of energy metrics in this study. We obtained mean monthly temperature data that were derived from meteorological recording station readings for the period 1961–1990 using surface interpolation techniques (Barrow et al. 1993). We obtained NDVI data from the NOAA/NASA Pathfinder AVHRR Land Data Set (see http://www.ciesin.org/). These data were collected between 1981 and 2001 at a spatial resolution of a 0.1° latitude/longitude grid, approximately equivalent in Britain to an 8×8 km quadrat. Daily readings are converted to maximum values for each 10-day period, which markedly reduces the effects of cloud cover. From these we calculated mean monthly NDVI values and then used a geographical information system to re-project these data at a 10 km resolution, compatible with our avian distribution data. For both temperature and NDVI, we calculated a mean annual measure of energy availability and a mean summer value from the monthly averages for May, June and July.

(b) Species traits

Data on breeding population size and body mass were taken from the compilation in Gaston & Blackburn (2000), with additional data, for Columba livia, from Greenwood et al. (1996). Species were categorized as (i) long-distance migrants if most of their breeding populations wintered outside Britain, in sub-Saharan Africa for most species, or (ii) short-distance migrants if the majority of the population wintered within Britain but in a different location from that in which they bred, or (iii) residents. Species were classified by trophic level (herbivores, omnivores and predators) on the basis of combined adult and chick diets using data on the major food items in each species' diet from Cramp et al. (1977–1994). Niche breadth and niche position data were derived from a previous study, for 85 species, from a canonical correspondence analysis based on avian abundance data and 34 environmental variables (Gregory & Gaston 2000). To reduce the skew in their frequency distributions, data on population size, niche position and body mass were logarithmically transformed to base 10.

(c) Statistical analyses

Numerous methods are available to model species distributions, the suitability of each being dependent on the aims of the study. We used logistic regression, which focuses on recording the general trend of a species' presence/absence response to predictor variables (Segurado & Araújo 2004). Using SAS we constructed two binary logistic regression models for each species using either temperature (°C) or NDVI as measures of environmental energy availability.

From each ISER model, we obtained two measures of the strength of the relationship: (i) the linear slope of the relationship, and (ii) the proportional change in deviance relative to a null model, i.e. one that lacked predictors. To reflect the directionality of the ISER, the change in deviance was multiplied by −1 if the slope was negative (see the electronic supplementary material, figure 1). The first of these metrics is a measure of how rapidly a species' pattern of occurrence changes with energy, and the second of how closely the model fits the observed data. There is a weak correlation between the two (temperature-based models r2=5%; NDVI-based models r2=18%; see also the electronic supplementary material, figure 1a,b). The use of square energy terms in the logistic regressions resulted in relatively little change in the fit of the ISERs (the mean percentage change in deviance decreased by 2.2 and 1.7% for temperature and NDVI models, respectively). We thus present results based on ISERs that were constructed using only a linear energy term. We used summer (May–July), rather than annual temperature- and NDVI, as the former were better predictors, as measured by change in deviance, of species distributions for the majority of species (84 and 69% of species, respectively, with temperature- and NDVI-based models).

Data on species distributions tend to be influenced by spatial autocorrelation, which invalidates the assumption of independent errors made by classical statistical tests of association and distorts estimates of the slope of the relationship between the response and predictor (Cressie 1991; Lennon 2000; Legendre et al. 2002). When analysing species distributions, autologistic regression is one of the more frequently used approaches for taking spatial autocorrelation into account (Augustin et al. 1996; Segurado & Araújo 2004). This uses autocovariate terms, which summarize the presence/absence data for a species in neighbouring cells, as additional predictors of the presence/absence of that species in a focal cell. As environmental factors will influence the presence/absence data in neighbouring cells, autocovariate terms contain information on both spatial and environmental factors. Therefore, in the context of determining the strength of the relationship between occurrence and the environment, their use provides a conservative approach.

We implemented autologistic models following the methodology of Gu & Zhu (2001) and He et al. (2003). We used a second-order autologistic regression model, in which the response at a focal site depends upon the responses at the focal site's eight immediate neighbouring sites (those to the north, northeast, east, southeast, south, southwest, west and northwest). Thus, a site in the interior of the study area will have eight neighbours, whereas sites on the boundary will have fewer neighbours. The autologistic regression model was estimated using the two-stage Markov chain Monte Carlo stochastic approximation method (MCMC-SA; see the electronic supplementary material, appendix A). The key feature of this model is that it does not use standard methods of parameter estimation. The latter depends upon being able to evaluate the likelihood function and its derivatives, which is not possible given the intractability of the computed conditional probabilities (see the electronic supplementary material, appendix A). In stage I, a sequence of large gain constants is used to estimate quickly the feasible range of parameter estimates, which are refined in stage II using an optimal stochastic approximation procedure. The models were run using C with codes obtained from F. L. He and H. T. Zhu.

Autologistic models converged for 155 and 146 species, respectively, when temperature and NDVI were used as predictors. We ranked species by the strength of their responses to environmental energy, from the strongest negative response to the strongest positive response, and compared the rankings that were obtained from independent error and autologistic models. There is a significant correlation between the spatial and non-spatial parameter estimates (P<0.0001 for both temperature-based, n=155, and NDVI-based models, n=146). However, the correlation coefficients are sufficiently low (0.7 in both cases), and the mean differences in ranks are sufficiently large (21 in both cases), to warrant a conservative approach that uses data from both spatial and non-spatial models when analysing the correlates of the slopes of ISERs. We did not calculate the change in deviance from autologistic models due to the time-consuming nature of this process, which was estimated to take nine days per species on a Pentium 4 PC with a 2.8 GHz processor, giving a total processing time of approximately four and a half years for all the species whose models converged.

We used our six measures of the strength of ISERs (the slope and deviance from non-spatial models, and the slope from autologistic models, with each metric being available for both temperature- and NDVI-based ISERs) as response variables in multiple regression models. The latter used ecological traits (log population size, log niche position, niche breadth, log body mass, and trophic level) to predict the strength of ISERs. Species that occupy either a small or large proportion of quadrats are less likely to show stronger responses to energy than those occupying an intermediate number, for purely statistical reasons. We therefore calculated an ‘information index’ for each species as p(1−p), where p is the proportion of quadrats it occupies, which reaches the highest value when an intermediate proportion of quadrats is occupied, and used this as an additional predictor of the strength of ISERs. Initial exploration of the data indicated the presence of simple non-linear relationships between the strength of ISERs and log-transformed population size, log-transformed niche position and niche breadth, so the square terms of these variables were also included in analyses.

We used general linear modelling to construct the multiple regression models of the strength of ISERs. Following Johnson & Omland (2004), we used an information-theoretic approach to model building. We calculated the probability that each model provided the best fit to the data, the model weight, using the Akaike information criterion (AIC). This was first done by confining the dataset to the 85 species for which niche position and niche breadth data were available and constructing all possible models given the above list of predictor variables (nine variables and 511 models). We repeated this analysis using data from all species, but excluding niche breadth, niche position and their square terms as predictors, to test if additional predictors were included in the best-fitting models. We present results from analyses that included species with both negative and positive responses to environmental energy availability, but restricting analyses to the latter gave qualitatively similar results.

Phylogenetic generalized least squares (PGLS) was used to test for phylogenetic non-independence in our data (Grafen 1989; Martins 1999; Garland & Ives 2000). PGLS explicitly incorporates the expected covariance among species into a statistical model fit by generalized least squares. The correlation between error terms is altered to reflect the degree of phylogenetic relatedness among the species to which they relate. The PGLS approach was implemented in R (Ihaka & Gentleman 1996) using the APE (analysis of phylogenetics and evolution) package (Paradis et al. 2004) and code written by R. P. Duncan. PGLS requires a hypothesis about the phylogenetic relatedness of the species analysed. We assumed that the species in our dataset were related according to the phylogeny of Sibley & Ahlquist (1990), with generic relationships as described in Sibley & Monroe (1990). Since we did not know all the branch lengths in the phylogeny, we repeated the analysis with two different assumptions about branch length models. In the first, all branches in the model were set equal. In the second, branch lengths were set to be proportional to the number of taxa below each node in the phylogeny. Branch lengths were calculated using TreeEdit v. 1.0a10. The PGLS approach generates a metric (λ), which varies from 0 (no phylogenetic autocorrelation) to 1 (strong phylogenetic autocorrelation).

3. Results and discussion

There was marked variation in the strength of ISERs, with the majority of species responding positively to energy availability (65 and 60% for temperature and NDVI, respectively). Most species responded more strongly to temperature than to NDVI, with the respective mean changes in deviance being 20 and 8% (see the electronic supplementary material, figures 1 and 2). The ability of temperature to explain spatial variation in bird distributions is, on average, twice as great as its ability to explain spatial variation in the distribution of British ground beetles measured at the same spatial scale (Eyre et al. 2005). Assuming that greater under-recording of beetles than of birds does not influence this result, it is interesting that a group of endotherms appears to respond to temperature more strongly than a group of ectotherms. This pattern may, however, arise because relationships between ground beetle occurrence in Britain and temperature have been more strongly disturbed by human activities than is the case for breeding birds; ground beetles are particularly sensitive to factors such as habitat availability and landscape structure (Rainio & Niemelä 2003), but it is unclear if they are more sensitive to these factors than birds. Variation in the effective spatial scale of analyses may also influence our findings; while both sets of analyses are conducted at a resolution of 10 km quadrats, such a spatial grain is, compared to the focal species home ranges, much larger for beetles than for birds.

PGLS analyses of each of our six response variables demonstrated that, in all cases, λ was not significantly different from zero, thus indicating that phylogenetic correlation is low in these data and that analyses that treat species as independent data points are justified, and we thus adopted the latter approach. Multiple regression analyses may be influenced by collinearity between predictors, which can be assessed by calculating tolerance values (for predictor xi, tolerance=1−r2, with the latter being calculated from the regression of xi against all remaining predictors). Tolerance values lower than 0.1 indicate that collinearity has a marked influence on the results of multiple regression analyses (Quinn & Keough 2002). The lowest tolerance value in our data was 0.28, for log niche position, and therefore we consider our analyses to be free of major concerns regarding collinearity.

Two general trends emerge from the multiple regression analyses. The ability of ecological predictors to explain interspecific variation in the strength of ISERs is weaker when such relationships are constructed using autologistic models rather than logistic models that ignore spatial autocorrelation (tables 1 and 2). This was expected given the conservative nature of the autologistic approach, which results from the autocovariate terms containing information on both space and the environment, and concurs with other studies of the relationships between species occurrence patterns and environmental factors (Betts et al. 2006). Also, while the information index is often retained in the best-fitting models, species for which more information is available do not have stronger ISERs.

View this table:
Table 1

Predictors of the strength of individual species–energy relationships among the 85 species of breeding British birds for which data on niche breadth and niche position were available. (Species–energy models were constructed using either temperature or NDVI as a measure of environmental energy availability, either without taking spatial autocorrelation into account (logistic) or with taking it into account (autologistic). The strength of such relationships is measured either by the slope of the relationship or the change in deviance, relative to a null model. An information-theoretic approach to model building is adopted, with smaller Akaike information criterion (AIC) values indicating a better fit to the data. The model weight is the probability that the model provides the best fit to the data out of all possible models that could be constructed, given the predictors used (those in the table together with the square of niche breadth). Positive effects, ++++P<0.0001, +++P<0.001, ++P<0.01, +P<0.05, +n.s. P>0.05; negative effects, −−−−P<0.0001, −−−P<0.001, −−P<0.01, P<0.05, −n.s. P>0.05; factors **P<0.01, *P<0.05, *n.s. P>0.05.)

View this table:
Table 2

Predictors of the strength of individual species–energy relationships among the entire assemblage of British breeding birds. (The methods used were identical to those described in the legend for table 1, except that niche breadth and niche position data were not used as predictors. Positive effects, ++++P<0.0001, +++P<0.001; negative effects, −−−−P<0.0001, −−−P<0.001, −−P<0.01, P<0.05, −n.s.P>0.05; factors *n.s.P>0.05.)

Before discussing the specific results arising from the multiple regression analyses, it is useful to discuss their interpretation. We do so using niche breadth as an example, although similar statements could arguably be constructed with other predictors in mind. According to the niche breadth hypothesis, a specialist species with a small niche breadth will only occur in high-energy areas. In contrast, generalist species will occur in areas with a much wider range of energy availability. Their distribution may be weakly correlated with energy availability, or even totally uncorrelated with energy if they occur at all energy levels; alternatively, generalists could occur in all cells, except those with the very lowest levels of energy. Generalist species may thus exhibit a mixture of strong and weak ISERs and specialist species will always exhibit strong ISERs. Such patterns will give rise to a constraint relationship between the strength of ISERs and niche breadth, albeit one with an overall negative slope. Therefore, while a negative relationship between the strength of ISERs and niche breadth is strong evidence for the niche breadth hypothesis, the lack of such a relationship arguably provides less conclusive evidence against such a hypothesis).

When data from all species are used, and thus niche breadth and position are not included as predictors, the best-fitting models of the strength of ISERs generally do not include additional predictors to those selected in the best-fitting models constructed using data from the smaller number of species for which niche breadth and position data are available (cf. tables 1 and 2). The exception to this is a negative relationship between the strength of single species–temperature relationships and body size when data from all species are analysed (table 2). This relationship is expected from the thermoregulatory load hypothesis (Turner et al. 1988), which predicts that, because smaller species are more susceptible to heat loss than larger ones, the occurrence of smaller species should exhibit stronger responses to temperature than larger ones. The relationship between the strength of ISERs and body size is, however, fairly weak. Moreover, it does not occur when spatial autocorrelation or niche position and breadth are taken into account, although this may partly be a consequence of the variation in body mass being smaller in species for which niche data are available (inter-quartile ranges are, respectively, 364 and 491 g for species with niche data and all species). It is possible that an explanation for the lack of strong support for the thermoregulatory load hypothesis is that adaptations enable smaller animals to survive in environments that would otherwise be thermally challenging to them, but such adaptations seem unlikely to remove the negative relationship between body size and rate of heat loss.

Regardless of the method used to measure the strength of ISERs, there is almost invariably a positive correlation between their strength and population size (tables 1 and 2; electronic supplementary material, figure 3). Many studies, including those of British breeding birds, that have grouped species into either abundance or range size classes have also found that the commonest species exhibit the strongest species–energy relationships (Jetz & Rahbek 2002; Ruggiero & Kitzberger 2004; Vázquez & Gaston 2004; Evans et al. 2005b,c). The MIH describes the relationship between extinction risk and local population size, which may appear to contrast with our use of national population estimates as a measure of abundance. However, in British breeding birds, local and national population sizes are strongly positively correlated with national range size and each other, and there is little evidence for marked simple spatial patterns in local population densities, such as lower densities at range edges (Gaston et al. 1997; Blackburn et al. 1999). Therefore, while we assume that national population sizes are a surrogate measure of local population size, there is strong evidence to support this assumption. Consequently, the finding that the commonest species exhibit the strongest species–energy relationships provides some evidence against the MIH, as this predicts the converse pattern.

Trophic level is rarely retained in the best-fitting models of the strength of ISERs. In a small number of models, however, predators exhibit significantly weaker ISERs than either herbivores or omnivores. The total number of avian predator species in 10 km squares also responds less strongly to environmental energy availability than the total numbers of herbivore or omnivore species (Evans et al. 2005c). This provides evidence against the trophic levels hypothesis that species–energy relationships arise because greater availability of environmental energy increases the length of the food chain, thus promoting the occurrence of predators (Oksanen et al. 1981; Fretwell 1987). Indeed, it has been suggested that energy availability will only limit food chain length in environments where it is available at much lower levels than in the UK (Post 2002).

The strength of ISERs is almost invariably negatively correlated with logarithmically transformed niche position (table 1), indicating that species that use specialized resources respond to energy more weakly than those that use more widespread resources (specialists have high niche positions). While the mechanisms behind this pattern are unclear, it provides evidence against the niche-position specialization hypothesis that greater environmental energy increases the availability of scarce resources, enabling specialized species to occur and thus promoting species richness (Abrams 1995; Evans et al. 2005a).

In four of the six types of multiple regression models, the strength of ISERs is negatively correlated with niche breadth (table 1; see also the electronic supplementary material, figure 3). These results fit the suggestion that high levels of environmental energy enable species to specialize on a few resource types, leading to reduced competition, greater co-existence and higher species richness (Abrams 1995; Evans et al. 2005a). Such an effect is not observed when species are classed into two groups with above- and below-average niche breadths (Evans et al. 2005c), perhaps because in the approach adopted here other factors such as population size are taken into account or because treating niche breadth as a continuous variable provides more statistical power.

The canonical correspondence analysis, generating the niche breadth measure, used temperature as one of the 34 base variables, so niche breadth may contain a small proportion of information on the range of temperature values within a species range, which, in turn, may be related to the strength of the ISERs. To assess the extent to which such non-independence may influence our results, we measured the strength of the relationship between niche breadth and the range of temperature and NDVI values within each species' distribution (respectively, r2=13.1 and 3.0%), and between temperature range and the strength of the single species–temperature relationship (deviance r2=7.8%; slope r2=15.7%). Such correlations are much weaker than the univariate relationships between the strength of ISERs and niche breadth (electronic supplementary material, table 3). Therefore, while non-independence may partly contribute to our finding that species with smaller niches tend to exhibit stronger ISERs, it cannot fully explain it.

We are not aware of any other evidence supporting the suggestion that species–energy relationships arise through availability of environmental energy promoting the occurrence of niche breadth specialists. Despite this, it is often suggested that narrower niches in the tropics (where availability of environmental energy is generally high) can partly explain latitudinal gradients in species richness. MacArthur (1955, 1972) suggested that greater environmental stability and lower seasonality in the tropics would allow population size to be more stable in the tropics, enabling narrower niches to evolve. In a detailed meta-analysis of this suggestion, Vázquez & Stevens (2004) found that there was no general trend for reduced niche breadth in the tropics, but statistical power of the meta-analysis was low and most studies reported trends in the expected direction. Vázquez & Stevens (2004) concluded that niche breadths are narrower in the tropics, but that this is a consequence of gradients in species richness rather than its cause. They suggested that high levels of environmental energy promote species richness and, assuming that species interactions are structured in an asymmetrically specialized nested fashion, that this promotes greater specialization. We currently lack the data to assess the causal direction of the correlation that we observe between the strength of the ISERs and reductions in niche breadth. Regardless of its cause, if our results generalize, then species that use a narrow range of resources are more likely to occur in high-energy areas and contribute to positive relationships between total species richness and environmental energy availability.

Acknowledgments

K.L.E. was supported by The Leverhulme Trust, and S.F.J. by the UK Population Biology Network (UKPopNet) funded by the Natural Environment Research Council (Agreement R8-H12-01) and English Nature. We thank the thousands of volunteers who collected the ornithological data, which was supplied by the British Trust for Ornithology. T. M. Blackburn helped with the PGLS analyses. Autologistic regression codes were supplied by F. He, who generously helped with the autologistic analyses and supplied the codes written and developed by H. Zhu.

Footnotes

References

View Abstract