Modeling top height growth of red alder plantations

Forest Ecology and Management 258 (2009) 323–331

Contents lists available at ScienceDirect

Forest Ecology and Management journal homepage: www.elsevier.com/locate/foreco

Modeling top height growth of red alder plantations Aaron R. Weiskittel a,*, David W. Hann b, David E. Hibbs c, Tzeng Yih Lam b, Andrew A. Bluhm c a

School of Forest Resources, University of Maine, 5755 Nutting Hall, Orono, ME 04469, United States Department of Forest Engineering, Resources, and Management, 204 Peavy Hall, Oregon State University, Corvallis, OR 97331, United States c Department of Forest Ecosystems and Society, 313 Richardson Hall, Oregon State University, Corvallis, OR 97331, United States b

A R T I C L E I N F O

A B S T R A C T

Article history: Received 17 November 2008 Received in revised form 9 April 2009 Accepted 26 April 2009

Height growth equations for dominant trees are needed for growth and yield projections, to determine appropriate silvicultural regimes, and to estimate site index. Red alder [Alnus rubra Bong.] is a fastgrowing hardwood species that is widely planted in the Pacific Northwest, USA. However, red alder dominant height growth equations used currently have been determined using stem analysis trees from natural stands rather than repeated measurements of stand-level top height from plantations, which may cause them to be biased. A regional dataset of red alder plantations was complied and used to construct a dynamic base-age invariant top height growth equation. Ten anamorphic and polymorphic Generalized Algebraic Difference Approach (GADA) forms were fit using the forward difference approach. The Chapman–Richards anamorphic and Schumacher anamorphic model forms were the only ones with statistically significant parameters that yielded biologically reasonable predictions across a full range of the available data. The Schumacher model form performed better on three independent datasets and, therefore, was selected as the final model. The resulting top height growth equations differed appreciably from tree-level dominant height growth equations developed using data from natural stands, particularly at the younger ages and on lower site indices. Both the rate and shape parameters of the Schumacher function were not influenced by initial planting density. However, this analysis indicates that the asymptote, which is related to site index, may be reduced for plantations with initial planting density below 500 trees ha1. The final equation can be used for predictions of top height (and thus) site index for red alder plantations across a range of different growing conditions. ß 2009 Elsevier B.V. All rights reserved.

Keywords: Top height Dominant height Red alder Plantations Planting density Pacific Northwest Generalized Algebraic Difference Approach Site index

1. Introduction Although varying definitions of stand dominant height exist, most measures of it are well correlated with site potential productivity (Sharma et al., 2002) and thus, the resulting equations can be used to estimate site index. The three elements of modeling stand dominant height are the choice of a model form, the structure of the modeling data, and the statistical procedures used to estimate model parameters. Dominant height models have taken many different forms and significant differences in performance can exist between them (e.g., Krumland and Eng, 2005; Cieszewski and Strub, 2008). The modeling data can also come from a variety of sources including measurements of top height taken on temporary plots, from stem analysis of individual dominant trees, or from repeated measurements of top height taken on permanent plots. These different methods can have important implications on the resulting site index equation (Raulier et al., 2003). Given a basic model form and a data set, the guide curve method, the parameter prediction

* Corresponding author. Tel.: +1 207 581 2857; fax: +1 207 581 2875. E-mail address: [email protected] (A.R. Weiskittel). 0378-1127/$ – see front matter ß 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.foreco.2009.04.029

method, or the difference equation method can be used to construct a dominant height growth equation (Clutter et al., 1983; Wang et al., 2008a) through the application of alternative statistical procedures such as fixed parameter, dummy variable, or mixed parameter nonlinear regression. The best statistical method for modeling dominant height growth will most likely be species dependent and is still debated in the literature (Cieszewski and Strub, 2007; Wang et al., 2007a, 2008b). A form of the difference equation method, Generalized Algebraic Difference Approach (GADA; Cieszewski and Bailey, 2000), has become a common means for fitting dominant height growth equations. GADA (for simplicity, we include the simpler Algebraic Difference Approach equations in this class of equations) dominant height growth equations have several desirable properties such as base-age invariance, generally fewer parameters to be estimated, prediction of dominant height forward or backward in age given any initial conditions, and the ability to identify site dependent curve changes from mean height-age trends (Cieszewski and Bailey, 2000). GADA dominant height growth equations have been developed for several conifer species including interior Douglas-fir [Pseudotsuga menziesii var. glauca (Mirb.) Franco] (Cieszewski, 2001),

324

A.R. Weiskittel et al. / Forest Ecology and Management 258 (2009) 323–331

