International Journal of Ecology

Volume 2012 (2012), Article ID 947103, 12 pages

doi:10.1155/2012/947103

Research Article

Trait-Environment Relationships and Tiered Forward Model Selection in Linear Mixed Models

Tahira Jamil,¹ Wout Opdekamp,² Rudy van Diggelen,² and Cajo J. F. ter Braak¹

¹Biometris, Wageningen University and Research Centre, P.O. Box 100, 6700 AC Wageningen, The Netherlands

²Ecosystem Management Research Group, Department of Biology, University of Antwerp, Universiteitsplein 1c, 2610 Antwerpen, Belgium

Received 23 December 2011; Revised 1 April 2012; Accepted 2 April 2012

Academic Editor: Jean Clobert

Copyright © 2012 Tahira Jamil et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

To understand patterns of variation in species biomass in terms of species traits and environmental variables a one-to-one approach might not be sufficient, and a multitrait multienvironment approach will be necessary. A multitrait multienvironment approach is proposed, based on a mixed model for species biomass. In the model, environmental variables are species-dependent random terms, whereas traits are fixed terms, and trait-environment relationships are fixed interaction terms. In this approach, identifying the important trait-environment relationship becomes a model selection problem. Because of the mix of fixed and random terms, we propose a novel tiered forward selection approach for this. In the first tier, the random factors are selected; in the second, the fixed effects; in the final tier, nonsignificant terms are removed using a modified Akaike information criterion. We complement this tiered selection with an alternative selection method, namely, type II maximum likelihood. A mesocosm experiment on early community assembly in wetlands with three two-level environmental factors is analyzed by the new approach. The results are compared with the fourth corner problem and the linear trait-environment method. Traits related to germination and seedling establishment are selected as being most important in the community assembly in these wetland mesocosms.

1. Introduction

Understanding the processes that drive community assembly has been and still is a major challenge in community ecology [1, 2]. Many studies have already shown the importance of environmental factors in controlling and shaping species composition [2, 3].

