Issue 
Ann. For. Sci.
Volume 67, Number 3, May 2010



Article Number  307  
Number of page(s)  13  
Section  Original articles  
DOI  https://doi.org/10.1051/forest/2009113  
Published online  18 February 2010 
Original article
Estimating growth in beech forests: a study based on long term experiments in Switzerland
Estimation de la croissance dans les hêtraies : une étude basée sur des expérimentations à long terme en Suisse
^{1}
Universidad de Santiago de Compostela, Unidad de Gestión Forestal
Sostenible, Departamento de Ingeniería Agroforestal, Lugo, Spain
^{2}
Swiss Federal Institute for Forest, Snow and Landscape Research
WSL, Birmensdorf,
Switzerland
^{3}
GeorgAugustUniversity Göttingen, Faculty of Forest Sciences and Forest
Ecology, BurckhardtInstitute, Göttingen, Germany
^{*} Corresponding author: Escuela Politécnica Superior, Campus
Universitario s/n. 27002 Lugo, Spain; juangabriel.alvarez@usc.es
Received:
13
May
2009
Accepted:
15
September
2009
• This contribution presents a dynamic stand growth model for Beech (Fagus sylvatica L.) forests, based on a dataset provided by the Swiss Federal Institute for Forest, Snow and Landscape Research WSL, Birmensdorf. The dataset includes 143 research plots, covering a wide range of growing sites and providing up to 16 interval measurements per research plot.
• The objective of this research is to complement the range of existing beech growth models by bridging the gap between the historical yield tables and the single tree growth models. The specific aim is to develop transition functions which will project three state variables (dominant height, basal area and number of trees per hectare) at any particular time, in response to any arbitrary silvicultural treatment.
• Two of the transition functions were derived using the generalized algebraic difference approach (GADA), the third one was derived with the algebraic difference approach (ADA). All the functions were fitted simultaneously using iterative seemingly unrelated regression and a baseageinvariant method. The influence of thinnings on basal area growth was included by fitting different transition functions for thinned and unthinned stands.
• The overall model provides satisfactory predictions for time intervals up to 20 years. The new model is robust and its relatively simple structure makes it suitable for economic analysis and decision support.
Résumé
• Cette contribution présente un modèle dynamique de croissance des peuplements de hêtres (Fagus sylvatica L.), basé sur un ensemble de données fournies par l’Institut Fédéral Suisse de Recherche sur la Forêt, la Neige et le Paysage, WSL à Birmensdorf. L’ensemble des données comprend 143 parcelles de recherche, couvrant un large éventail de sites et fournissant jusqu’à 16 intervalles de mesures par parcelle de recherche.
• L’objectif de cette recherche est de compléter la gamme de modèles de croissance du hêtre existants, en jetant un pont entre les tables de production historiques et les modèles de croissance d’arbre. L’objectif spécifique est de développer des fonctions de transition qui projeterons trois variables d’état (hauteur dominante, surface terrière et nombre d’arbres par hectare) à n’importe quel moment déterminé, en réponse à n’importe quel traitement sylvicole arbitraire.
• Deux des fonctions de transition ont été calculées en utilisant l’approche différence algébrique généralisée (GADA), la troisième a été dérivée de l’approche différence algébrique (ADA). Toutes les fonctions ont été ajustées en utilisant simultanément une régression itérative sans lien apparent et une méthode basée sur l’invariance de l’âge. L’influence des éclaircies sur la croissance de la surface terrière a été inclue en ajustant différentes fonctions de transition pour les peuplements éclaircis et les peuplements non éclaircis.
• Le modèle général fournit des prédictions satisfaisantes pour des intervalles de temps jusqu’à 20 ans. Le nouveau modèle est robuste et sa structure relativement simple fait qu’il est convient pour l’analyse économique et l’aide à la décision.
Key words: thinning effect / simultaneous fitting / GADA / Fagus sylvatica
Mots clés : effet de l’éclaircie / ajustement simultané / GADA / Fagus sylvatica
© INRA, EDP Sciences, 2010
1. INTRODUCTION
The European Beech (Fagus sylvatica L.) is one of the most important tree species in Europe. Its natural distribution range is defined by a maritime climate with moderate fluctuations of temperature and precipitation. The species avoids areas with extended very cold and dry periods and prefers moist sites. The natural range where beech is dominant and the phytogeographic range where it occurs, have been described in numerous research papers and more than 20 maps of the range have been published (Bolte et al., 2007).
In Swiss forests, European beech is the second species by growing stock with more than 18% (about 75 million m^{3}) and the annual gross increments of this species between 1985 and 1995 represents the 20% of the total Swiss forests (SAEFL/WSL, 2005). This species will very possibly play an important role in Swiss forestry, because is the main natural species in lower altitudes where natural forests are politically demanded and it can produce highly valuable timber on good sites.
The yield tables developed by Paulsen (1795) possibly represent the first serious approach to beech growth modeling. These and subsequent generations of yield tables portray the development of pure beech stands in fixed (usually 5yearly) time intervals for a given silvicultural treatment and growing site.
Due to improved computer technology and increasing information needs, the yield tables are being gradually replaced by flexible growth models. Several single tree growth simulators were developed during the past two decades to evaluate alternative silvicultural strategies (e.g. Nagel et al., 2002; Sterba and Monserud, 1997). Nevertheless, to date, the only growth model for the species in Switzerland is a yield table developed by Badoux (1983). This model provides limited information about the forest stand and does not reflect accurately the evolution under different stand density conditions. During the past years, there has been considerable interest in testing the applicability of individual tree models for a wide variety of Swiss forests in order to provide a potential alternative to the commonly used yield tables. However, for beech forests, the diameter increment was underestimated (Hallenbarter et al., 2005) and the lack of detailed stand structure information, not always available, had a strong effect on the quality of the simulation results (Schmid et al., 2006).
Our research has shown that it is possible to develop growth models which predict population mean parameters (e.g. dominant height) and areabased quantities (e.g. basal area per ha) with a comparatively high degree of accuracy. Thus, the objective of this contribution is to complement the range of existing beech growth models with a populationbased approach, thus bridging the gap between yield tables and single tree models. We are using the extensive data base provided by the Swiss Federal Institute for Forest, Snow and Landscape Research to illustrate our statespace modeling approach.
2. MATERIAL AND METHODS
2.1. Data
The data used in this study originate form the growth and yield database of the Swiss Federal Institute for Forest, Snow and Landscape Research WSL in Birmensdorf Switzerland. The data were collected on 143 longterm pure and evenaged beech research plots (Fig. 1) on 103 locations, most of them on altitudes below 800 m in the Jura region, the Swiss plains and the PreAlps.
Figure 1
Map of all beech growth and yield plots (black dots = beech plots with tree height measurements, white dots = plots without height measurements). Map reproduced with permission of swisstopo (JA082265). 
Some of the plots which are concentrated at the same location are thinning experiments where different thinning types were applied. In these experiments, the treatments included A, B, C, D and Hgrade thinnings. These data reflect the response of beech stands of different ages growing on different sites to a variety of thinning regimes. The dataset used in this study covers a wide range of treatments, including extreme weights and types of removals. It is the best available source of information on beech growth in Switzerland.
The plots were remeasured between one and 17 times from 1889 to 2001. The shortest observation interval between plot remeasurements is two years, the longest 111 years. During the field enumeration, breast height diameters of all trees in the plot were measured crosswise with a caliper at a height of 1.3 m and the exact place of measurement was marked. Total heights and heights to crown base were assessed for selected trees including dominant trees, but not in all plots enumerations. Dominant height, H (m) was defined as arithmetic mean height of the 100 thickest trees per hectare and it was calculated for each plot where total heights were measured. In addition, descriptive variables such as the relative social position within the plot were assessed. Five of the plots are still being remeasured. The mean, maximum and minimum values and the standard deviation of the main stand variables are shown in Table I.
Summary statistics of the dataset used for model development.
2.2. Modeling approach
The method used in this analysis is based on the statespace approach described by García (1994), which makes use of state variables to characterize the system at an initial stage, and transition functions to project all or some of the state variables to future states.
2.2.1. State variables
In selecting the state variables, the principle of parsimony must be taken into account (García, 1994; Vanclay, 1995): the model should be the simplest one that describes the biological phenomena and that remains consistent with the structure and function of the actual biological system. A two dimensional vector including dominant height (H) and stand basal area (G) as explanatory variables may be sufficient to describe the state of the stand at a given age (t), and the advantage of including additional variables may not be justified (García, 1988). Nevertheless, in situations covering a wide range of silvicultural regimes, it may be necessary to include an additional variable such as the number of trees per unit area (N) (CastedoDorado et al., 2007a; DieguezAranda et al., 2006; García, 1994). Taking into account that the effect of thinning on stand basal area growth was analyzed in this study, the stand conditions at any point in time are defined by three state variables (dominant height, number of trees per hectare and stand basal area) and the model include three transition functions. The outputs of these three transition functions provide the inputs for a static function, which predicts the total volume of the stand for a given age.
The three transition functions must possess some important properties to provide consistent estimates (García, 1994): (i) consistency, meaning no change for zero elapsed time; (ii) pathinvariance, where the result of projecting the state first from t_{0} to t_{1}, and then from t_{1} to t_{2}, must be the same as that of the onestep projection from t_{0} to t_{2} and (iii) causality, in that a change in the state can only be influenced by inputs within the relevant time interval. Dynamic equations generated by integration of differential equations (or summation of difference equations when using discrete time) automatically satisfy these conditions.
Fulfilment of the above mentioned requirements depends on both the construction method and the specific mathematical function. A useful theory is the algebraic difference approach (ADA) developed by Bailey and Clutter (1974) and its generalization known as GADA (Cieszewski and Bailey, 2000). Dynamic equations have the general form (omitting the vector of model parameters) of Y = f(t, t_{0}, Y_{0}), where Y is the value of the function at age t, and Y_{0} is the reference variable defined as the value of the function at age t_{0}. The ADA essentially involves replacing a basemodel sitespecific parameter with its initialcondition solution. The GADA allows expanding the base equations according to various theories about growth characteristics (e.g., asymptote, growth rate), thereby allowing more than one parameter to be sitespecific and allowing derivation of more flexible dynamic equations (see Cieszewski, 2002; 2003).
In the present study, the transition functions for dominant height and stand basal area growth were derived with the generalized algebraic difference approach (GADA) and the transition function for reduction in tree number was derived with the algebraic difference approach (ADA).
An effective system for parameter estimation needs to be based on identifying individual trends represented in the data (Bailey and Clutter, 1974; Cieszewski, 2003). These trends represented in dominant height, number of trees per hectare and stand basal area can be modeled by considering that individuals’ responses all follow a similar functional form with parameters that vary among individuals (local parameters) and parameters that are common for all individuals (global parameters). In practice both baseage specific (BAS) and baseage invariant (BAI) estimation methods can be used. However, only BAI methods assume that the data measurements always contain measurement and environmental errors (both on the left and right hand sides of the model) that must be modeled (Cieszewski, 2003). Therefore, a BAI parameter estimation technique was used and from among the different BAI methods available the dummy variables method described in Cieszewski et al. (2000) was selected.
Based on the idea that one model parameter is specific for each individual (sample plot), and that the remaining parameters are shared by all individuals, Cieszewski (2003) proposed a simplified approach of mixedeffects model for transition functions development. The method of Cieszewski (2003), which we employed in this study, results from an adaptation of the nonlinear mixedeffects model definition by Lindstrom and Bates (1990), whose general formulation can be written as: Y_{ij} = f(A_{i}β + B_{i}b_{i}, t_{ij}) + e_{ij}, where Y_{ij} is the jth response on the ith individual (i.e., dominant height or basal area of plot i at age measurement j); t_{ij} is the predictor for the jth response on the ith sample plot; f is a nonlinear function of the predictor and parameter vectors β (a pvector of fixed population parameters –global parameters common to all sample plots) and b_{i} (a qvector of the random effects associated with sample plot i – sitespecific parameters); A_{i} is a design matrix of size r × p for the global parameters; B_{i} is a design matrix of size r × q for the sitespecific parameters; r is the length of the vector of global parameters; and e_{ij} is a normally distributed random noise with zero mean.
Finding a model form that can describe the data well using only one parameter that varies with site is the first step of the method by Cieszewski (2003). If more than one parameter has to be sitespecific, the varying parameters must be correlated with each other and defined as functions of fixed parameters and just one common varying parameter. After such a reparameterization, all except the varying parameter are uncorrelated parameters of the base model, and can be defined as pure fixed effects. The varying parameter is the unique parameter related to the random site effects (Lindstrom and Bates, 1990). The sitespecific parameter of each sample plot can be defined as mixedeffects (both fixed and random effects) comprising the sum of the global mean for this parameter (ψ) plus a sitespecific effect (ξ_{i}).
The next steps of the method by Cieszewski (2003) involve: (i) applying the GADA to ψ + ξ_{i} defined as X (a theoretical vector variable defining site quality with mean ψ and sitespecific adjustments ξ_{i} with unknown distributional properties and an unknown magnitude of values); (ii) solving for X; (iii) substituting the solution for X in the explicit threedimensional site equation (Y = f(t,X)) for initial conditions t_{0} and Y_{0} (Y = f(t,t_{0},Y_{0}), where Y_{0} comprises for each sample plot the sum of the observed variable (dominant height or basal area) at age t_{0} plus a sitespecific effect e_{0i} interpreted as the random measurement and environmental errors); and (iv) fitting in one stage using a stochastic regression approach, or a varying parameter regression with desired fitting criteria. In the last fitting procedure, an estimation of the variables (i.e., Y_{0} at t_{0}) resulting from fitting a curve to the observed trend in the data for each individual growth series are used as predictors instead of the observed variables. This direct use of data, as constants, does not violate regression assumptions because the measurement and environmental errors associated with these data are estimated simultaneously with all the other parameters of the model (Cieszewski, 2003).
As Cieszewski et al. (2000) pointed, the described baseageinvariant approaches are similar to mixedeffect modeling. The socalled local, or sitespecific, parameters are in fact the random effects with a nonzero mean equal to the fixed effect representing average reference dominant height, basal area or number of stems per hectare at the given base age. This parameter is not modeled explicitly, but it can be easily recovered after the model fitting. All the other fixed effects are simply the global model parameters.
2.2.2. Transition functions for height and basal area
Dominant height and stand basal area growth models are key components of stand models, since they are directly related to economic variables such as stand volume. A number of researchers have pointed out desirable attributes of these two transition functions (Bailey and Clutter, 1974; Cieszewski and Bailey, 2000; Clutter et al., 1983). The most frequently listed criteria included the three general properties of transition functions described above and some specific requirements: polymorphism (for dominant height growth curves), sigmoid growth pattern with an inflexion point, horizontal asymptote at old ages (multiple asymptote for dominant height growth curves), logical behaviour (height should be zero at age zero and equal to site index at the reference age; the curve should never decrease), theoretical basis or interpretation of model parameters derived by analytically tractable algebraic operations and baseage invariance. The fulfilment of these attributes depends on both the model construction method and the statistical procedures used for model fitting, and cannot always be achieved.
Three different base equations were evaluated for modelling dominant height and stand basal area growth pattern: Bertalanffy (Bertalanffy, 1957); Loglogistic and Korf (1939). These functions have been widely used in modelling stand height and stand basal area growth (e.g., Álvarez González et al., 2005; Barrio Anta et al., 2006; Cieszewski, 2002; DiéguezAranda et al., 2006; Monserud, 1984).
As general notational convention, a_{1}, a_{2},…, a_{n} were used to denote parameters in base models, while b_{1}, b_{2},…, b_{m} were used for global parameters in subsequent GADA formulations. All the GADA based models have the general implicit form of Y = f (t, t_{0}, Y_{0}, b_{1}, b_{2}, ..., b_{m}).
The first GADA formulation was tested on the basis of the base model proposed by Bertalanffy (1957), whose integral form can be, for simplicity, represented as (1)where a_{1} is an asymptote or limiting value, a_{2} is often called a “rate parameter”, and a_{3} is often referred to as “an initial pattern parameter”. To derive models with both polymorphism and variable asymptotes from Equation (1), more than one parameter has to be sitespecific. Thus, both the asymptote a_{1} and the shape parameter a_{3} are assumed to be dependent on X. To facilitate such derivation the base equation is reparameterized into a form more suitable for manipulation of these two parameters (using exp (a′_{1}) instead of a_{1} and taking the natural logarithm of the function) as follows: The site parameters can be conditioned to be consistently proportional to each other’s inverse over the site productivity dimension by the following definitions: which defines the following relationship: (2)For this formulation, the solution for X involves finding roots of a quadratic equation and selection of the most appropriate root to substitute into the dynamic equation. The selection of the most appropriate expression for X may depend on the equation parameters that in turn depend on the data and the domain of the applicable ages (Cieszewski and Bailey, 2000). The solution for X in Equation (2) with initial condition values t_{0} and Y_{0} is where Selecting the root more likely to be real and positive (i.e., the one involving addition rather than subtraction of the square root), and substituting it into Equation (2) results in the first GADA formulation used in this study: (E1)The next two GADA formulations were derived on the base of the loglogistic model: (3)where a_{1} is the asymptote parameter, a_{2} is the halfsaturation parameter that defines the value of exp( − u) at which Y(u) = a_{1} / 2, and u is usually a linear combination of logtransformations of the independent variable(s) (Monserud, 1984; Cieszewski, 2002).
Within the different forms investigated, the GADA based formulation by Cieszewski (2002) was found to perform particularly well. It considers the abovementioned u as a logarithmically transformed function of age, such that u = a_{3}lnt. Therefore, Equation (3) can be simplified as: (4)In a first approach, Cieszewski (2002) replaced a_{1} with a constant plus the unobserved site variable X, and a_{2} by b_{2} / X, becoming the next formulation (5)The solution for X involves again finding roots of a quadratic equation. With initial condition values t_{0} and Y_{0} it is Selecting the root more likely to be real and positive (i.e., the one involving addition rather than subtraction of the square root), and substituting it into Equation (5) results in the second GADA formulation: (E2)A second GADA formulation can be obtained based on Equation (4) by replacing a_{1} with a constant plus the unobserved site variable X, and a_{2} by b_{2}X(6)The solution for X in Equation (6) with initial condition values t_{0} and Y_{0} is Substituting it into Equation (6) results in the third GADA formulation: (E3)The fourth transition function is based on the base model proposed by Korf (1939): where a_{1} is the asymptote parameter, and the other two parameters (a_{2} and a_{3}) are related to the inflexion point and the growth rate.
Cieszewski (2002) replaced a_{1} with the exponential of the unobserved site variable X, and a_{2} by b_{1} / X(7)The solution for X involves taking the natural logarithm of the functions and again finding roots of a quadratic equation. With initial condition values t_{0} and Y_{0} it is Selecting the root more likely to be real and positive (i.e., the one involving addition rather than subtraction of the square root), and substituting it into equation (7) results in the fourth GADA formulation: (E4)Different methods have been used to take into account the thinning effects in basal area growth. Some authors (Barrio et al., 2006; Knoebel et al., 1986; Zarnoch et al., 1991) have developed different equations for thinned and unthinned stands that have the same mathematical structure but that have been parameterized using different data sets. In the other hand, other studies (e.g. Hasenauer et al., 1997; Hynynen, 1995) have analyzed the inclusion of a thinning response function to express the basal area growth of a thinned stand as a product of a reference growth and the thinning response function: the reference growth accounts for the factors affecting stand growth in unthinned stands while the thinning response function predicts the relative growth response following thinning. In this work, we used the first approach to take into account the effect of thinning on stand basal area growth.
Mathematical expression of some mortality models derived from differential equations (8) to (10). Where N and N_{0} are the number of trees per hectare present at ages t and t_{0}, respectively.
2.2.3. Transition function for reduction in number of trees
One of the most important transition functions of a dynamic growth and yield model is a mortality model that estimates the natural decline in number of trees caused by stand density and other environmental factors. However, mortality remains one of the least understood components of natural processes growth. Many functions have been used to model empirical mortality equations. Some of them are mathematical relationships between stem number reduction and stand variables (e.g. Álvarez González et al., 2004; Rennolls and Peace, 1986), others are biologically based functions derived from differential equations. These biologically based functions have properties that are essential in a mortality model but are not always present in a pure mathematical model (Clutter et al., 1983): consistency, pathinvariance and asymptotic limit of stocking approaching zero as age becomes very large. Also, for evenaged stands it is reasonable to assume that ingrowth is negligible, so if age at time two is greater than age at time one, the density at time two will be less than density at time one. Specification of a difference equation model that possesses these properties and properly expresses the population relationship is often a difficult task (Clutter et al., 1983). In this work, the equation for estimating reduction in tree number was developed on the basis of the following differential function: (8)(9)(10)where N is the number of trees per hectare present at age t, ΔN / Δt is the instantaneous mortality rate operating at age t, and α, β, γ and δ the model parameters that regulate the mortality rate.
The three equations (Eqs. (8)–(10)) consider the mortality rate related to number of trees per hectare and age. Five different transition functions derived from these three differential functions have been used in this study, which are shown in Table II.
2.2.4. Model fitting
The statespace approach has been widely used to develop stand growth models (e.g. CastedoDorado et al., 2007a; DiéguezAranda et al., 2006; NordLarsen and Johannsen, 2007). Usually, each transition function is fitted separately, however the transition functions for dominant height, basal area and number of trees reduction together define a system of three regression equations and the random errors of the equations could be correlated. In order to improve the efficiency of the estimation the iterative seemingly unrelated regression (ISUR) method, that takes into account the crossequation correlations, was used. The estimate of the crossequation error covariance matrix to initiate the iterative procedure was obtained using Ordinary Least Squares. A similar approach was used by NordLarsen and Johannsen (2007) to develop a stand growth model for European beech in Denmark.
In the general formulation of the dynamic equations, the error terms e_{ij} are assumed to be independent and identically distributed with zero mean. Nevertheless, because of the longitudinal nature of the data sets used for model fitting, correlation between the residuals within the same plot may be expected. To overcome any possible autocorrelation, we modelled the error terms using a 1order continuous time autoregressive error structure (CAR(1)), which allows the model to be applied to irregularly spaced, unbalanced data (Gregoire et al., 1995): where e_{ij} is the jth ordinary residual on the ith sample plot, d_{1} = 1 for j > 1 and d_{1} = 0 for j = 1, e_{ij − 1} is the j − 1th ordinary residual on the ith sample plot, ρ_{1} is the autoregressive parameter to be estimated, and t_{ij} − t_{ij − 1} is the time that separates the jth from the j − 1th observations within each sample plot. In such cases ε_{ij} now includes the error terms under conditions of independence.
To evaluate the presence of autocorrelation, graphs representing residuals plotted against lagresiduals from previous observations within each plot were examined visually. The heteroscedasticity was also examined visually, by plotting residuals as a function of predicted values, in order to investigate possible weighting factors.
The fitting of the equations were accomplished by SAS/ETS^{®} MODEL procedure (SAS Institute Inc., 2004). The correction for autocorrelation was programmed in the MODEL procedure by adding to the model formulation the term di*(rho**dist)*zlag1(resid.V), where di is a dummy variable which values is 0 for the first observation of each sample plot and 1 in other case; rho is the autoregressive parameter (ρ_{1}); dist is the time that separates the the jth from the j − 1th observations within each sample plot and zlag1(resid.V) is a SAS function referring to the previous residual of the stand variable V (dominant height, basal area or number of trees per hectare).
2.2.5. Model comparison and evaluation
The comparison of the estimates for the different equations fitted for each stand variable was based on numerical and graphical analyses of the residuals (e_{ij}). Three statistical criteria obtained from them were examined: bias (), which tests the systematic deviation of the model from the observations; root mean square error (RMSE), which analyses the accuracy of the estimates; and the model efficiency (MEF), which shows the proportion of the total variance that is explained by the model, adjusted for the number of model parameters and the number of observations. The expressions of these statistics are as follows: where y_{i}, and are the observed, predicted and average values of the dependent variable, respectively; n is the total number of observations used to fit the function; and p is the number of model parameters.
Other important steps in evaluating the functions fitted were graphical analysis of the residuals and examination of the appearance of the fitted curves overlaid on the trajectories of the dependent variables for each plot. Visual or graphical inspection is an essential point in selecting the most appropriate model because curve profiles may differ drastically, even though fit statistics and residuals are similar (Huang et al., 2003).
Once the complete model has been constructed, assessment of their validity is often needed to ensure that the predictions represent the most likely outcome in the real world (Yang et al., 2004). Moreover, although the behaviour of the functions for each stand variable within the model plays an important role in determining the overall outcome, the validity of each individual component does not guarantee the validity of the overall outcome, which is usually considered more important in practice. The overall model outcome must therefore also be evaluated. The only method that can be regarded as “true” validation involves the use of a new independent data set (Pretzsch et al., 2002; Yang et al., 2004) but the scarcity of such data forces the use of alternative approaches. The common method of splitting the data set in two portions does not provide additional information (Huang et al., 2003) and neither is recommended from the point of view of parameter estimation.
Moreover, other techniques such as double crossvalidation or statistical testing provide very limited information about the predictive ability of the models (Kozak and Kozak, 2003; Yang et al., 2004). Therefore, we decided to defer model validation until a new data set is available for assessing the quality of the predictions.
The functions of each stand variable were evaluated to confirm their usefulness and reliability by the goodnessoffit statistics and the graphical analysis described above (refer also to Fig. 2). In addition, plots of observed against predicted values of each state variable were inspected.
Figure 2
HResiduals versus HLag1Residuals for model E1 (first row), GResiduals versus GLag1Residuals for model E1 (second row) and NResiduals versus NLag1Residuals for model E5 (third row) fitted using a firstorder continuous time autoregressive error structure. 
Parameter estimates and goodnessoffit statistics for the transition functions for dominant height considered.
To evaluate the maximum time interval for accurate projections, we used the critical error (E_{crit}). In this study, the critical error is expressed as a percentage of the observed mean basal area and is computed by rearranging Freese’s statistic (Reynolds, 1984): where n is the total number of observations in the data set, y_{i} is the observed basal area, the basal area predicted from the fitted model, is the average of the observed basla area, τ is a standard normal deviate at the specified probability level (τ = 1.96 for α = 0.05), and is obtained for α = 0.05 and n degrees of freedom.
3. RESULTS AND DISCUSSION
At first, the transition functions for H, G and N were fitted separately to select the best equation. The models were fitted using nonlinear least squares without expanding the error terms to account for autocorrelation. However, a slight trend in residuals as a function of lag1 residuals within the same sample plot was apparent in all the models. After applying a firstorder continuous time autoregressive error structure, all the parameters associated with the error term were significant at the 5% level, and the error trends in residuals disappeared (Fig. 2). The sole purpose of the autocorrelation correction was to improve the interpretation of the statistical properties of the models, and it has no use in practical applications unless one is working repeatedly on the same sample plot.
The parameter estimates for each model, including its standard errors and p values, and the corresponding goodnessoffit statistics are shown in Tables III to V for dominant height, basal area and reduction in number of trees, respectively. After considering the statistical criteria and analysing the graphical output, model E4, was selected as the best transition function for dominant height and basal area and model E5 for reduction in number of trees.
Although model E4, derived from the Korf function, does not achieve the best values based on the statistical criteria which were used to compare the transition functions for dominant height (Tab. III), in the graphical analysis this model showed the best performance.
Parameter estimates and goodnessoffit statistics for the transition functions for basal area considered.
For basal area growth, all the models explained more than 98% of the total variance, and model E4, derived from the Bertalanffy equation obtained the best values of the statistical criteria used (Tab. IV). However, as in the case of dominant height growth, the model E4, derived from the Korf equation showed a better performance.
Parameter estimates and goodnessoffit statistics for the transition functions for reduction in tree number considered.
Differences in basal area growth between thinned and unthinned sample plots were detected, therefore, different transition functions are proposed for each treatment. These results suggest that, for our data set, the stand basal area growth rate following a thinning exceeded the stand basal area growth rate of a stand with similar initial stand conditions which had not been recently thinned. This “thinning effect”agrees with several other studies (e.g., Amateis, 2000; Hasenauer et al., 1997; Pienaar and Shiver, 1984; Pienaar et al., 1985; Sharma et al., 2006). These studies involving different species have demonstrated that there is a difference, however slight, in the basal area growth in thinned and unthinned stands of the same initial age, site index and stand basal area.
This result, however, contradicts the findings of other studies. Barrio et al. (2006) and CastedoDorado et al. (2007b) proposed the same equation for thinned and unthinned stands assuming that the thinning effect was an inherent part of the model. The authors concluded that any contradictory results relating to the “thinning effects theory”may be at least partly attributed to the specific structure of the experimental data sets. Some data sets include only a narrow range of thinning trials where mainly low or moderate thinnings were carried out.
Due to the fact that the data used in this study do not provide enough information about the time interval during which the thinning effect is relevant for basal area growth, we propose that the transition function for basal area in thinned stands should not be used for projection intervals longer than 5 years.
Table V shows the parameter estimates for each of the 5 equations for the reduction in stem number and the statistics used to compare them. Models E5 to E7 are derived from a differential function considering that the relative rate of change in the stem number is proportional to a power function of age. Models E8 and E9 are derived from differential equations considering a hyperbolic and an exponential function of age, respectively.
Taking into account the statistics used to compare the models, there is no clear evidence of how to include age as an explanatory variable. However, only for model E5 all the parameters were significant. Therefore, this model was selected for fitting all the equations simultaneously. This equation implies that the relative rate of change in the number of trees is proportional to a power function of age. In the past, similar results were obtained for different species (e.g. Clutter and Jones, 1980; DieguezAranda et al., 2005; Woollons, 1998).
Figure 3
Upper left: dominant height curves for site indices of 18, 22, 26, 30 and 34 m at a reference age of 80 years for the GADA equation E1. Upper right: reduction in tree number curves for 400, 300, 2 200, 3 100, 4 000 trees per ha at a reference age of 40 years for the model (E5). Down: basal area growth rates for thinned (left) and unthinned stands (right) as affected by initial stem number and age overlaid on the trajectories of basal area. 
In a second approach, the selected models for H, G and N were fitted simultaneously to minimize the total sum of square errors in the system and to take the crossequation error correlations into account. When ordinary least squares was used to fit the models simultaneously, the correlation value of the errors was only significant for the dominant heightbasal area models (ρ = 0.3252). This value was reduced by 21.7% using iterative seemingly unrelated regression. The graphical analysis of residual H, G and N against the corresponding predicted values showed no trends.
The transition function for dominant height proposed for beech stands finally obtained by the simultaneous fitting is the following: where Y is the predicted dominant height (m) at age t (years), and Y_{0} is the predictor dominant height at age t_{0}.
This model explained 98.1.5% of the total variance of the data, and its RMSE was 1.01 m. In selecting the optimum base age, it was found that 80 years was superior for predicting height for other ages. The curves for dominant height development for 18, 22, 26, 30, 34 m at a reference age of 80 years overlaid on the profile plots of the data set are shown in Figure 3.
The proposed transition functions for basal area for unthinned and thinned beech stands, respectively, obtained by the simultaneous fitting are the following: where Y is the predicted basal area (m^{2}/ha) at age t (years), and Y_{0} is the predictor basal area at age t_{0}.
The basal area transition function accounted for more than 98% of the total variance of the data with a RMSE of 1.27 m^{2}/ha. The curves of basal area growth rates for thinned and unthinned stands as affected by initial basal area and age, superimposed over the trajectories of observed values over time are shown in Figure 3.
Finally, the proposed equation for estimating the natural stem number decline for beech stands is: where N is the predicted number of trees per hectare at age t (years), and N_{0} is the predictor number of trees per hectare at age t_{0}.
The above model explained approximately 96.4% of the total variance of the data and the RMSE was 83.4 trees/ha. The trajectories of observed and predicted number of trees over time for 400, 1 300, 2 200, 3 100 and 4 000 trees per hectare at a reference age of 40 years are shown in Figure 3.
To achieve a satisfactory overall evaluation of the model, the observed basal area from any one inventory of the sample plots was projected forward for different time intervals using the transition function. All possible ascending growth intervals were considered and the thinnings were simulated at the same age that was carried out using the real rate G_{removed} / G_{total}. The root mean square error and critical error in basal area estimations were calculated for different time intervals and the results are presented in Figure 4.
Figure 4
Plots of RMSE in percent (up) and critical error in percent (down) for basal area predictions against different prediction intervals. The rectangle indicates the maximum time interval for a required accuracy of 15% in basal area estimations. 
According to Huang et al. (2003) an accuracy within 10–20% of the observed mean value at 95% confidence intervals is reasonable as limit for the acceptance level. The 15% error is reached at a projection interval of 25 years and 18 years for RMSE and critical error, respectively. Therefore, we conclude that the model provides satisfactory predictions for time intervals up to 20 years.
4. CONCLUSIONS
This paper presents the most important results of a dynamic growth modeling study involving European Beech (Fagus sylvatica) forests, based on an extensive dataset provided by the Swiss Federal Institute for Forest, Snow and Landscape Research WSL, in Birmensdorf/Switzerland. The dataset includes 143 research plots, which had been remeasured between one and seventeen times, providing up to 16 interval measurements per research plot. In this model, the initial stand conditions at any point in time are defined by the three state variables dominant height, basal area and number of trees per hectare. According to this basic structure, the wholestand model requires five standlevel inputs for simulation: the stand ages at the beginning and at the end of the projection interval as well as the initial dominant height, number of trees per hectare and stand basal area.
The model uses transition functions to project the state variables at any particular time. Two of these transition functions were derived using the generalized algebraic difference approach (GADA) and the third one was derived with the algebraic difference approach (ADA). These techniques ensure that baseage and path invariance properties provide consistent predictions. All the functions were fitted simultaneously using iterative seemingly unrelated regression and a baseageinvariant method that accounts for sitespecific and global effects (Cieszewski et al., 2000).
The global wholestand growth model proved to be robust for medium term projections of basal area. Considering the required accuracy in forest growth modelling, we can state that, based on the critical errors and relative RMSE, the model provides satisfactory predictions of basal area for time intervals up to 20 years.
Several single tree growth simulators have been tested during the past years for beech stands in Switzerland. However, their applicability depends on detailed stand structure information, which is not always available (Schmid et al., 2006). The wholestand model presented in this study only requires five standlevel inputs for simulation: the age of the stand at the beginning and at the end of the projection interval as well as the initial dominant height, number of trees per hectare and stand basal area.
This model bridges the gap between the traditional yield tables and the desitable single tree models. It is preferable to yield tables in evaluating alternative silvicultural strategies in evenaged beech forests with varying densities and species combinations. One particular advantage is the fact that the weight and the type of a thinning may be specified quite easily using, for example, the basal area and stem number proportions removed in a thinning (refer to Gadow, 2006, p. 133). Moreover, the relatively simple structure of the growth model makes it suitable for decision support systems that enable forest managers to generate optimal management strategies and developing a reasonable design of a forested landscape in space and time (Gadow and Pukkala, 2008).
Acknowledgments
We are grateful to be able to use the extensive data base provided by the Swiss Federal Institute for Forest, Snow and Landscape Research.
References
 ÁlvarezGonzález J.G., CastedoDorado F., Ruiz González A.D., López Sánchez C.A. and Gadow K.V., 2004. A twostep mortality model for evenaged stands of Pinus radiata D. Don in Galicia (Northwestern Spain). Ann. For. Sci. 61: 439–448 [Google Scholar]
 ÁlvarezGonzález J.G., Ruiz González A.D., Rodríguez Soalleiro R. and Barrio Anta M., 2005. Ecoregional site index models for Pinus pinaster in Galicia (northwestern Spain). Ann. For. Sci. 62: 115–127 [CrossRef] [EDP Sciences] [Google Scholar]
 Amateis R.L., 2000. Modeling response to thinning in loblolly pine plantations. South. J. Appl. For. 24: 17–22 [Google Scholar]
 Badoux E., 1983. Ertragstafeln Buche. Eidg. Amt. Forstl. Versuchswes., 3rd ed. [Google Scholar]
 Bailey R.L. and Clutter J.L., 1974. Baseage invariant polymorphic site curves. For. Sci. 20: 155–159 [Google Scholar]
 Barrio Anta M., Castedo Dorado F., DiéguezAranda U., Álvarez González, J.G., Parresol B.R. and Rodríguez R., 2006. Development of a basal area growth system for maritime pine in northwestern Spain using the generalized algebraic difference approach. Can. J. For. Res. 36: 1461–1474 [CrossRef] [Google Scholar]
 Bertalanffy L.V., 1957. Quantitative laws in metabolism and growth. Q. Rev. Biol. 32: 217–231 [Google Scholar]
 Bolte A., Czajkowski T. and Kompa T., 2007. The northeastern distribution range of European beech – a review. Forestry 80: 413–429 [CrossRef] [Google Scholar]
 CastedoDorado F., DiéguezAranda U. and ÁlvarezGonzález J.G., 2007. A growth model for Pinus radiata D. Don stands in northwestern Spain. Ann. For. Sci. 64: 453–465 [Google Scholar]
 CastedoDorado F., DiéguezAranda U., Barrio Anta M. and ÁlvarezGonzález J.G., 2007. Modelling stand basal area growth for radiata pine plantations in Northwestern Spain using the GADA. Ann. For. Sci. 64: 609–619 [CrossRef] [EDP Sciences] [Google Scholar]
 Cieszewski C.J., 2002. Comparing fixed and variablebaseage site equations having single versus multiple asymptotes. For. Sci. 48: 7–23 [Google Scholar]
 Cieszewski C.J., 2003. Developing a wellbehaved dynamic site equation using a modified Hossfeld IV function Y3 = (axm)/ (c + xm_1), a simplified mixedmodel and scant subalpine fir data. For. Sci. 49: 539–554 [Google Scholar]
 Cieszewski C.J. and Bailey R.L., 2000. Generalized algebraic difference approach: theory based derivation of dynamic site equations with polymorphism and variable asymptotes. For. Sci. 46: 116–126 [Google Scholar]
 Cieszewski C.J., Harrison M. and Martin S.W, 2000. Practical methods for estimating nonbiased parameters in selfreferencing growth and yield models. University of Georgia, PMRCTR 20007. [Google Scholar]
 Clutter J.L. and Jones E.P., 1980. Prediction of growth after thinning in oldfield slash pine plantations, USDA For. Serv. Pap. SE217. [Google Scholar]
 Clutter J.L., Fortson J.C., Pienaar L.V., Brister G.H. and Bailey R.L., 1983. Timber management – A quantitative approach, John Wiley & Sons, 333 p. [Google Scholar]
 DiéguezAranda U., CastedoDorado F., Álvarez González J.G. and RodríguezSoalleiro R., 2005. Modelling mortality of Scots pine (Pinus sylvestris L.) plantations in the northwest of Spain. Eur. J. For. Res. 124: 143–153 [CrossRef] [Google Scholar]
 DiéguezAranda U., Castedo F., Álvarez González J.G. and Rojo A., 2006. Dynamic growth model for Scots pine (Pinus sylvestris L.) plantations in Galicia (northwestern Spain). Ecol. Model. 191: 225–242 [CrossRef] [Google Scholar]
 Gadow K.V., 2006. Forsteinrichtung – Adaptive Steuerung und Mehrpfadprinzip. Universitätsdrucke Göttingen. [Google Scholar]
 Gadow K.V. and Pukkala T., 2008. Designing Green Landscapes. Managing Forest Ecosystems Vol. 15, Springer Verlag, Dordrecht. [Google Scholar]
 García O., 1988. Growth modelling – a (re)view. N. Z. J. For. Sci. 33 (3): 14–17. [Google Scholar]
 García O., 1994. The StateSpace Approach in Growth Modeling. Can. J. For. Res. 24: 1894–1903 [CrossRef] [Google Scholar]
 Gregoire T.G., Schabenberger O. and Barrett J.P., 1995. Linear modelling of irregularly spaced, unbalanced, longitudinal data from permanentplot measurements. Can. J. For. Res. 25: 137–156. [CrossRef] [Google Scholar]
 Hallenbarter D., Hasenauer H. and Zingg A., 2005. Validierung des Waldwachstumsmodells MOSES für Schweizer Wälder. Schweiz. Z. Forstwes. 156, 5: 149–156. [CrossRef] [Google Scholar]
 Hasenauer H., Burkhart H.E. and Amateis R.L., 1997. Basal area development in thinned and unthinned loblolly pine plantations. Can. J. For. Res. 27: 265–271 [CrossRef] [Google Scholar]
 Huang S., Yang Y. and Wang Y., 2003. A critical look at procedures for validating growth and yield models. In: Amaro A., Reed D. and Soares P. (Eds.). Modelling forest systems. CAB International, Wallingford, UK, pp. 271–293. [Google Scholar]
 Hynynen J., 1995. Predicting the growth response to thinning for Scots pine stands using individualtree growth models. Silva Fenn. 29: 225–247 [Google Scholar]
 Knoebel B.R., Burkhart H.E. and Beck D.E., 1986. A growth and yield model for thinned stands of yellowpoplar. For. Sci. Monogr. 27. [Google Scholar]
 Korf V., 1939. Pøíspìvek k matematické definici vzrùstového zákona lesních porostù. Lesnická práce 18: 339–356 [Google Scholar]
 Kozak A. and Kozak R., 2003. Does cross validation provide additional information in the evaluation of regression models? Can J. For. Res. 33: 976987 [CrossRef] [Google Scholar]
 Lindstrom M. and Bates D., 1990. Nonlinear mixed effects models for repeated measures data. Biometrics 46: 673–687 [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
 Monserud R.A., 1984. Height growth and siteindex curves for inland Douglasfir based on stem analysis data and forest habitat type. For. Sci. 30: 943–965 [Google Scholar]
 Nagel J., Albert M. and Schmidt M., 2002. Das waldbauliche Prognose und Entscheidungsmodell BWINPro 6.1. Forst Holz 57 (15/16): 486–493. [Google Scholar]
 NordLarsen T. and Johannsen V.K., 2007. A statespace approach to stand growth modelling of European beech. Ann. For. Sci. 64: 365–374 [CrossRef] [EDP Sciences] [Google Scholar]
 Paulsen J.C., 1795. Praktische Anweisung zum Forstwesen, Detmold. [Google Scholar]
 Pienaar L.V. and Shiver B.D., 1984. An analysis and models of basal area growth in 45yearold unthinned and thinned slash pine plantation plots. For. Sci. 30: 933–942 [Google Scholar]
 Pienaar L.V., Shiver B.D. and Grider G.E., 1985. Predicting basal area growth in thinned slash pine plantations. For. Sci. 31: 731–741 [Google Scholar]
 Pretzsch H., Biber P., Ïurský J., Gadow K.V., Hasenauer H., Kändler G., Kenk G., Kublin E., Nagel J., Pukkala T., Skovsgaard J.P., Sodtke R. and Sterba, H., 2002. Recommendations for standardized documentation and further development of forest growth simulators. Forstw. Cbl. 121 (3): 138–151. [CrossRef] [Google Scholar]
 Rennolls K. and Peace A., 1986. Flow models of mortality and yield for unthinned forest stands. Forestry 59: 47–58 [CrossRef] [Google Scholar]
 Reynolds M.R., 1984. Estimating the error in model predictions. For. Sci. 30: 454–468 [Google Scholar]
 SAEFL/WSL, 2005. Forest Report 2005 – Facts and Figures about the Condition of Swiss Forests. Swiss Agency for the Environment, Forest and Landscape and Swiss Federal Research Institute, Berne/Birmensdorf. [Google Scholar]
 SAS Institute Inc., 2004. SAS/ETS1 9.1.2. User’s Guide. SAS Institute Inc., Cary, NC. [Google Scholar]
 Schmid S., Zingg A., Biber, P. and Bugmann H., 2006. Evaluation of the forest growth model SILVA along an elevational gradient in Switzerland. Eur. J. For. Res. 125: 43–55 [CrossRef] [Google Scholar]
 Sharma M., Smith M., Burkhart H.E. and Amateis R.L., 2006. Modeling the impact of thinning on height development of dominant and codominant trees, Ann. For. Sci. 63: 349–354 [CrossRef] [EDP Sciences] [Google Scholar]
 Sterba H. and Monserud R.A., 1997. Applicability of the forest stand growth simulator Prognaus for the Austrian part of the Bohemian Massif. Ecol. Model. 98: 23–34 [CrossRef] [Google Scholar]
 Vanclay J.K., 1995. Growth models for tropical forests: a synthesis of models and methods. For. Sci. 41: 7–42 [Google Scholar]
 Woollons R.C., 1998. Evenaged stand mortality estimation through a twostep regression process, For. Ecol. Manage. 105: 189–195 [CrossRef] [Google Scholar]
 Yang Y., Monserud R.A. and Huang S., 2004. An evaluation of diagnostic tests and their roles in validating forest biometrics models. Can. J. For. Res. 34: 619–629 [CrossRef] [Google Scholar]
 Zarnoch S.J., Feduccia D.P., Baldwin V.C. and Dell T.R., 1991. Growth and yield predictions for thinned and unthinned slash pine plantations on cutover sites in the West Gulf region. USDA Forest Service Res. Pap. SO264. [Google Scholar]
All Tables
Mathematical expression of some mortality models derived from differential equations (8) to (10). Where N and N_{0} are the number of trees per hectare present at ages t and t_{0}, respectively.
Parameter estimates and goodnessoffit statistics for the transition functions for dominant height considered.
Parameter estimates and goodnessoffit statistics for the transition functions for basal area considered.
Parameter estimates and goodnessoffit statistics for the transition functions for reduction in tree number considered.
All Figures
Figure 1
Map of all beech growth and yield plots (black dots = beech plots with tree height measurements, white dots = plots without height measurements). Map reproduced with permission of swisstopo (JA082265). 

In the text 
Figure 2
HResiduals versus HLag1Residuals for model E1 (first row), GResiduals versus GLag1Residuals for model E1 (second row) and NResiduals versus NLag1Residuals for model E5 (third row) fitted using a firstorder continuous time autoregressive error structure. 

In the text 
Figure 3
Upper left: dominant height curves for site indices of 18, 22, 26, 30 and 34 m at a reference age of 80 years for the GADA equation E1. Upper right: reduction in tree number curves for 400, 300, 2 200, 3 100, 4 000 trees per ha at a reference age of 40 years for the model (E5). Down: basal area growth rates for thinned (left) and unthinned stands (right) as affected by initial stem number and age overlaid on the trajectories of basal area. 

In the text 
Figure 4
Plots of RMSE in percent (up) and critical error in percent (down) for basal area predictions against different prediction intervals. The rectangle indicates the maximum time interval for a required accuracy of 15% in basal area estimations. 

In the text 