loblolly pine [Pinus taeda L.] (Wang et al., 2008a), radiata pine [Pinus radiata D. Don] (Die´guez-Aranda et al., 2005), subalpine fir [Abies lasiocarpa (Hook.) Nutt.] (Cieszewski, 2003), Scots pine [Pinus sylvestris L.] (Cieszewski et al., 2007), and 10 northern California conifer species or species groups (Krumland and Eng, 2005). This approach has also been used to develop dominant height growth equations for hardwood species such as European beech [Fagus sylvatica L.] (Nord-Larsen, 2006), Pyrenean oak [Quercus pyrenaica Willd.] (Carvalho and Parresol, 2005), Spanish cork oak [Quercus suber L.] (Sa´nchez-Gonza´lez et al., 2008), and six species or species groups in northern California (Krumland and Eng, 2005). The application of GADA to fast-growing hardwood species has been limited to Eucalyptus species (Amaro et al., 1998; Calegario et al., 2005; Wang et al., 2007b) and to naturally regenerated red alder [Alnus rubra Bong.] trees (Krumland and Eng, 2005). Development of a dynamic dominant height growth equation for fast-growing plantation hardwood species may be particularly problematic because of their rapid peak in height growth, potentially greater sensitivity of height growth to stand density, and higher variability in annual height growth patterns. The Alnus genus is comprised of approximately 30 species that are widely distributed throughout the northern temperate zone and several species are of high commercial importance including red, black [Alnus glutinosa (L.) Gaertn.], and Italian [Alnus cordata (Loisel) Duby] alders. It is our observation that most members of the Betulaceae family, such as the Alnus species, exhibit common properties of shade intolerance, relatively short lifespans, and fast early growth. In the Pacific Northwest (PNW) region of the United States, planting of red alder has significantly increased in recent years and this trend is expected to continue (Briggs, 2007). On high quality sites, red alder plantations can reach a mean diameter of 30 cm in approximately 22 years and even greater productivities have been reported (Bluhm and Hibbs, 2006). Until about age 30, red alder height growth greatly exceeds the height growth of the other native conifers in the region (Harrington et al., 1994). Several tree-level dominant height growth equations exist for natural, unmanaged red alder stands including those of Porter and Wiant (1965), Worthington et al. (1960), Harrington and Curtis (1986), Nigh and Courtin (1998), and Krumland and Eng (2005). However, a tree-level dominant height growth or a stand-level top height growth equation for managed red alder plantations does not exist and is much needed in the region. In addition, all existing dominant height growth equations for red alder were developed using treelevel stem analysis data as opposed to repeated measurements of stand top height, which can significantly influence the behavior of the site curves. For example, stand-level top height growth equations often predict taller heights below the base age and shorter heights above the base age when compared to tree-level dominant height growth equations developed from stem analysis (Smith, 1984; Monserud, 1985; Raulier et al., 2003). Raulier et al. (2003) attributed this difference to changes over time in the composition of the dominant trees. Stem analysis studies, such as Harrington and Curtis (1986) and Nigh and Courtin (1998), require actual measurements of site index in their analysis so only dominant trees with ages equal to or greater than the base age are sampled. As a result, their measured heights at younger ages may not be characteristic of the height of true dominants found on the plot at those younger ages due to changes in social status (Oliver and Larson, 1996; Raulier et al., 2003). However, the result of this dynamic behavior would be adequately characterized by using repeated measurements of top height on permanent plots (Raulier et al., 2003). A key assumption in using dominant height or top height as an index of site productivity potential is that height growth of dominant trees is not influenced by stand density. Several studies have shown this to be generally the case (Skovsgaard and Vanclay,

2008). However, higher levels of stand density has been found to reduce dominant height growth for a number of predominantly intolerant conifer species with strong epinastic control (Gaiser and Merz, 1951; Lynch, 1958; Cieszewski and Bella, 1993; Scott et al., 1998; MacFarlane et al., 2000; Flewelling et al., 2001). As a result, stand density has been incorporated into dominant height growth equations for ponderosa pine (Lynch, 1958), lodgepole pine (Cieszewski and Bella, 1993), and Douglas-fir (Flewelling et al., 2001). Less attention has been given to hardwood species with weak epinastic control that can reduce dominant height growth at low levels of stand density. In red alder, DeBell and Giordano (1994) reported that the average height growth rate was lower for trees growing on plots with densities below 700 trees ha1. They also reported that average height growth rate was reduced on plots with a density of over 4000 trees ha1 (DeBell and Giordano, 1994). Maximum height growth of red alder at intermediate densities has also been reported in several additional studies (Bormann and Gordon, 1984; Cole and Newton, 1987; Giordano and Hibbs, 1993). All of these studies were conducted using very young trees on a very limited number of installations. The overall goal of this research project was to model top height growth of plantation grown red alder using an extensive, regional database. Specific objectives were to: (1) fit several different GADA top height growth model forms using three alternative statistical techniques; (2) compare the predictive performance of the resulting equations on both the fitting data set and an evaluation data set; (3) for the final model, compare predicted development of top height in plantations to the development of dominant height in natural stands; (4) for the final model, assess the influence of planting density on estimated top height growth. 2. Methods 2.1. Study area Sixty existing research installations of the Hardwood Silviculture Cooperative (HSC; http://www.cof.orst.edu/coops/hsc) and the Weyerhaeuser Company were used in this study. The installations were located in Oregon, Washington, and British Columbia and covered the full range of growing conditions for the species, including sites where commercial plantation management would not be implemented (Table 1). All installations were located between the Pacific Coast and the Cascade Mountains. The overall climate for the region is humid oceanic, with a distinct dry summer and a cool, wet winter. Regional 25-year annual mean rainfall ranged from 125 to 354 cm (http://www.daymet.org). The mean minimum temperature in January ranged from 8.3 to 2.9 8C and the mean maximum temperature in August ranged from 17.9 to 25.9 8C (http://www.daymet.org). Elevation ranged from 38 to 550 m above sea level and all aspects were represented. Soils varied from shallow sandy loams to very deep silty clay loams. The Table 1 Attributes of the research installations used in this analysis (n = 60). Attribute Mean minimum temperature in January (8C) Mean maximum temperature in August (8C) Growing degree days (>5 8C) Annual precipitation (cm) Elevation (m) Slope (%) Aspect (8) Soil water holding capacity (mm)

Mean

Minimum

Maximum

1.8

8.3

2.9

22.9

1.8

17.9

25.9

4185.8 214.1 238.5 16.7 167.1 331.9

556.5 48.3 149.8 17.8 101.0 152.6

1885.1 125.1 38.0 0.0 0.0 27.5

4864.7 354.1 550.0 90.0 360.0 611.0

0.49

S.D.

A.R. Weiskittel et al. / Forest Ecology and Management 258 (2009) 323–331 Table 2 Measurement attributes for the modeling (n = 797) and evaluation (n = 790) plots used in the analysis. Attribute

S.D.

Minimum

Maximum

Modeling plot/measurements (n = 797) Top height (m) 10.9 Total age (years) 8.9 Planting density (# ha1) 1346.7 Trees per ha 1278.1 Basal area (m2 ha1) 10.1

Mean

4.6 4.0 889.7 826.2 8.2

1.9 3.0 168.0 138.4 0.1

24.6 18.0 3765.8 3743.6 30.9

Evaluation plot/measurements (n = 790) Top height (m) 9.7 Total age (years) 7.5 Planting density (# ha1) 1719.7 Trees per ha 1200.1 Basal area (m2 ha1) 6.2

4.0 3.4 813.6 744.5 5.4

2.1 3.0 370.6 252.0 0.1

20.4 17.0 3654.6 3617.5 27.0

plantations were established between 1989 and 1997 at varying planting densities and levels of vegetation control. 2.2. Data Several permanent measurement plots were established at each installation (Tables 2 and 3). Size of the measurement plots ranged from 0.097 to 0.202 ha, excluding their buffers. The HSC installations included four different blocks with initial targeted planting densities of 248, 569, 1298, and 2966 trees ha1. Each block was subdivided into several randomly assigned treatment plots including a control as well as a range of thinning and pruning options. Individual tree-level diameter and heights were first measured at age 4 and have been remeasured every 3–5 years since. The Weyerhaeuser Company installations used a similar design as they had blocks with initial, targeted planting densities of 494, 988, 1730, and 2965 trees ha1. Individual tree-level diameter and heights were first measured between ages 3 and 5 and have been remeasured every 1–5 years since. All trees on each plot were measured for diameter at 1.37 m (DBH; to the nearest 0.1 cm) and a subsample of trees were measured for total height. Missing tree heights were imputed using a regional height-diameter equation developed for plantation grown red alder and then scaled to each plot’s height measurements using simple linear regression through the origin similar to the approach of Hann and Hanus (2004). Top height (HD; m) at each age was then calculated as the mean height of the 98.8 largest trees ha1 (based on diameter) that did not have height damage or other severe damage for each measurement period. Stand age (A) was defined as the number of years since seed. Because all stands were planted with 1-year-old seedlings, A was determined by adding one to the known plantation age.

Table 3 Number of plot and measurement combinations by stand age and Harrington and Curtis (1986) natural stand site index for the modeling and evaluation plots. Stand age

Site index (m) 10–15

16–20

21–25

>25

Modeling plot/measurements (n = 797) 0–5 12 156 6–10 16 151 11–15 11 90 16–20 2 21

39 149 126 16

– 1 6 1

Evaluation plot/measurements (n = 790) 0–5 8 225 6–10 2 189 11–15 6 84 16–20 2 3

40 137 91 2

– 1 – –

325

The data was divided into a modeling data set and an evaluation data set. The modeling data set consisted of all control plot data. The evaluation data set consisted of all measurements taken on the treatment plots before the treatments were applied. The modeling data set was restricted to stand ages between 3 and 18 years. Because red alder has been known to live to 100 years or longer (Harrington et al., 1994), examination of the extrapolative behavior of alternative models to older ages was important in deciding a final model. Two other HSC datasets, one with very young trees and one with older trees, were used to evaluate model behavior. In the first dataset, plot-level tree heights were measured on these installations (153 total plots) at ages of 1, 2, and 3 years. The age one data was ignored in this analysis because that was the seedling size at planting and, therefore, differences in site productivity were not manifested in the measurements. The HSC also had data older than 18 years from four installations in natural stands of pure red alder. Plot size ranged from 0.101 to 0.202-ha, excluding their buffers. Each of the control plots had five measurements of DBH and height with initial stand ages ranging from 15 to 18 years, and final stand ages ranging from 29 to 32 years. Forty heights were subsampled across the diameter range on each of the plots. Again, missing tree heights were imputed using the developed regional heightdiameter equation that was scaled to each plot’s height measurements using simple linear regression through the origin. Although the older data is from natural stands, we wanted to ensure the equations extrapolated in a biologically realistic fashion and this was the only data available to achieve this aim. 2.3. Data analysis Ten different GADA type model forms were evaluated as the first step of the analysis: (1) Chapman–Richards anamorphic, (2) Chapman–Richards polymorphic, (3) Schumacher anamorphic, (4) Schumacher polymorphic, (5) log-logistic, (6) King–Prodan (King, 1966), (7) McDill and Amateis (1992), (8) Cieszewski (2001), (9) Cieszewski (2003) and (10) Cieszewski et al. (2006). Formulations for model forms 1–6 were from Krumland and Eng (2005). The equations were fit under the assumption that all parameters are fixed using generalized nonlinear regression with the gnls package in R (R Development Core Team 2007). The continuous first-order autoregressive (CAR1) procedure was used to account for residual autocorrelation, while heteroscedasticity was addressed using a variance power function of total age (A). Nonlinear regression with the non-overlapping forward difference data structure (NLDIF) was used in this initial screening of alternative model forms because it has been previously found to perform well (Wang et al., 2007a). The resulting parameterizations of the equations had to meet the following criteria for further consideration: 1. All parameters had to be significantly different from zero with a probability of 95% (p = 0.95). 2. Predictions from the equation had to exhibit a sigmoidal form over the entire range of possible ages and site indices. 3. Measured top height (HDM; m) had to equal predicted top height (HDP) when the age of measurement (AM) equaled the age of prediction (AP) (Bailey and Clutter, 1974). Models meeting these criteria were then fit using both the nonlinear dummy variable (NLDV) estimation method and the nonlinear mixed effects (NLME) estimation methods, and the data structure that contained all possible pairs of HDM–AM measurements (Wang et al., 2008a). NLME was fit using the NLME package in R, while the NLDV parameters were estimated with SAS. Although the use of R2 in NLME has been debated, a generalized form of it was used in this analysis (Xu, 2003). This measure

326

A.R. Weiskittel et al. / Forest Ecology and Management 258 (2009) 323–331

effectively partitioned the variance explained by the fixed effects and the random effects at each level. Since plots were nested within installations, multi-level random effects were estimated for installations and plots within an installation. Preliminary analysis indicated that most of the variation in the parameters was between installations-rather than within. These alternatives and NLDIF were then evaluated based upon: 1. Their fit to the modeling data using the coefficient of determination (R2), Akaike’s information criteria (AIC), and the residual standard error (RSE) of the weighted residuals. The weighted residuals were also examined over A, HDP, and predicted site index (SI as defined as HDP at A of 20) using both graphs and lack-of-fit statistics on different strata of the data. 2. Their ability to accurately predict the evaluation data, by specific ages and overall, using both the mean of the unweighted differences (where the difference is defined as HDP–HDM) as a measure of bias (i.e., accuracy) and the standard deviation of the unweighted differences, which removes the bias (as a measure of precision). To conduct this analysis, HDM and A from the last measurement on each plot were used to localize the equation to the plot by applying the least squares procedure described in Wang et al. (2008a) to the NLDIF, NLDV, and NLME parameterizations for each of the remaining candidate model forms (we will abbreviate these combinations as NLDIF-LS, NLDV-LS, and NLME-LS, respectively). In addition, the Wang et al. (2008a) mixed model BLUP estimator for localizing the NLME estimator (NLME-ME) was also applied to each remaining candidate model form. The resulting localized equations were then used to calculate HDP for each A on the plot. 3. The parameter estimation and localization method judged superior from steps 1 and 2 was then used in each of the remaining candidate model forms to evaluate the predictive behavior when extrapolated to A values not found in the modeling data set. This analysis used the younger and older evaluation data sets not utilized in modeling or evaluation. The analysis procedures were the same as used in step 2 with the exception that the first measurement of HDM and A was used to localize the equation to the older, natural plot data, and the measurement of HDM at A of 2 was used to localize the equation to the younger plantation data.

and upper confidence bands that fall below zero indicates the presence of under prediction bias. The greater the distance between these particular confidence bands and zero, the greater the severity of the bias. The final GADA model form (and associated parameters) was also used to explore the possible impacts of stand density using the following approach: 1. The graph of the weighted residuals (from the fitting process) over stand density was examined for trends. 2. The graph of plot-level relative SI over planting density was also examined for trends. To calculate relative SI for a plot, SI was calculated for each plot on an installation, the installation’s largest value of SI was identified, and then all of the plot SI values on the installation were divided by this largest value. The rationale for this graph recognizes that there often is variation in SI between plots on an installation due to differences in microsite and random events. As a result, any plot in an installation may exhibit the largest SI value for the installation, absent a treatment effect. If there is no impact of planting density upon predicted SI for the 60 installations, then the graph of relative density should exhibit values at or near one across all planting densities. Deviation from this expectation could indicate that planting density has an impact upon predictions from the base model. For a given planting density, the greater the deviation of the largest relative SI value is from one, the greater the implied impact of planting density upon predicted SI. 3. If trends were detected, then appropriate indicator variables based on stand density were formulated to statistically test whether the trend(s) indicated a significant impact upon the rate and shape parameters of the base model. Because the GADA approaches used in this study replaces the scale parameter, SI in this case, with actual measurements of A and H40, we could not statistically test whether any trends detected in step 1 had a significant impact upon scale parameter (i.e., predicted SI). 4. If the trend(s) was significant, then appropriate modifications to the rate and/or shape parameter(s) of the base model were formulated to account for the impact of stand density upon HDP over time. 3. Results

The evaluation data set and the following process was then used to compare the final GADA model form (and associated parameters) to the Harrington and Curtis (1986) system of equations (consisting of a dominant height growth rate equation and a site index equation): 1. The final GADA equation was localized to each plot using the least squares procedure described in Wang et al. (2008a) and the HDM and A values from the last measurement on each plot. For the Harrington and Curtis (1986) equations, the localization consisted of calculating site index using the HDM and A values from the last measurement on each plot and the Harrington and Curtis (1986) site index equation. 2. The resulting localized equations were then used to calculate HDP for each A on the plot. The dominant height growth rate equation of Harrington and Curtis (1986) was used in this step of the analysis. 3. For each equation and each value of A in the evaluation data set, the mean of the unweighted differences (where the difference is defined as HDP  HDM) was calculated as a measure of bias (i.e., accuracy) and the 99% confidence interval for the mean difference was also calculated. It is desirable that the confidence interval include zero for each value of A. Lower confidence bands that exceed zero indicates the presence of over prediction bias,

Only two model formulations fit using the NLDIF method had statistically significant parameters and expected behavior. These included the Chapman–Richards (Eq. (1)) and Schumacher (Eq. (2)) anamorphic model formulations: HDP ¼ HDM



1  eb1 AP 1  eb1 AM b2 b AM2 Þ

HDP ¼ HDM eb1 ðAP

b2

(1)

(2)

These two models were then fit using the NLDV and the NLME procedures described in Wang et al. (2008a). For both formulations/equations, the resulting parameter estimates and standard errors are similar across the three alternative fitting methods (Table 4). The absolute values of the NLME parameters were consistently larger than the NLDIF estimates, which were generally larger than the NLDV estimates. The fit statistics indicated that the NLDV method generally provided the ‘‘best’’ statistics for all three criteria (with values of RMSE and AIC closer to zero being better and values of R2 closer to one being better). The NLME method had the poorest statistics for RMSE and AIC, while the NLDIF method gave the poorest R2 statistics. The fit statistics also indicate a mixed view of the two

A.R. Weiskittel et al. / Forest Ecology and Management 258 (2009) 323–331

327

Table 4 Parameter estimates, standard errors, and fit statistics [residual stand error (RSE, m), Akaike’s information criteria (AIC), and R-square (R2)] for the Chapman–Richards anamorphic and Schumacher anamorphic equations by fitting technique. Parameter

NLDIF

NLDV

NLME

Estimate

S.E.

Estimate

S.E.

Estimate

S.E.

Chapman–Richards b1 b2 RSE (m) AIC R2

0.14335 1.49193 3.2489 1919.067 0.8715

0.00663 0.05208 – – –

0.15253 1.60800 2.8894 1543.324 0.9620

0.00781 0.06501 – – –

0.16020 1.70087 4.7928 4344.233 0.9422

0.00571 0.05469 – – –

Schumacher b1 b2 RSE (m) AIC R2

4.48127 0.65888 3.1336 1941.234 0.8665

0.09841 0.03443 – – –

4.32340 0.59956 1.8582 1245.324 0.9873

0.04545 0.02973 – – –

5.42200 0.86032 5.0706 4379.991 0.9409

0.16488 0.02987 – – –

model forms. The Chapman–Richards formulation had the best fit statistics for the NLME method, while the Schumacher formulation had the best fit statistics for the NLDV method. The NLDIF method had mixed fit statistics for the two models. In general, the fit statistics between the two models were quite similar in size for a given fitting method. The NLDIF-LS and the NLDV-LS methods had very similar bias and precision statistics, with the average of the statistics across the two model forms being slightly better for the NLDIF-LS method (Table 5). Of the two remaining methods, NLME-LS had higher precision and lower bias than NLME-ME for both equations. Evaluation of the fitting-localization methods across A was done by examining the evaluation statistics at those ages with 50 or more observations (i.e., A values of 4, 7, 10, and 12 years), ranking the bias and precision statistics for each A, and summing the rankings for each fitting-localization combination. This showed that the NLDIF-LS method performed the best across all A values, followed by NLDV-LS, NLME-LS, and NLME-ME. From these result, NLDIF-LS parameters were used for the remaining analysis. Examination of the NLDIF-LS evaluation statistics indicated that the Schumacher model had slightly lower bias and the same precision as the Chapman–Richards model. Graphs of HDP over age for three SI values (10, 20, and 30 m at 20 years) for the two models

are presented in Fig. 1). Graphs of HDP PAI for the two models are also presented in Fig. 1. For comparative purposes, predictions from the tree-level dominant height growth equation of Harrington and Curtis (1986) are also presented in these graphs. These graphs illustrate differences in predictions from the two models at ages that are outside the range of the modeling data. The Schumacher model had higher peak PAI values than the Chapman–Richards model and also predicted higher PAI values for ages above 18 years (Fig. 1). The result of the latter dissimilarity is the more strongly asymptotic behavior of the Chapman–Richards model than the Schumacher model (Fig. 1). The two models were applied to the younger, independent data set by localizing them with the 2-year-old measurements and then predicting the 3-year-old measurements. Both models underpredicted HD; the Chapman–Richards model by an average of 0.72 m and the Schumacher model by an average of 0.34 m. Precision of the predictions was similar for the two models, with the standard deviation of the difference being 0.53 m for the Chapman–Richards and 0.56 m for the Schumacher model. The results of applying the two models to the older independent data set are found in Table 6. The models were localized with the first measurements on the four plots (A of 15 for three of the plots and A of 19 for the other plot) and then used to predict HDM for all of the measurements on each plot. Overall, the Chapman–Richards

Table 5 Unweighted mean bias (HDP  HDM) and unweighted standard deviation of the differences for the evaluation dataset by total stand age for the Chapman–Richards anamorphic (CR) and Schumacher anamorphic (SR) equations fit using NLDIF-LS, NLDV-LS, NLME-LS, and NLME-ME. Total stand age (years)

Number of observations

3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 All

26 231 16 48 168 34 19 60 48 68 39 15 11 4 3 790

NLDIF-LS

NLDV-LS

NLME-LS

Mean bias (m)

Bias S.D. (m)

NLME-ME

Mean bias (m)

Bias S.D. (m)

Mean bias (m)

Bias S.D. (m)

Mean bias (m)

Bias S.D. (m)

CR

SR

CR

SR

CR

SR

CR

SR

CR

SR

CR

SR

CR

SR

CR

SR

0.80 0.00 0.18 0.41 0.27 0.05 0.19 0.22 0.01 0.05 0.03 0.01 0.01 0.01 0.00 0.13

0.95 0.07 0.08 0.35 0.16 0.01 0.17 0.24 0.03 0.03 0.02 0.00 0.00 0.00 0.00 0.09

1.28 1.10 0.72 0.93 0.82 0.59 0.42 0.55 0.22 0.20 0.13 0.01 0.01 0.02 0.00 0.83

1.25 1.09 0.74 0.92 0.80 0.57 0.44 0.55 0.23 0.16 0.08 0.01 0.01 0.02 0.00 0.83

1.00 0.18 0.29 0.47 0.34 0.10 0.20 0.23 0.02 0.06 0.03 0.00 0.00 0.00 0.00 0.21

0.80 0.16 0.03 0.34 0.17 0.00 0.18 0.25 0.03 0.02 0.01 0.01 0.01 0.01 0.01 0.06

1.27 1.11 0.71 0.94 0.82 0.58 0.42 0.56 0.22 0.20 0.13 0.02 0.02 0.00 0.00 0.84

1.27 1.09 0.73 0.92 0.80 0.57 0.44 0.55 0.23 0.16 0.07 0.02 0.02 0.02 0.02 0.83

1.14 0.30 0.37 0.51 0.37 0.12 0.21 0.22 0.01 0.06 0.03 0.01 0.01 0.00 0.01 0.26

1.54 0.34 0.28 0.43 0.20 0.03 0.17 0.20 0.01 0.05 0.02 0.00 0.01 0.00 0.00 0.24

1.26 1.11 0.72 0.96 0.83 0.59 0.42 0.56 0.22 0.20 0.13 0.01 0.01 0.02 0.02 0.85

1.21 1.10 0.73 0.92 0.81 0.59 0.42 0.54 0.22 0.19 0.11 0.01 0.01 0.02 0.00 0.85

1.24 0.20 0.51 0.07 0.29 0.07 0.35 0.03 0.40 0.07 0.18 0.05 0.02 0.31 1.74 0.19

1.63 0.25 0.40 0.09 0.12 0.06 0.33 0.04 0.48 0.01 0.13 0.11 0.07 0.34 1.89 0.17

1.09 1.25 0.86 1.12 1.12 0.89 0.66 0.91 0.86 0.73 0.94 1.01 1.01 0.74 0.73 1.09

1.07 1.26 0.89 1.12 1.12 0.92 0.70 0.95 0.93 0.79 1.02 1.09 1.10 0.79 0.79 1.13

A.R. Weiskittel et al. / Forest Ecology and Management 258 (2009) 323–331

328

Fig. 1. Plots of the GADA formulations of the Chapman–Richards anamorphic (left) and Schumacher anamorphic (right) stand-level top height (HD) growth equations developed using nonlinear regression with the non-overlapping forward difference data structure (NLDIF) fitting procedures compared to the tree-level dominant height growth equation of Harrington and Curtis (1986) for site indices of 10, 20, and 30 m at base age 20.

model under-predicted HDM by 1.53 m, while the Schumacher model under-predicted HDM by 0.63 m. Precision of the predictions was 1.49 m for the Chapman–Richards model and 1.09 m for the Schumacher model. Examination of the results across the measurement ages shows that the Schumacher equation was more accurate (exhibited lower bias), while its precision was always greater. This result is in contrast to the overall statistics. Based upon the fit statistics, the evaluation statistics, and the extrapolation statistics for the younger and older data sets, we concluded that the Schumacher model fit with the NLDIF procedure was the best for predicting

Table 6 Unweighted mean bias (HDP  HDM) and unweighted standard deviation of the bias for the older independent dataset by total stand age for the Chapman–Richards anamorphic (CR) and Schumacher anamorphic (SR) equations fit using nonlinear regression with the non-overlapping forward difference data structure (NLDIF-LS). Total stand age (years)

Number of observations

Mean bias (m)

Bias S.D. (m)

CR

SR

CR

SR

18 19 21 22 24 25 28 29 32 All

3 1 3 1 3 1 1 3 1 20

0.44 0.00 0.81 2.53 1.77 2.67 4.06 2.53 4.70 1.53

0.12 0.00 0.06 2.09 0.53 1.75 2.65 0.45 2.61 0.63

0.72 0.00 1.30 0.00 0.68 0.00 0.00 0.53 0.00 1.49

0.76 0.00 1.36 0.00 0.85 0.00 0.00 0.84 0.00 1.09

HDM in red alder plantations. Therefore the final model with estimated parameters was: 0:65888 A0:65888 Þ M

HDp ¼ HDM e4:48127ðAp

(4)

The mean differences (bias) and confidence intervals about the mean differences across all values of A in the evaluation data set are shown in Fig. 2 for the final Schumacher equation and the Harrington and Curtis (1986) equations. For the Schumacher equation, many of the A specific confidence bands incorporate zero. The mean differences for those ages in which the confidence bands do not include zero are all within 0.25 m of zero, which is within the measurement accuracy of the height measurements. For the Harrington and Curtis (1986) equations, all confidence bands exclude zero, indicating the equations are biased throughout the A range of the data. Furthermore, the size of the bias increases to around 1 m as the A values get younger. These results indicate that the dominant height growth of trees in natural red alder stands differ significantly from the top height growth of red alder plantations. A graph of the residuals over planting density (Fig. 3) illustrated slight under-prediction for stands with a planting density less than 500 trees ha1 and a slight over-prediction for stands with a planting density greater than 1500 trees ha1. Furthermore, the graph of relative SI over planting density (Fig. 4) showed a trend indicating that the predicted SI for plots with planting densities under 500 trees ha1 were always lower than the largest SI value for a given installation and that the size of the difference increased as the planting density decreased. The following model form was

A.R. Weiskittel et al. / Forest Ecology and Management 258 (2009) 323–331

329

Fig. 2. Mean residual (HDP  HDM) and the corresponding 99% confidence interval about the mean resulting from predicting the top height (HD) measurements on the evaluation data set for measured values of total age from seed (A) using the final Schumacher GADA equation for red alder plantations and the site-index/dominant height growth equations of Harrington and Curtis (1986) for dominant red alder trees growing in natural stands. Both equations were localized using the last measurement of HD and A on the plot and then the localized equations were used to predict the HD for the remaining values of A on each plot.

then formulated to test whether planting density affected the rate and shape parameters of the Schumacher model: HDP ¼ HDM eðb1;0 þb1;1 IND1 þb1;2 IND2 ÞðAP

b2;0 þb2;1 IND1 þb2;2 IND2 A

M

b2;0 þb2;1 IND1 þb2;2 IND2 Þ

(5) where IND1 = 1.0 if planting density was 720 trees ha1, 0.0 otherwise; IND2 = 1.0 if planting density was 360 trees ha1, 0.0 otherwise. The parameters b1,1, b1,2, b2,1, and b2,2 were tested with a t-test under the null hypothesis that all of these parameters were zero. Results of the tests were that the null hypothesis was true for all of

Fig. 4. Plot of relative Schumacher site index over planting density for the modeling and evaluation data sets combined.

these parameters with p = 0.95. Therefore, using this method, there is no evidence that planting density significantly affects the b1 and b2 rate and shape parameters of Eq. (2). However, Fig. 4 does show that predicted site indices are reduced for planting densities below approximately 500 trees ha1. 4. Discussion

Fig. 3. Plot of residuals (HDP  HDM) from the GADA formulation of the Schumacher anamorphic equation fit using nonlinear regression with the non-overlapping forward difference data structure (NLDIF) method over initial plant density (trees ha1) with a locally weighted polynomial (LOESS) regression line.

The results of this analysis indicated that the GADA technique produced an unbiased and precise equation for predicting top height growth of red alder. In contrast to several other recent studies in conifer species (Cieszewski, 2001, 2003; Cieszewski and Strub, 2008), the anamorphic model forms performed better than the polymorphic model forms. This finding is consistent with

330

A.R. Weiskittel et al. / Forest Ecology and Management 258 (2009) 323–331

existing tree-level dominant height growth curves for red alder. Krumland and Eng (2005) tested several anamorphic and polymorphic GADA forms for red alder trees in northern California and found that the Schumacher anamorphic form provided the best fit to their stem analysis data. Nigh and Courtin (1998) found the polymorphic parameter in their tree-level dominant height growth equation was not significantly different from zero which resulted in an anamorphic final equation. Although polymorphic forms were not explicitly tested, both Porter and Wiant (1965) and Worthington et al. (1960) used anamorphic model forms for red alder tree-level dominant height growth. Interestingly, while the Harrington and Curtis (1986) tree-level dominant height growth equation is clearly polymorphic in the extreme values of site index, its behavior over the majority of the data range is closely anamorphic. The preference for an anamorphic model form might also be the result of the relatively narrow range in site index found for the species (Table 3). Similar to the findings of Wang et al. (2007a), the NLDV-LS and NLDIFF-LS methods provided very similar accuracy and precision statistics for HDP on the evaluation data set. Furthermore, the finding that the NLME-ME method provided poorer bias and precision statistics for HDP than the NLDV-LS method for the anamorphic Chapman–Richards formulation is in agreement with the findings of Wang et al. (2008a). However, that the NLME-LS method provided poorer accuracy and precision statistics than the NLDV-LS method contradicts the findings of Wang et al. (2008a). While not in total agreement with previous work, these results add to the literature exploring GADA fitting technique and should, ultimately, lead to consensus on the best modeling practices for developing dominant height growth or top height growth equations. Overall, we found that the differences between the alternative model fitting techniques were relatively small, emphasizing the importance of an appropriate initial model form. Both the Chapman–Richards anamorphic and Schumacher anamorphic model forms produced appreciably different growth curves than Harrington and Curtis (1986) (e.g., Fig. 1), and these differences are statistically significant (Fig. 2). This is likely a function of at least two factors. First, all previous red alder dominant height growth equations are from individual trees rather than stand-level top height data. Previous work has shown that height growth curves based on tree-level data can differ significantly from those based on stand-level data (e.g., Monserud, 1985). For example, Raulier et al. (2003) showed that height growth curves based on tree-level data generally gave predictions lower than those based on stand-level data below the base age and predictions greater than the stand-level predictions above the base age, a trend repeated in this study. This is likely a function of tree social status dynamics—as the status of dominant trees change with time and is not stable (Oliver and Larson, 1996; Raulier et al., 2003). In the case of red alder, the difference in predictions between the top height growth equation developed in this study and the Harrington and Curtis (1986) dominant height growth equation were relatively large, particularly at younger ages and on lower site indices. This suggests that tree dominance is relatively instable for red alder grown in plantations. Second, the top height growth curves developed in this study are based solely on plantation data, while the previously developed equations only used data from naturally regenerated stands. The Harrington and Curtis (1986) dominant height growth equation predicted a peak in dominant height growth sooner than either the Chapman–Richards or the Schumacher top height growth equations developed in this study (Fig. 1). In addition, the Harrington and Curtis (1986) dominant height growth equation had a lower peak in dominant height growth than the Schumacher top height growth equation developed here. The results of our evaluation using the younger independent data set indicates that the

Schumacher’s higher and later peak PAI is more congruent with the top height data than the Harrington and Curtis (1986) equation or the Chapman–Richards equation. Though some of these differences can be attributed to the use of tree dominant height versus stand top height, significant differences in developmental patterns occur between naturally regenerated stands and plantations; possibly explaining some of the observed trends. Red alder is an aggressive pioneer species that regenerates strongly on sites with full sunlight and exposed mineral soil (Harrington et al., 1994). Naturally regenerated stands are generally very dense and often contain aggressive competitors like vine maple [Acer circinatum], salmonberry [Rubus spectabilis], and salal [Gaultheria shallon]. As suggested by the Harrington and Curtis (1986) equation, this may result in lower and earlier peak PAI values. In contrast, plantations are generally established with varying levels of vegetation control and with a more uniform spacing. These management activities decrease the level of both intra- and inter-competition at younger ages and may result in the observed a higher peak PAI values. Tree-level equations generally show a higher asymptote than equations based on stand-level data (Monserud, 1984; Raulier et al., 2003), which was supported by comparisons between the equations developed in this present study and the tree-level dominant height growth equation of Harrington and Curtis (1986). Examination of how well the Chapman–Richards and the Schumacher equations characterized the older, independent, natural stand data showed that both equations under-predicted HDM in the natural stand data, with the under-prediction being larger for the Chapman–Richards equation. We have concluded that the Schumacher equation likely provides more realistic predictions in older ages. However, the equations developed here are from a dataset with the maximum age value of 18 years (Table 3). Only the continued remeasurement and reanalysis of these permanent plots will definitively determine: (1) if red alder plantations reach a lower asymptotic top height than natural stands and (2) whether the Schumacher or the Chapman–Richards equation best characterizes the asymptotic behavior of plantation grown red alder. The results of this study indicate that planting density had no statistically significant influence on the parameters of the Schumacher GADA formulation, which determine the rate and shape of the top height growth equation. Because the GADA procedure used in this study eliminated the scaling parameter associated with SI that determines the asymptotic behavior of the two equations, there is no direct method to statistically test the potential impact of planting density on the scaling parameter. However, Fig. 4 indicates that planting densities below approximately 500 trees ha1 do affect the asymptote of the model. For the traditional Schumacher model form, the asymptote is defined by the plot’s SI, which is a multiplier on the exponential portion of the equation. In the anamorphic GADA formulation of the Schumacher model (Eq. (2)), the scaling parameter associated with SI is replaced by a function of AM and HDM. Thus, in this study, it is apparent that the impact of low planting density on HDP (and therefore SI) is adequately expressed through AM and HDM. The implication of this is that the impact of low planting density on HDP is the same as would occur by lowering the SI estimate for the plot of interest. Another implication of this finding is that predicted SI for stands below a planting density of 500 trees ha1 would require some type of correction if SI is to be used as a measure of site productivity in red alder plantations. Acknowledgements Thanks to the Hardwood Silviculture Cooperative and its supporting members for maintaining the network of plots and

A.R. Weiskittel et al. / Forest Ecology and Management 258 (2009) 323–331

providing access to the data. Special thanks to the Weyerhaeuser Company for also providing access to their data. Comments from Greg Johnson, David Marshall, Hubert Hasenauer, and two anonymous reviewers helped improved a previous version of this manuscript. References Amaro, A., Reed, D., Tome´, M., Themido, I., 1998. Modelling dominant height growth: Eucalyptus plantations in Portugal. Forest Science 44, 37–46. Bailey, R.L., Clutter, J.L., 1974. Base-age invariant polymorphic site curves. Forest Science 20, 155–159. Bluhm, A.A., Hibbs, D.E., 2006. Red alder: its management and potential. In: Deal, R.L., Harrington, C.A. (Eds.), Red Alder: A State of Knowledge. USDA Forest Service, Pacific Northwest Research Station, Portland, OR, pp. 73–86. Bormann, B.T., Gordon, J.C., 1984. Stand density effects in young red alder plantations: productivity, photosynthate partitioning, and nitrogen fixation. Ecology 65, 394–402. Briggs, D., 2007. Management practices on Pacific Northwest west-side industrial lands, 1991–2005: with projections to 2010. SMC Working Paper Number 6. Stand Management Cooperative, University of Washington, Seattle, WA, 79 pp. Calegario, N., Daniels, R.F., Maestri, R., Neiva, R., 2005. Modeling dominant height growth based on nonlinear mixed-effects model: a clonal Eucalyptus plantation case study. Forest Ecology and Management 204, 11–21. Carvalho, J., Parresol, B.R., 2005. A site model for Pyrenean oak (Quercus pyrenaica) stands using a dynamic algebraic difference equation. Canadian Journal of Forest Research 35, 93–99. Cieszewski, C.J., 2001. Three methods of deriving advanced dynamic site equations demonstrated on inland Douglas-fir site curves. Canadian Journal of Forest Research 31, 165–173. Cieszewski, C.J., 2003. Developing a well-behaved dynamic site index equation using a modified Hossfeld IV function Y3 = (axm)/(c + xm1), a simplified mixedmodel and scant subalpine fir data. Forest Science 49, 539–554. Cieszewski, C.J., Bailey, R.L., 2000. Generalized algebraic difference approach: theory based derivation of dynamic site equations with polymorphism and variable asymptotes. Forest Science 46, 116–126. Cieszewski, C.J., Bella, I.E., 1993. Modeling density-related lodgepole pine height growth, using Czarnowski’s stand dynamics theory. Canadian Journal of Forest Research 23, 2499–2506. Cieszewski, C.J., Strub, M.R., 2007. Parameter estimation of base-age invariant site index models: which data structure to use? A discussion. Forest Science 53, 552–555. Cieszewski, C.J., Strub, M.R., 2008. Generalized algebraic difference approach derivation of dynamic site equations with polymorphism and variable asymptotes from exponential and logarithmic functions. Forest Science 54, 303–315. Cieszewski, C.J., Strub, M.R., Zasada, M., 2007. New dynamic site equation that fits best the Schwappach data for Scots pine (Pinus sylvestris L.) in Central Europe Forest. Ecology and Management 243, 83–93. Cieszewski, C.J., Zasada, M., Strub, M.R., 2006. Analysis of different base models and methods of site model derivation for Scots pine. Forest Science 52, 187–197. Clutter, J.L., Fortson, J.C., Pienarr, L.V., Brister, G.H., Bailey, R.E., 1983. Timber Management: A Quantitative Approach. Wiley, New York, NY, 333 pp. Cole, E.C., Newton, M., 1987. Fifth-year responses of Douglas-fir to crowding and nonconiferous competition. Canadian Journal of Forest Research 17, 181–186. DeBell, D.S., Giordano, P.A., 1994. Growth patterns of red alder. In: Hibbs, D.E., DeBell, D.S., Tarrant, R.F. (Eds.), The Biology and Management of Red Alder. Oregon State University Press, Corvallis, OR, pp. 116–130. Die´guez-Aranda, U., Burkhart, H.E., Rodrı´guez-Soalleiro, R., 2005. Modeling dominant height growth of radiata pine (Pinus radiata D. Don) plantations in northwestern Spain. Forest Ecology and Management 215, 271–284. Flewelling, J.W., Collier, R., Gonyea, B., Marshall, D., Turnblom, E., 2001. Height–age curves for planted stands of Douglas-fir, with adjustments for density. Stand Management Cooperative Working Paper 1. University of Washington, College of Forest Resources, Seattle, WA, 25 pp. Gaiser, R.N., Merz, R.W., 1951. Stand density as a factor in estimating white oak site index. Journal of Forestry 49, 572–574. Giordano, P.A., Hibbs, D.E., 1993. Morphological response to competition in red alder: the role of water. Functional Ecology 7, 462–468.

331

Hann, D.W., Hanus, M.L., 2004. Evaluation of nonspatial approaches and equation forms used to predict tree crown recession. Canadian Journal of Forest Research 34, 1993–2003. Harrington, C.A., Curtis, R.O., 1986. Height growth and site index curves for red alder. Research Paper PNW-358. USDA Forest Service Pacific Northwest Research Station, Portland, OR, 14 pp. Harrington, C.A., Zasada, J.C., Allen, E.A., 1994. Biology of red alder (Alnus rubra Bong.). In: Hibbs, D.E., DeBell, D.S., Tarrant, R.F. (Eds.), The Biology and Management of Red Alder. Oregon State University Press, Corvallis, OR, pp. 3–22. King, J.E., 1966. Site index curves for Douglas-fir in the Pacific Northwest. Weyerhaeuser Forestry Paper 8. Western Forestry Research Center, Weyerhaeuser Company, Centralia, WA, 49 pp. Krumland, B., Eng, H., 2005. Site index systems for major young-growth forest and woodland species in northern California. California Forestry Report 4. Department of Forestry and Fire Protection, State of California Resources Agency, Sacramento, CA, 219 pp. Lynch, D.W., 1958. Effects of stocking on site measurement and yield of secondgrowth ponderosa pine in the Inland Empire. Research Paper 56. USDA, Forest Service, Intermountain Forest and Range Experiment Station, Ogden, UT, 36 pp. MacFarlane, D.W., Green, E.J., Burkhart, H.E., 2000. Population density influences assessment and application of site index. Canadian Journal of Forest Research 30, 1472–1475. McDill, M.E., Amateis, R.L., 1992. Measuring forest site quality using the parameters of a dimensionally compatible height growth function. Forest Science 38, 409– 429. Monserud, R.A., 1984. Height growth and site index curves for inland Douglas-fir based on stem analysis data and forest habitat type. Forest Science 30, 943–965. Monserud, R.A., 1985. Comparison of Douglas-fir site index and height growth curves in the Pacific Northwest. Canadian Journal of Forestry 15, 673–679. Nigh, G.D., Courtin, P.J., 1998. Height models for red alder (Alnus rubra Bong.) in British Columbia. New Forests 16, 59–70. Nord-Larsen, T., 2006. Developing dynamic site index curves for European beech (Fagus sylvatica L.) in Denmark. Forest Science 52, 173–181. Oliver, C.W., Larson, B.C., 1996. Forest Stand Dynamics, Update edition. Wiley, New York, NY, 544 pp. Porter, D.R., Wiant Jr., H.V., 1965. Site index equations for tanoak, pacific madrone, and red alder in the Redwood region of Humboldt Co., California. Journal of Forestry 63, 286–287. Raulier, F., Lambert, M.C., Pothier, D., Ung, C.H., 2003. Impact of dominant tree dynamics on site index curves. Forest Ecology and Management 184, 65–78. ˜ ellas, I., Montero, G., 2008. Base-age invariant cork Sa´nchez-Gonza´lez, M., Can growth model for Spanish cork oak (Quercus suber L.) forests. European Journal of Forest Research 127, 173–182. Scott, W., Meade, R., Leon, R., Hyink, D., Miller, R., 1998. Planting density and tree-size relations in coast Douglas-fir. Canadian Journal of Forest Research 28, 74–78. Sharma, M., Amateis, R.L., Burkhart, H.E., 2002. Top height definition and its effect on site index determination in thinned and unthinned loblolly pine plantations. Forest Ecology and Management 168, 163–175. Skovsgaard, J.P., Vanclay, J.K., 2008. Forest site productivity: a review of the evolution of dendrometric concepts for even-aged stands. Forestry 81, 13–31. Smith, V.G., 1984. Asymptotic site-index curves, fact or artifact? Forestry Chronicle 60, 150–156. Wang, M., Borders, B.E., Zhao, D., 2007a. Parameter estimation of base-age invariant site index models: which data structure to use? Forest Science 53, 541–551. Wang, M., Borders, B.E., Zhao, D., 2008a. An empirical comparison of two subjectspecific approaches to dominant heights modeling: the dummy variable method and the mixed model method. Forest Ecology and Management 255, 2659–2669. Wang, M., Borders, B.E., Zhao, D., 2008b. Parameter estimation of base-age invariant site index models: which data structure to use? A reply. Forest Science 54, 129–133. Wang, Y., LeMay, V.M., Baker, T.G., 2007b. Modelling and prediction of dominant height and site index of Eucalyptus globulus plantations using a nonlinear mixed-effects model approach. Canadian Journal of Forest Research 37, 1390–1403. Worthington, N.P., Johnson, F.A., Staebler, G.R., Lloyd, W.J., 1960. Normal yield tables for red alder. Research Paper 36. USDA Forest Service, Pacific Northwest Forest and Range Experiment Station, Portland, OR, 32 pp. Xu, R.H., 2003. Measuring explained variation in linear mixed effect models. Statistics in Medicine 22, 3527–3541.