The use of species-traits instead of species identity in community ecology research has many advantages as the latter reflect a species adaptation to its environment [4]. Hence, a trait-based approach not only allows a comparison of the same process in different vegetation types (e.g., [5–7] but also gives insight into the mechanisms responsible for such patterns [8] and allows predictions about possible future changes.

In the last decade, plant trait data have become more easily available, especially in Western Europe (LEDA, BIOLFLOR, etc.). This has further strengthened the growing interest of ecologists to study the responses of plant functional traits to environmental conditions [9, 10].

Despite this increasing interest, our knowledge of plant community assembly is still hampered as the quantification of the effect of plant traits on community assembly stays a real statistical challenge [11]. Most studies on trait-environment relationship, especially model-based ones, are limited to single species. Knowledge about the effect of traits on plant community assembly stays limited. Empirical evidence is limited as most of the studies are observational and correlative [12, 13]. Therefore, there is an immense need for randomized multispecies experiments that study trait-environmental links and for statistical methods that can link the experimental environmental factors to traits in such multispecies studies. This paper proposes such methods and applies them to a factorial three-year mesocosm study of plant communities with three environmental factors, each on two levels. The experimental measurements are the biomass per species in each of the three years.

The linkage between the traits and the environment is expressed differently in different statistical models. It is a Pearson’s correlation in the fourth corner problem [11] and an interaction term in the mixed model approach, as we will show in this paper. Environment-trait relationships are usually very complex due to the high number of interacting environmental variables and traits. To understand patterns of variation in species density, a multivariate approach will be necessary as a one-to-one approach might not be sufficient. This way, selecting important trait-environment interactions becomes a model selection problem in mixed models.

Mixed models are extremely flexible and form a computationally attractive tool to model complex and large datasets. Their potential applications in ecology are numerous. The resulting flexibility and model complexity make model selection even more vital [14]. In mixed models, model selection not only includes selecting the best mean structure but also the most optimal variance-covariance structure [15, 16]. Despite the fact that mixed models have been available for a few decades, there is surprisingly little literature available concerning model specification, that is, “which set of candidate models should be considered?” “How to select a model?” and “What is the best model to use?” These are the critical questions in making valid inference from data. This also includes the variable selection problem in mixed model analysis.

One of the goals of model selection is a trade-off between model complexity and accuracy. Depending on the modeling objective, different procedures to select an optimal model subject to a particular criterion are available. However, it is important to adopt a model selection procedure that reflects the ultimate objective of the modeling process [17]. Model selection is often done through sequential testing either step-up (forward) or step-down (backward) regression methods [18]. For simple problems, the outcomes of model selection using these two approaches might happen to be similar; however, in more complex situations, with many candidate models, the results of the two approaches may be quite different. Yet, selecting the best model from all possible models with different fixed and random effect factors is computationally forbidding as the number of models grows exponentially with the number of factors.

This paper develops a novel model selection method called “tiered forward selection.” The method uses a modified Akaike information criterion. In the first tier, the random factors are selected; in the second, the fixed effects; in the final tier, nonsignificant terms are removed. In our case study, the random factors are the environmental factors, while the fixed effects are related to traits and trait-environment interactions. We complement this tiered selection with an alternative one-shot method, namely, type II maximum likelihood (type II ML) Jamil and ter Braak [19].

The paper is structured as follows. After a short description of the mesocosm example data and association data screening, the linear mixed model, the tiered forward selection, and Type-II ML are presented. Next, we describe two simple existing methods for detecting trait-environment relationships that are not based on mixed models. These methods, the fourth corner method and the linear trait-environment method (LTE), use permutation for determining statistical significance. After presenting the results, we discuss statistical issues and shortly interpret the results in biological terms.

2. Material and Methods

Data

Data are from an outdoor mesocosm experiment investigating early community assembly from a pool of floodplain species covering a wetness gradient. The mesocosm experiment was constructed with 80 PVC-containers (W × L × H :  111 × 91 × 61 cm), comprising of a 25 cm drainage layer (Argex clay aggregate, Argex NV, Zwijndrecht, Belgium) covered with a geotextile boundary and then 25 cm soil. The soil consisted of a homogenised mixture of fen peat and alluvial clay (3 : 1 v/v).

Experimental treatments were (i) waterlogging (with and without waterlogging, w/nw), (ii) canopy presence (with and without canopy, c/nc) at initial inoculation, and (iii) mowing (with and without summer mowing, m/nm) for a full-factorial design (2 × 2 × 2) with 10 replicates.

We selected 34 plant species that occur frequently in floodplains of temperate Europe and represent the entire moisture gradient from mesic meadows to reed beds and swamps (see Table 5 for species list). Seeds were sown in all mesocosms at a rate of approximately 1000 seeds species⁻¹ m⁻² on May 15, 2006. The seeding rate was similar to seed production on wet meadows [20]. The canopy (c/nc) consisted of Poa pratensis, Lolium perenne, and Alopecurus pratensis. These grasses were used as matrix vegetation in the canopy treatment but were also part of the species pool introduced to all containers. They were pregrown for 6 weeks and cut at 10 cm height when the other seeds were added. This grass species mixture ensured the persistence of a grass canopy in both dry and wet conditions. In other mesocosms, soil was kept bare until the sowing date. Waterlogging (w/nw) was maintained at 5 cm below soil surface. Nonwaterlogged soils were allowed to drain freely. Mowing (m/nm) involved annual mowing of vegetation to 2 cm in June-July.

A more detailed description of the experimental setup can be found in Kotowski et al. [3]. Aboveground plant biomass was harvested in August of 2006, 2007, and 2008. The harvest was sorted to species level with the exception of Poa, Lolium, and Alopecurus, which were grouped together. Dry mass was determined after 72 h of drying at 70°C.

Traits were either measured (for details see [3]) or extracted from several plant trait databases:  BiolFlor [26], CloPla [22], and LEDA [23]. Where a database contained several entries per species-trait combination, the following aggregation scheme was used. For categorical traits, we discerned two cases: entries to different trait-attributes are likely (e.g., a species can have several clonal growth organs) and species have only one possible entry per trait (traits: HEM, STA, DUR, and SEX). In the first case (traits: AB0, AB2, GR2, BE0, BE1, GRS, BES, LAT2, and LAT3) the number of entries per trait-attribute was expressed relatively. For the trait “clonal growth organs” this would for example, look like “no clonal growth organs” , and “aboveground runners” , “aboveground plant parts”/“belowground runners” . In case more than one entry per species-trait combination is found in a database, these values are weighted by the number of entries per trait-attribute and expressed relatively to the total number of entries for that trait. For numerical traits (traits: CH, LDMC, and SLA), the median value of all possible entries was used. A description of the different traits can be found in Table 1.

Table 1: Traits used for analysis with code and description and the data source (Kotowski et al. [3], BiolFlor [21], CloPla [22], and LEDA [23]).

Germination traits
Source

Z1
SW
Seed weight (mg)
Kotowski et al. [3]

Z2
GP
Total germination percentage in full light (%)
Kotowski et al. [3]

Z3
T50
Time of 50% germination in light (days)
Kotowski et al. [3]

Z4
WGR
Wet germination ratio (germination in wet mesocosms/germination in dry mesocosms)
Kotowski et al. [3]

Z5
DGR
Dark germination ratio (dark germination/light germination)
Kotowski et al. [3]

Seeding traits

Z6
H7
Average height of seedlings at the 7th day from germination (mm)
Kotowski et al. [3]

Z7
LWR7
Mean leaf weight ratio at the 7th day ( leaves and cotyledons)
Kotowski et al. [3]

Z8
LAR7
Mean leaf area ratio,that is, mean quotient of the total leaf area per plant and the total weight per plant at the 7th day (mm² mg⁻¹)
Kotowski et al. [3]

Z9
AGR
Mean actual growth rate of seedlings between the 7th and 22nd days (mg day⁻¹)
Kotowski et al. [3]

Z10
RGR
Mean relative growth rate of seedlings between the 7th and 22nd days
Kotowski et al. [3]

Z11
LWR
Mean leaf weight ratio; increase between the 7th and 22nd days
Kotowski et al. [3]

Z12
LA7
Mean leaf area of seedlings at the 7th day (mm²)
Kotowski et al. [3]

Adult traits

Z13
CH
Canopy height (maximum value) (m)
LEDA (Rothmaler only)

Z14
LDMC
Leaf dry matter content (mg g⁻¹)
LEDA

Z15
SLA
Specific leaf area (adult plants only) (mm² mg⁻¹)
LEDA

Z16
HEM
Leaf distribution along stem (hemirosettes)
Biolflor

Z17
STA
Flowering start (1: January–April; 2: May-June; 3: July–September)
Biolflor

Z18
DUR
Flowering duration (1: short (1-2 months); 2: medium (3-4 months); 3: long (> 4 months))
Biolflor

Z19
SEX
Reproductive type (sexual)
Biolflor

Z20
AB0
BudBank vertical distribution—aboveground: 0 (0 buds)
CloPla

Z21
AB2
BudBank vertical distribution—aboveground: 2 (>10 buds)
CloPla

Z22
GR2
BudBank vertical distribution—groundlevel: 2 (>10 buds)
CloPla

Z23
BE0
BudBank vertical distribution—belowground: 0 (0 buds)
CloPla

Z24
BE1
BudBank vertical distribution—belowground: 1 (1–10 buds)
CloPla

Z25
GRS
BudBank seasonality—groundlevel-seasonal
CloPla

Z26
BES
BudBank seasonality—belowground: seasonal
CloPla

Z27
LAT2
Lateral spread (2: 0.01–0.25 m year−1)
CloPla

Z28
LAT3
Lateral spread (3: >0.25 m year−1)
CloPla

Data Screening

Prior to analysis, trait data was screened for zero-variance predictors and for multicollinearity among predictors. Predictors with a single unique value (also known as “zero variance predictors”) and near zero-variance predictors (for details see [27]) can cause numerical problems and lead to misinterpretation. Both near zero-variance and zero-variance predictors were removed from the dataset. Predictor detection was performed using the caret package for R [28]. Multicollinearity among predictors is a problem in mixed models. From each pair of trait predictors with correlation greater than 0.80, one predictor was removed using the “find Correlation” function in the caret package in R [28]. This function removes the predictor that has highest mean pairwise correlation with the other predictors.

Trait predictors were checked for normality by making histograms. Predictors departing from normality (WGR, H7, AGR, and LA7) were log-transformed (Table 1). Once the final set of predictors was determined, predictors were centered and scaled using their mean and standard deviations [27]. Species with small average biomass (<0.1) were removed from the analysis, and 23 species remain. Species biomass was log-transformed (log()).

The three experimental factors are indicated by c, w, and m, and the levels of each are coded numerically as −1 (nc/nw/nm) and 1 (c/w/m) without loss of generality. In this coding, c, w, and m are orthogonal and also orthogonal to the interactions cw, cm, and wm.

2.1. Linear Mixed Models for Trait-Environment Relations

2.1.1. Single Trait-Environment Relationships

The model for , the species biomass for the jth species in the th site (mesocosm) is

where is a known environmental factor, is the site effect, and is the error term. This model is of interest for our case study as each of the factors in the mesocosm experiment was coded as if it were quantitative, with the levels of the factors coded as −1 and +1. See [29] in case X is a multilevel factor. We assume that the intercept and slope for the jth species depend on the known value of a particular trait [29]:

and also assume and . This is a random intercept and random slope model, where trait is a predictor for the random intercept and slope. Inserting , and in (1) gives

where represents the trait-environment interaction. The above model is fitted for each combination of trait and environmental variable. To test the trait-environment interaction (with null hypothesis: = 0), we fit the model without this term:

and then compare the two models by an analysis of variance resulting in a value for the likelihood ratio (LR) test of model with trait-environment interaction term against the null model without this term.

2.1.2. Multiple Trait-Environment Relationship

In community assembly, traits and environmental variables should not be considered in isolation as they influence and often coordinate each other. The development of a multivariate framework, in which multiple traits can be linked to multiple environmental variables is needed [30]. In this section, we show how to link environmental variables and species traits in a multivariate framework. We use a multitrait and multienvironment version of the mixed model to select the species traits, environments, and trait-environment interactions that significantly contribute to the species biomass distribution model. Equation (1) for a one-to-one model can readily be extended to cover multi-trait and multi-environmental variables:

where Greek letters are used for random terms. Environmental variables enter the model both as fixed terms and as random terms, and traits enter as fixed terms. In this formulation, the selection of the best traits and environmental factors is a model selection problem

Tiered Forward Models Selection

When the number of fixed effects and random effects is large, it is computationally very expensive and time-consuming to compute all possible candidate models [31] due to the presence of random terms and variance-covariance structure [17, 32]. Furthermore the number of candidate models increases exponentially with increased number of fixed effects and random effects [33]. Different protocols for model selection have been developed, in particular stepup (forward) [34] and top-down (backward) protocols [35]. Most stepwise functions take a start model and according to some criteria iteratively add or delete a predictor at each step, to get to the best parsimonious regression model. Backward selection starts from the model with all possible terms included and is only feasible if that model can be fitted. In our case, the number of environmental variables should be less than the number of sites and the number of traits should be less than the number of species. In our data, the former holds true, but the latter does not. Here, we develop an approach called tiered forward selection. The analysis was done for the three years separately and for the combined data of all three years.

Model Selection Criteria

Different Information criteria, such as AIC [36], AICC [37], CAIC [38], and BIC [39], can be used for model selection in linear mixed model [40]. Generally, these information criteria are a function of the likelihood for a given model and a penalty term based on the number of parameters in the model. The general form of information criterion (IC) is

where is the log likelihood derived from fitting the mixed model to the data using either ML or REML.

The use of these criteria is somewhat arbitrary, and no formal inference can be made based on these values. Comparison of the values of the criteria for a set of candidate models simply indicates if a superior model exist within the given candidate models [40]. For an extensive review and discussion on the theoretical aspects of model selection criteria and procedures see [41] and [17]. The most widely used information criterion is

We use the variant SigAIC defined as

where is the number of parameter estimates. The variant, SigAIC, which multiplies df by instead of by 2 [42], guarantees that the addition of a single parameter to a model will result in a lower SigAIC value if and only if that parameter is significant at the 5% level as judged by the LR test.

Phase Tier I: Selection of Environmental Variance Components

The start or null model is the model with crossed random effects for species and sites. In the R package lme4 [43], it can be represented as

start.model <- lmer( y ~1 + (1|sp)+(1|site), data),

where represents the vectorized response data, while sp and site indicate species and sites, respectively. REML is the default estimation method. In each consecutive step, the environmental predictor for which the species-dependent random effect term increases the log-likelihood most is added. This means that all models with one extra term have to be fitted each step, and the best term is retained for the next step. For example, for “w”:

lmer( y0 ~1+(1+w|sp)+(1|site),data).

In lme4, such model can be fitted as an update of the start model with the statement:

update(start.model,.~. + (1 + w|sp) + (1|site) - (1|sp), data).

To generate all models needed in this tier, we used the statement

update(start.model,as.formula (paste(". ~.+(1+",block[j],"|sp) + (1|site)-(1|sp)"))).

Here, block is a vector of the candidate predictors, that is, the environmental factors. All models, after fitting the predictors in the block, are arranged in the the order of the predictive criterion. The best candidate model is compared with the null model. If the best candidate model is statistically significant, it becomes the null model for the next step. This process is repeated until the increase in log-likelihood is no longer statistically significant as judged by LR test.

After this tier, the selected predictors are added as fixed effect. Thus, environmental variables in the model are now the component of both fixed terms and random terms.

Tier II: Selection of Fixed Trait Effects

In our case study, the start model for the next phase has a specification:

start.model<-lmer(y~c + w + cw + (1+c+w+cw|sp) + (1|site),method = “ML”, data).

Now traits and trait-environment interactions can be added as fixed effects to the model. It is important to note that REML is the default estimation method for mixed models. Generally REML estimates of variance components are preferred to ML estimates. However, in REML it is not legitimate to compare models with different sets of fixed effects as the contrast used to develop the restricted maximum likelihood depends on the fixed effect design matrix [44]. Therefore, the ML estimation method is used in this tier.

Given a starting model and a set (here block) of variables to evaluate, the starting model is updated by adding every single trait variable and trait-environment interaction. To evaluate the importance of a trait, the current model is updated with the statement:

update(start.model,.~.+ Z + c:Z).

A generic way to update the current model with any single trait is

update(start.model,as.formula(paste(". ~.+ ",block1[j],"+",block2[j]))).

Here, block1 is for trait main effects and block2 for the trait-environment interaction. Models fitted in this way are then ordered based on the chosen predictive criterion, here SigAIC, after which the best fitting model is retained for the next step. The procedure continues until the addition of new traits, and trait-environment interactions do not significantly improve the model.

Next, the same procedure is repeated for the three-way interaction, but keeping the marginal effect of a trait variable and its two-way interaction:

update(start.model,as.formula(paste(". ~. + ", block1[j],"+", block2[j], "+", block3[j],"+", block4[j]))).

Block1 is for trait main effects, block2 and block3 for two-way interactions, and block4 for three-way interaction. This structure ensures the marginality principle [21] which entails a model can include an interaction term (high-order term) only if it includes the main effects (and all lower-order terms) that compose the interaction. The condition requires that a mixed model with an interaction, say X : Z, must also include the main effects X and Z. In general, we neither test nor interpret the main effects of explanatory variables that are also included in an interaction. The procedure continues until addition of new predictors does not significantly improve the model.

Tier III: Removal of Nonsignificant Interaction Terms

In this tier nonsignificant interaction terms are sequentially removed. An example statement in this tier is

update(start.model, as.formula(paste(". ~. - ", block1[j]))).

The final model was obtained by sequential removal of nonsignificant interaction terms. Now we refit the final model using REML estimation, this is needed to obtain unbiased estimates of the covariance parameters. As the ML estimation leads to biased covariance parameter estimates. This can result in smaller estimated standard errors for the estimates of fixed effects in the model. The REML estimate for the fixed effects is not identical to its ML version and differs more from ML as the number of fixed effects in the model increases.

2.2. Model Selection by Type-II Maximum Likelihood

When a system is described by a statistical model, model complexity leads to a very large computing time and poor estimation, especially if the number of predictors is large relative to data size. As an alternative to and improvement over stepwise methods, shrinkage methods have been proposed. One of these is the relevance vector machine (RVM) which has gained popularity within the Bayesian framework [45–47]. RVM introduces a Gaussian prior to the regression coefficients, with one variance component for each predictor. The variance component or hyperparameter controls the degree of sparseness. These parameters are commonly adjusted by cross validation. In Bayesian framework, the hyperparameters are estimated by using empirical Bayes or, equivalently, Type-II maximum likelihood [48].

The Type-II maximum likelihood shrinks the coefficients to zero and readily sets some of the coefficients to zero, namely, if their variance component is estimated to be zero Jamil and ter Braak [19]. These zero variance coefficients are equivalent to pruning the corresponding predictors from the model. Hence, this method readily helps in pruning predictors from the model and does variable selection and model estimation. A prototype statement in lme4 in linear regression is Jamil and ter Braak [19]

lmer(y ~ (0 + x1 | v) + (0 + x2 | v), data, REML = FALSE),

where is the response and and predictors, is an all ones N-vector, and train is a data frame containing these vectors. In our multi-trait-multi-environment context a prototype statement is

lmer(y ~ (0 + z1:x1 | v)+(0 + z1: x2|v)+(0 + z2: x1 | v)+(0 +z2: x2|v) +(1+x1+x2|sp)+ (1|site), data, REML = FALSE),

where and indicate two environmental factors, and two traits, and z1 : x1 and related term indicate the vector that is the product of the corresponding trait and environmental factor.

2.3. Fourth Corner Method

The fourth corner method, developed by Legendre et al. [49] and extended in Dray and Legendre [11], is the oldest integrated trait-environment method. The species data in this method must be nonnegative and can thus be presence-absence data, abundance data, or species biomass. The fourth corner statistic measures the link between trait and environment via the species data table with a weighted Pearson’s correlation coefficient between trait and environmental variable, each vectorized as in the mixed model approach [26]. The weights are the species data (presence-absence, abundance, or biomass). The significance of the trait-environment relationship is tested by a permutation test. Dray and Legendre [11] offer different permutation scenarios. We used the combined approach implemented in the combine. 4th corner function in the ade4 package in R. It combines the values of two fourth corner models, namely, Model 2 (site permutation) and model 4 (species permutation), as proposed by Dray and Legendre [11] by taking their maximum [50]. This method controls the type I error [51].

2.4. The Linear Trait-Environment Method (LTE)

The linear trait-environment (LTE) starts with a two-step analysis Cormont et al. [51]. In the first step, regressions per species biomass to each environmental variable give a species-specific regression coefficient:

where and are intercept and slope of kth species. In the second step, these regression coefficients are correlated to each trait:

where and are the intercept and slope for trait Z, and is a species-specific error term with mean 0. LTE integrates both steps in a single model. LTE achieves this integration based on a linear model with main effects for the trait and environmental variable and their interaction. The interaction between a trait and an environmental variable in this model captures the trait-environment relationship, in particular the trait-dependent effect of environment on species biomass. The significance of this interaction is tested by a permutation test with the same permutation strategy as in the fourth corner problem [51].

Summary of Data Analysis Strategy

We applied the different methods that link functional traits to the experimental environmental factors to data on 23 floodplain species in a factorial mesocosm experiment with three two-level factors (c, w, and m for canopy, waterlogging, and mowing). We estimated the one-to-one interactions for each trait-environment combination by mixed models, the fourth corner method [11, 49] and the linear trait-environment method [51]. Then, we applied the tiered forward selection method to select the traits and environmental variables that significantly contribute to the explanation of species biomass. The environmental variables were c, w, m, and their first-order interactions. Furthermore, we performed Type-II maximum likelihood analysis (Type-II ML) to select trait-environment interactions. The analysis has been carried out on the three-year mean log-transformed biomass, unless stated explicitly otherwise.

3. Results

In this study, we explored different methods that link the species traits to the environmental factors in the mesocosm experiment. Table 2 summarizes the results of the one-to-one analyses by mixed models, the fourth corner method and the linear trait-environment method, and the multitrait-to-multienvironment analyses by mixed models and Type-II maximum likelihood. Only significant trait-environment interactions obtained from these analyses are reported together with the sign of the relationship. The factor mowing is excluded from Table 2 as its variance component in the mixed model approach was very small, smaller than that of the interaction “cw” between the other two factors, canopy “c” and waterlogging “w.” Mowing thus has a similar effect on all species. The sign of coefficient for trait-environment interactions is almost consistent between different methods. Mixed model search for trait-environment combinations (interactions) predicts the response of species to the environmental variables and traits variables.

Table 2: Results of analysis from different methods using species biomass, species traits and environmental variables. The sign (+/−) represents the positive or negative significant Pearson’s correlation/coefficients between the environmental variable and the species trait. For ease of interpretation, only significant relationships are shown.

One-to-one mixed model
Linear trait-environment
Fourth corner method
Multivariate mixed model
Type-II ML

c
w
cw
c
w
cw
c
w
cw
c
w
cw
c
w
cw

SW

−

−
+

−
+
+
−

+
−
+

T50

+
−

−

WGR

+

+
−

+

+

DGR
−
−

−
+

−

+

H7
−
−
+
−
−
+
−
−
+
−

−

LWR7

+

+

−
+

AGR
−

−
−
+

+
−

−

LA7

−

−
+

SLA

−

HEM

+

+

STA

+

+
−

DUR

−
+

−
−

BE0
+

+

+

+

+

BE1

−

GRS

+
−

+

+
−

+
−

BES

−

Type-II ML also identifies the trait-environment interactions which are also identified by other methods, but the selection is optimistic. It encountered too many trait-environment interactions, and many of them were not common to other methods as Type-II ML is very tolerant in allowing predictors to stay in the model Jamil and ter Braak [19]. In case of correlated predictors, Type-II ML prunes the weaker of the two.

H7 and SW appeared important traits that are consistently significant with canopy or/and waterlogging for all methods (Table 2). Plants with small seedlings performed best when waterlogging or canopy acted as a sole stress factor. When both factors interacted, large seedlings with heavy seeds had an advantage. Although less consistent across methods, DGR and AGR seem important as well. Species with a later flowering onset (STA) and a seasonal budbank (GRS) were significantly linked to waterlogging in all models, except for the multivariate mixed model. Type-II ML gave many more trait-environment interactions than other methods (Tables 2 and 4).The coefficients estimated by Type-II ML are plotted in Figure 1. The length of bar is proportional to magnitude of coefficient, and direction of bar indicates the sign of coefficient.

Figure 1: Estimated effects of trait-environment interaction using Type-II ML.

Table 3 summarizes the results of the tiered forward model selection in mixed models using SigAIC. The analysis was done for each year and by combining the data for all three years. All significant traits, except BE0, are germination and seedling establishment traits. Canopy captures more interactions with traits in first and second years of germination and establishment as compared to waterlogging.

Table 3: Results of the multitrait and multienvironment analysis using tiered forward model selection in mixed models. The following treatments c: canopy, w: water logging, and cw: interaction c × w were always in the model as a fixed effect and random effects.

Year 1
Year 2
Year 3
All

(Intercept)
−2.66
−2.24
−1.67
−2.20

c
−3.28
−1.97
−1.27
−2.17

w
0.45NS
−0.06NS
0.28NS
0.22

cw
−0.78
−1.02
−0.45NS
−0.75

SW

0.58
0.34
0.11

c : SW

0.86

w : SW

−0.59
−1.24
−0.72

WGR

−0.11NS
−0.12

w : WGR

0.88
0.46

H7
0.83
0.56
0.63
0.77

c : H7
−1.53
−1.49
−1.16
−1.47

w : H7
−0.75
−0.90

cw : H7
0.60
0.63

LWR7

−0.51
0.08NS
−0.28

w : LWR7

0.66
1.09
0.67

AGR
0.67

0.47

c : AGR
−1.26

−1.38

SLA

−0.45
−0.13

w : SLA

−0.44
−0.29

HEM

0.36

BE0

−0.26NS
0.38
0.08

w : BE0

0.87
0.56

c : BE0

0.99

NS is for nonsignificant.

Table 4: Trait-environment interactions from Type-II that were not common to other methods.

Type-II ML

c
w
cw

GP

+

RGR
−

LWR

−

CH

+

LDMC
−

+

AB0

+

AB2

+

GR2

+
+

LAT2
−

LAT3
+
−

Table 5: Species selected for the experiments. F: Ellenberg moisture indicator [24]. SM: mass of 1000 seeds (amount introduced to each mesocosm); names follow Tutin et al. 2001.

Mesic habitat species
F
SM (g)
Wetland species
F
SM (g)

Achillea millefolium
4
0.17
Carex ovalis
7
0.49

Daucus carota
4
0.87
Cirsium oleraceum
7
1.90

Galium verum
4
0.54
Eupatorium cannabinum
7
0.22

Leucanthemum vulgare
4
0.43
Geum rivale
7
1.21

Angelica sylvestris
5
2.46
Cirsium palustre
8
1.11

Centaurea jacea
5
2.33
Epilobium hirsutum
8
0.11

Heracleum sphondylium
5
7.35
Filipendula ulmaria
8
0.60

Leontodon autumnalis
5
0.80
Lythrum salicaria
8
0.06

Lolium perenne
5
0.90
Myosotis palustris
8
0.28

Plantago lanceolata
5
2.63
Phalaris arundinacea
8
1.02

Poa pratensis
5
0.31
Senecio aquaticus
8
0.24

Prunella vulgaris
5
0.86
Valeriana officinalis
8
0.72

Rumex acetosa
5
0.81
Carex acuta
9
0.90

Alopecurus pratensis
6
0.71
Lychnis flos-cuculi
9
0.16

Ranunculus acris
6
1.79
Lycopus europaeus
9
0.27

Urtica dioica
6
0.16
Peucedanum palustre
9
1.35

Veronica longifolia
6
0.07
Typha latifolia
10
0.04

See [24, 25].

We used three methods to relate each single trait to each single environmental factor. Two of the methods use an explicit linear model for the data (the linear mixed model and LTE), whereas the third method uses the data are weights. Also, two of the methods use statistical significance testing by a Monte Carlo permutation strategy (LTE and the fourth corner method), whereas the mixed model uses LR testing. Overall, the results produced by the different methods that were investigated are concordant regarding the key trait-environment relationships.

What effects do these theoretical similarities and dissimilarities have on the results? The numbers of one-to-one relationships found (Table 2) were 17 in the mixed model, 11 in LTE, and 16 in the fourth corner method. The mixed model disagreed in its significance judgment in 14 cases (out of the 48 shown in Table 2) with LTE and in 15 cases with the fourth corner method. LTE and the fourth corner disagreed in 13 cases. The methods, thus, are about equally dissimilar among one another. The multi-trait multi-environment mixed model yielded 8 significant interactions, so fewer than the one-to-one methods. This is to be expected if one trait can replace another because of their mutual correlation. In our data, traits show some, but no high correlation. Apparently the other aspect of a multivariate model, namely, that is it can reduce the error variance and thereby can find more significant effects, is less important.

4. Discussion

4.1. Statistical

For analysis of trait-environment relationships, here we presented and compared multivariate methods that combine three matrices (i) species biomass × sites/plots, (ii) sites × environmental variables, and (iii) species × traits.

The relationship between species traits and the environment is generally assessed indirectly using a two-step approach. Species biomass is first linked to environmental conditions, and species responses to environmental variation are then related to the biological and/or physiological traits of species. In such analyses, the relationship between the environment and the species traits is thus assessed indirectly. The fourth corner method [11] calculates a weighted Pearson’s correlation between the trait and the environmental variable by using all species-site combinations as cases, the measure of abundance as a weight and by assigning to each case the trait and the environmental value of the combination. This generates two matrices: species-trait, and species-environment, weighted by abundances and then analyzed simultaneously. It sacrifices some information at the species level, as zero abundance implies zero weight. The later approach may also be called direct functional analysis. Direct assessment of the traits-environment relationships requires a statistical method that takes into account simultaneously the information stored in three matrices to link species traits and environmental conditions through species responses.

This paper developed a novel multi-trait and multi-environmental variable model selection method called tiered forward selection. Here in the first tier, the random factors are selected; in the second, the fixed effects are selected; in the final tier, nonsignificant terms are removed based on a predictive modified Akaike information criterion. Here, random factors are the environmental variables, while the fixed effects are related to traits and trait-environment interactions. This is a direct functional approach which retains all the available information and provides guidelines for choosing key environmental variables, species trait, and also trait-environment interactions. Further, we compared the performance of mixed model with the fourth corner method, the linear trait-environment method (LTE) and the one-shot method, namely, Type-II maximum likelihood (Type-II ML).

Mixed model approach can be applied to assess how single species respond to environmental gradient and also to understand patterns of variation in species biomass. Fourth corner method provides trait environment relationships at the species level, and Fourth corner testing procedure [11] evaluates the statistical significance of the trait environment relationships. The role of the species data is thus rather different from that in mixed model and in LTE. In particular, absences or zeroes, therefore, do not carry any information in the analysis. The method works well for data stemming from unimodal response models [11].

The linear trait-environment (LTE) method Cormont et al. [52] was developed as an alternative to the fourth corner method to account for negative species data values. The method is linear and as such is more closely related to the liner mixed model than the fourth corner method. However, just like the fourth corner method, LTE has no variance components. The statistical analysis thus proceeds by advanced permutation testing. LTE estimates the parameters by least square. In contrast to LTE, linear mixed model (LMM) estimates the fixed effects by generalized least square estimation and their significance can be test by parametric bootstrap. The multi-trait multi environment method resulted in less number of significant trait-environment interactions. The fourth corner method tests the significance of trait-environment relationship by a permutation test. Dray and Legendre [11] offer different permutation scenarios but do not faithfully controlled the type I error ter Braak et al. [53]. The interaction between a trait and an environmental variable in LTE captures the trait-environment relationship. The significance of this interaction is tested by a Monte Carlo permutation test as in the fourth corner problem. The LTE differs from the fourth corner method by using multivariate linear regression [51].

In this study, we considered trait-environment relationships. In hindsight, perhaps the most important reason that it is difficult to quantify this relationship is that traits are measured on species and environmental variables on sites. So how can these be related? They can only be related via the sites × species data. So, there is a third entity involved, which complicates the issue. Statisticians are often keen to distinguish interactions from correlations. Two variables are correlated when a change in one variable is likely to be associated with a change in the other. By contrast, interaction involves a third variable and considers the effects of the variables on this third variable. Two variables are said to interact when the one variable modifies the effect of the other on the third variable. In terms of regression modeling, the third variable is the response variable, the other two are predictor variables, and the interaction might be represented by the product of the predictor variables.

Is the trait-environment relationship now best expressed as a correlation as in the fourth corner problem (and LTE) or as an interaction as in the LMM model? The fourth corner problem is able to express the relationship as a correlation by taking the individual organism as the statistical unit: the cases are the individual organisms. This trick gives the third variable (the sites × species data) another role; the elements become weights. This is the logical approach when individuals are (randomly) sampled rather than sites. However, in the practice of much ecological research, primarily sites are sampled and then individuals within sites. The sampling process is thus hierarchical, and hierarchical statistical models are thus a natural way to model it. We took this hierarchical approach. By giving the third variable the role of response variable, the trait-environment relationship becomes naturally an interaction. In contrast with the fourth corner problem, this approach does not ignore the information that some species are absent at a site. The advantage of the mixed is that it has the potential of predictive use: which species from a species pool are expected to occur under specified environmental conditions when we only know the trait values of the species in the pool. Of course, at the current stage, we are still ignoring any competition and successional processes that must also be important in community assembly, but that does not necessarily invalidate the prediction. It makes it less precise.

4.2. Biological Interpretation

A major objective of this study was to compare the effects of competition from canopy and waterlogging on assembly processes in a floodplain and how plant functional traits are related to the successful establishment of species. Both wetness (waterlogging) and light availability during initial community assembly (canopy presence) affected species germination in our experiment, but only few species were directly eliminated at this stage (one under oxic and three under dark conditions) [3]. Clearly, canopy presence was a much stronger filtering factor. Especially in the first year it almost totally disabled establishment, which is probably due to high light attenuation [3]. However, more severe root competition could possibly also play a role. Waterlogging is a major constraint on growth and establishment of plant species in wetlands [54]. Only species that are adapted to this environment can occur and thrive. It is relevant to understand which specific combination of traits makes the germination and establishment of plants resilient to waterlogging and which species traits showed the greatest tolerance to waterlogged conditions. The knowledge of which traits determine species seedling and germination might help for conservation planning [55].

Species able to establish successfully within the grass canopy showed high seed weight, combined with a small individual size and a relatively low actual growth rate. These traits are related to a stress-tolerant strategy, allowing plants to minimize resource requirements and survive in suboptimal conditions. Apparently even the largest seedlings had difficulties in reaching layers with sufficient light. This confirms that large seed size may contribute to seedling establishment in shade through various mechanisms [56]. Large seeds with a large nutrient stock are often thought to be advantageous in dense canopies as seedlings should possess enough resources to reach layers with higher light availability [57]. However, a significant fraction of resources in a large seed may remain in storage instead of being used for immediate seedling development [58, 59], a strategy characteristic for stress-tolerant species. When germinating on bare ground, seedlings of different species compete with each other for light. Yet, the intensity of this competition is much higher in the nonwaterlogged treatment. Traits responsible for rapid establishment and outcompeting neighbours appear more important here than those responsible for shade tolerance [60, 61]. A combination of fast growth and large-sized seedlings is prerequisite for success under dry conditions without imposed light stress. In waterlogged soils, specific leaf area (SLA) decreased as waterlogging induces an increased allocation to roots [54]. The ability to germinate in wet conditions is a main determinant of community assembly. This is in accordance with the habitat filter theory which states that the number of species in the local species pool is reduced by habitat constraints.

Hence, all traits (except one) that were selected by the tiered forward model selection describe germination and seedling establishment. This stresses the importance of both two stages as major bottle necks for species recruitment [52, 62] and how they may largely determine patterns of biological diversity [57, 63]. Moreover, because of their small stature, seedlings can be subject to a totally different light regime and soil resource availability than adult plants, even in the same site.

In conclusion, we have demonstrated different methods that link environmental factors (e.g., waterlogging and canopy) to species traits during early assembly process in a wetland mesocosm. Our results clearly stress how the choice of a particular statistical method to analyse the trait-environment link will have consequences for the ecological interpretation of this link. Of the studied methods, the multitrait multienvironmental mixed model is clearly best suited for predictive usage.

References

J. M. Diamond, “Assembly of species communities,” in Massachusetts, Harvard University Press, Cambridge, Mass, USA, 1975.

E. Weiher and P. A. Keddy, “The assembly of experimental wetland plant communities,” Oikos, vol. 73, no. 3, pp. 323–335, 1995.

W. Kotowski, O. Beauchard, W. Opdekamp, P. Meire, and R. van Diggelen, “Waterlogging and canopy interact to control species recruitment in floodplains,” Functional Ecology, vol. 24, no. 4, pp. 918–926, 2010.

S. Menezes, D. J. Baird, and A. Soares, “Beyond taxonomy: a review of macroinvertebrate trait-based community descriptors as tools for freshwater biomonitoring,” Journal of Applied Ecology, vol. 47, no. 4, pp. 711–719, 2010.

S. Díaz, I. Noy-Meir, and M. Cabido, “Can grazing response of herbaceous plants be predicted from simple vegetative traits?” Journal of Applied Ecology, vol. 38, no. 3, pp. 497–508, 2001.

S. Lavorel and E. Garnier, “Predicting changes in community composition and ecosystem functioning from plant traits: revisiting the Holy Grail,” Functional Ecology, vol. 16, no. 5, pp. 545–556, 2002.

S. Kahmen and P. Poschlod, “Effects of grassland management on plant functional trait composition,” Agriculture, Ecosystems & Environment, vol. 128, no. 3, pp. 137–145, 2008.

S. Kahmen and P. Poschlod, “Plant functional trait responses to grassland succession over 25 years,” Journal of Vegetation Science, vol. 15, no. 1, pp. 21–32, 2004.

E. Weiher, A. V. D. Werf, K. Thompson, M. Roderick, E. Garnier, and O. Eriksson, “Challenging theophrastus: a common core list of plant traits for functional ecology,” Journal of Vegetation Science, vol. 10, no. 5, pp. 609–620, 1999.

C. Violle, M.-L. Navas, D. Vile et al., “Let the concept of trait be functional!,” Oikos, vol. 116, no. 5, pp. 882–892, 2007.

S. Dray and P. Legendre, “Testing the species traits environment relationships: the fourth-corner problem revisited,” Ecology, vol. 89, no. 12, pp. 3400–3412, 2008.

B. J. McGill, B. J. Enquist, E. Weiher, and M. Westoby, “Rebuilding community ecology from functional traits,” Trends in Ecology and Evolution, vol. 21, no. 4, pp. 178–185, 2006.

D. Vile, B. Shipley, and E. Garnier, “Ecosystem productivity can be predicted from potential relative growth rate and species abundance,” Ecology Letters, vol. 9, no. 9, pp. 1061–1067, 2006.

S. Greven and T. Kneib, “On the behaviour of marginal and conditional AIC in linear mixed models,” Biometrika, vol. 97, no. 4, pp. 773–789, 2010.

R. Wolfinger, “Covariance structure selection in general mixed models,” Communications in Statistics, vol. 22, no. 4, pp. 1079–1106, 1993.

R. D. Wolfinger, “Heterogeneous variance-covariance structures for repeated measures,” Journal of Agricultural, Biological, and Environmental Statistics, vol. 1, no. 2, pp. 205–230, 1996.

J. A. Hoeting, R. A. Davis, A. A. Merton, and S. E. Thompson, “Model selection for geostatistical models,” Ecological Applications, vol. 16, no. 1, pp. 87–98, 2006.

W. D. Hosmer and S. Lemeshow, Applied Logistic Regression, Wiley, New York, NY, USA, 2nd edition, 2000.

T. Jamil and C. J. F. ter Braak, “Selection properties of type II maximum likelihood (empirical Bayes) in linear models with individual variance components for predictors,” Pattern Recognition Letters, vol. 33, no. 9, pp. 1205–1212, 2012.

A. Klimkowska, R. van Diggelen, S. den Held, R. Brienen, S. Verbeek, and K. Vegelin, “Seed production in fens and fen meadows along a disturbance gradient,” Applied Vegetation Science, vol. 12, no. 3, pp. 304–315, 2009.

J. A. Nelder, “Letter to the editors: what is the mixed-models controversy?” International Statistical Review, vol. 76, no. 1, pp. 134–135, 2008.

J. Klimešová and L. Klimeš, “CLO-PLA3. a database of clonal growth architecture of Central European plants,” 2006, http://clopla.butbn.cas.cz/

M. Kleyer, R. M. Bekker, I. C. Knevel et al., “The LEDA traitbase: a database of life-history traits of the Northwest European flora,” Journal of Ecology, vol. 96, no. 6, pp. 1266–1274, 2008.

H. Ellenberg, H. E. Weber, R. Düll, V. Wirth, W. Werner, and D. Paulissen, Zeigerwerte von Pflanzen in Mitteleuropa, vol. 18, Scripta Geobotanica, 2nd edition, 1992.

T. G. Tutin, V. H. Heywood, N. A. Burges, D. H. Valentine, S. M. Walters, and D. A. Webb, Flora Europaea, vol. 5 of Set and CD-ROM Pack, Cambridge University Press, Cambridge, UK, 2nd edition, 2001.

I. Kühn, W. Durka, and S. Klotz, “BiolFlor—A new plant—trait database as a tool for plant invasion ecology,” Diversity and Distributions, vol. 10, no. 5-6, pp. 363–365, 2004.

M. Kuhn, “Building predictive models in R using the caret package,” Journal of Statistical Software, vol. 28, no. 5, pp. 1–26, 2008.

M. Kuhn, “caret: Classification and Regression Training. R package version 4.98,” 2011, http://CRAN.R-project.org/package=caret

T. Jamil, W. Ozinga, and C. J. F. ter Braak, “Selecting traits that explain species-environment relationships: a generalized linear mixed model approach,” submitted

N. L. Poff, J. D. Olden, N. K. M. Vieira, D. S. Finn, M. P. Simmons, and B. C. Kondratieff, “Functional trait niches of North American lotic insects: traits-based ecological applications in light of phylogenetic relationships,” Journal of the North American Benthological Society, vol. 25, no. 4, pp. 730–755, 2006.

G. C. Fernandez, “Model selection in PROC MIXED—a user—friendly SAS® macro application,” in Proceedings of the SAS Global forum, Orlando, Fla, USA, April 2007.

R. C. Littell, G. A. Milliken, W. W. Stroup, R. D. Wolfinger, and O. Schabenberger, SAS for Mixed Models, SAS Institute, Cary, NC, USA, 2006.

M. Yuan and Y. Lin, “Efficient empirical bayes variable selection and estimation in linear models,” Journal of the American Statistical Association, vol. 100, no. 472, pp. 1215–1225, 2005.

B. T. West, K. B. Welch, and A. T. Galecki, Linear Mixed Models: A Practical Guide Using Statistical Software, Chapman & Hall/CRC, 2006.

P. J. Diggle, P. Heagerty, K. Y. Liang, and S. L. Zeger, Analysis of Longitudinal Data, Oxford University Press, Oxford, UK, 2002.

H. Akaike, “Information theory and an extension of the maximum likelihood principle,” in 2nd International Symposium on Information Theory, Tsahkadsor, Armenian SSR, 1973.

C. M. Hurvich and C. Tsai, Regression and Time Series Model Selection in Small Samples, Oxford University Press, 1989.

H. Bozdogan, “Model selection and akaike's information criterion (AIC): the general theory and its analytical extensions,” Psychometrika, vol. 52, no. 3, pp. 345–370, 1987.

S. Gideon, “Estimating the Dimension of a Model,” The Annals of Statistics, vol. 6, no. 2, pp. 461–464, 1978.

M. J. Gurka, “Selecting the best linear mixed model under REML,” The American Statistician, vol. 60, no. 1, pp. 19–26, 2006.

K. P. Burnham and D. R. Anderson, Model Selection and Multimodel Inference : A Practical Information-Theoretic Approach, Springer, New York, NY, USA, 2002.

K. W. Broman and T. P. Speed, “A model selection approach for the identification of quantitative trait loci in experimental crosses,” Journal of the Royal Statistical Society, vol. 64, no. 4, pp. 641–656, 2002.

D. Bates, M. Maechler, B. Bolker, et al., “lme4: Linear mixed-effects models using S4 classes. R package version 0.999375-39,” 2011, http://CRAN.R-project.org/package=lme4

G. Verbeke and G. Molenberghs, Linear Mixed Models for Longitudinal Data, Springer, New York, NY, USA, 2000.

S. Xu, “An empirical bayes method for estimating epistatic effects of quantitative trait loci,” Biometrics, vol. 63, no. 2, pp. 513–637, 2007.

C. Bishop, Pattern recognition and machine learning (information science and statistics), Springer, New York, NY, USA, 2006.

M. E. Tipping, “Sparse bayesian learning and the relevance vector machine,” Journal of Machine Learning Research, vol. 1, pp. 211–244, 2001.

J. Berger, Statistical Decision Theory and Bayesian Analysis, Springer, 1985.

P. Legendre, R. G. Galzin, and M. L. Harmelin-Vivien, “Relating behavior to habitat: solutions to the fourth-corner problem,” Ecology, vol. 78, no. 2, pp. 547–562, 1997.

D. Chessel, A. B. Dufour, and S. Dray, “ade4: Analysis of Ecological Data : Exploratory and Euclidean methods in Environmental sciences. R package version 1.4-17,” 2011, http://CRAN.R-project.org/package=ade4

A. Cormont, C. C. Vos, C. A. M. van Turnhout, R. P. B. Foppen, and C. J. F. Ter Braak, “Using life-history traits to explain bird population responses to changing weather variability,” Climate Research, vol. 49, no. 1, pp. 59–71, 2011.

B. Shipley, P. A. Keddy, D. R. J. Moore, and K. Lemky, “Regeneration and establishment strategies of emergent macrophytes,” Journal of Ecology, vol. 77, no. 4, pp. 1093–1110, 1989.

C. J. F. ter Braak, A. Cormont, and S. Dray, “Improved testing of species traits-environment relationships in the fourth corner problem,” Ecology. In press

J. P. M. Lenssen, F. B. J. Menting, and W. H. V. D. Putten, “Plant responses to simultaneous stress of waterlogging and shade: amplified or hierarchical effects?” New Phytologist, vol. 157, no. 2, pp. 281–290, 2003.

S. Dolédec, B. Statzner, and M. Bournard, “Species traits for future biomonitoring across ecoregions: Patterns along a human-impacted river,” Freshwater Biology, vol. 42, no. 4, pp. 737–758, 1999.

M. Westqby, E. Jurado, and M. Leishman, “Comparative evolutionary ecology of seed size,” Trends in Ecology & Evolution, vol. 7, no. 11, pp. 368–372, 1992.

J. P. Grime, Plant Strategies and Vegetation Processes, Wiley, Chichester, UK, 1979.

N. C. Garwood, “Functional morphology of tropical tree seedlings,” in The Ecology of Tropical Forest Tree Seedlings, M. D. Swaine, Ed., pp. 59–129, UNESCO, Paris, France, 1996.

P. T. Green and P. A. Juniper, “Seed mass, seedling herbivory and the reserve effect in tropical rainforest seedlings,” Functional Ecology, vol. 18, no. 4, pp. 539–547, 2004.

P. A. Keddy, L. Twolan-Strutt, and I. C. Wisheu, “Competitive effect and response rankings in 20 wetland plants: are they consistent across three environments?” Journal of Ecology, vol. 82, no. 3, pp. 635–643, 1994.

A. Stockey and R. Hunt, “Predicting secondary succession in wetland mesocosms on the basis of autecological information on seeds and seedlings,” Journal of Applied Ecology, vol. 31, no. 3, pp. 543–559, 1994.

P. J. Grubb, “Control of forest growth and distribution on wet tropical mountains: with special reference to mineral nutrition,” Annual Review of Ecology and Systematics, vol. 8, pp. 83–107, 1977.

M. Henry, H. Stevens, D. E. Bunker, S. A. Schnitzer, and W. P. Carson, “Establishment limitation reduces species recruitment and species richness as soil resources rise,” Journal of Ecology, vol. 92, no. 2, pp. 339–347, 2004.