Abstract
A new stand-level dynamics model based on observed stand growth trajectories is presented. This stand-level dynamics model uses the trajectories of observed plots through Reineke's quadratic mean diameter-density space to predict the change in quadratic mean diameter, ingrowth and mortality over time. The model uses the collection of observed trajectories as a probability distribution that is used to guide an informed random walk. An imputation model is used to select k nearest neighbors (bandwidth) which are then used to build joint kernel distributions. From these kernel distributions, m random samples (sample intensity) are averaged to predict the change in quadratic mean diameter, ingrowth and mortality. All levels of k tested (10, 20, 30) performed well as long as sampling intensity was above 1. Variability in predictions was reduced at sampling intensities above 1, but no significant differences were visible among sampling intensities 5 and above.
Introduction
Stand-level growth and yield models predict stand average summary attributes such as density (trees per unit area), basal area (sum of tree cross-sectional area per unit area), volume per unit area, and some mean tree attributes such as quadratic mean tree diameter and mean height (Weiskittel et al., 2011). Stand growth and yield models, in the form of yield tables, have been in use in North America for over 100 years (Husch et al., 2003; Weiskittel et al., 2011) and for over 200 years in Europe (Pretzsch, 2008). In the 1930s and 1940s, mathematical formulations of stand-level growth and yield models were developed as an alternative to the traditional tabular or graphical presentation (Schumacher, 1939). In the early 1960s, Buckman (1962) and Clutter (1963), using differential calculus, developed compatible growth and yield equations. Furnival and Wilson (1971) further developed the concept of compatible equations, expanding the components modeled and used systems of equations to solve for all parameters simultaneously. Moser and Hall (1969) applied the systems of equations approach to uneven-aged stands.
The past few decades have seen a shift toward the development of individual tree growth and yield models (Weiskittel et al., 2011); however, stand-level models still play an important role in many forestry and ecological applications. Long-term strategic forest management plans are often modeled at the stand level and thus require stand-level growth and yield projections. While these stand-level projections are often summaries of predictions from individual tree growth and yield models, improved model predictions can be achieved by using a stand-level model in conjunction with an individual tree model (Weiskittel et al., 2011).
One formulation of a stand-level model that has not received much application, but has interesting potential for use in conjunction with an individual tree model is the direction fields proposed by Leary and Stephans (1987). They modeled basal area per unit area over age as a series of short line segments representing the average direction (i.e. the slope of the basal area − age growth curve) stands moved through the basal area − age space. Direction fields could be developed for any pair of stand-level parameters (or even higher dimensions).
When the stand-level parameter pair is the quadratic mean diameter − density space, the direction fields would define the self-thinning space as proposed by Reineke (1933). While most research on self-thinning and self-thinning trajectories has focused on estimating the maximum boundary (e.g. Yang and Titus, 2002; VanderSchaaf, 2010), a few researchers have tried to derive growth and yield models from the self-thinning concept (e.g. Smith and Hann, 1984, 1986, Cao, 1994, Arisman et al., 2005). If the direction moved in the self-thinning space was defined by changes in quadratic mean diameter, ingrowth and mortality, then the direction fields approach would describe a system similar to the systems of equation approach to modeling stand-level growth and yield. Cao's (1994) approach was similar to this idea, modeling the changes in quadratic mean diameter and tree survival in even-aged loblolly pine plantations. The inclusion of ingrowth would generalize the system to be applicable to a broader range of stand conditions.
While Leary and Stephans (1987) represented their direction fields using the mean direction, the direction fields could be viewed as a probability distribution of directions and the direction a stand moves modeled as a random walk based on the direction field distributions. We propose the development of a random walk, stand-level growth and yield model based on empirical multidimensional copulas derived from nearest neighbor imputation. A relatively simple imputation model is developed to predict stand trajectories through the self-thinning space for mixed species stands composed of single- to multi-cohorts. The model predicts the changes in quadratic mean diameter, ingrowth and mortality. The influences of bandwidth (number of nearest neighbors used to develop the kernel density distributions) and sample size (number of random samples from the multivariate copula used to estimate the direction) are explored in the paper.
Methods
Data
Data for this analysis came from a network of permanent sample plots across New Brunswick (NB) and Nova Scotia (NS), Canada, and forest inventory and analysis (FIA) plots in Maine, USA (ME). Plots spanned NB, NS and ME, and consisted of a mixture of single- to multi-species, single- to multi-cohort stands. In NB and NS, plots were a single 400 m2 circular plot, while in ME, the plots were the FIA Phase 2/Phase 3 plot cluster design (Bechtold and Patterson, 2005). In NB, plot establishment began in 1987 and continues today, with plots typically remeasured every 3–5 years. Plot measurements in NS began as early as 1965, with an initial establishment period continuing until 1970. A renewed plot establishment effort was started in 1990. Plots are remeasured on a 5-year cycle in NS. Because of the changes in sampling grid and intensity, and problems linking older panels to the new panels, only plot measurements since 1998 were used for the ME FIA data.
Minimum threshold diameters for measurement varied across the regions. Because of issues with data quality, and for consistency across all regions, all measurements of trees with diameter at breast height <9.0 cm were dropped for this study. Table 1 shows the distribution of species by abundance and number of plots present.
Percent occurrence of common tree species by data source
| Species | Percent occurrence by number of plots | Percent occurrence by basal area per ha | ||||
|---|---|---|---|---|---|---|
| NS | NB | ME | NS | NB | ME | |
| Balsam fir (Abies balsamea) | 78.7 | 76.9 | 70.7 | 22.4 | 19.7 | 12.9 |
| Red spruce (Picea rubens) | 69.9 | 45.0 | 55.9 | 19.2 | 18.1 | 11.2 |
| Black spruce (Picea nigra) | 23.2 | 28.8 | 9.0 | 5.0 | 14.1 | 2.0 |
| White spruce (Picea glauca) | 34.0 | 33.8 | 21.3 | 7.6 | 6.2 | 2.3 |
| White pine (Pinus strobus) | 19.1 | 8.4 | 29.3 | 5.3 | 1.0 | 8.1 |
| Red maple (Acer rubrum) | 66.4 | 59.0 | 67.2 | 16.1 | 9.1 | 12.9 |
| Yellow birch (Betula alleghaniensis) | 26.6 | 27.4 | 41.1 | 5.3 | 3.2 | 5.5 |
| Sugar maple (Acer saccharum) | 13.5 | 21.1 | 25.7 | 4.4 | 5.7 | 5.8 |
| Paper birch (Betula papyrifera) | 36.6 | 46.6 | 49.0 | 3.6 | 5.2 | 5.5 |
| Quaking aspen (Populus tremuloides) | 10.4 | 23.8 | 20.6 | 1.5 | 6.8 | 2.6 |
| Species | Percent occurrence by number of plots | Percent occurrence by basal area per ha | ||||
|---|---|---|---|---|---|---|
| NS | NB | ME | NS | NB | ME | |
| Balsam fir (Abies balsamea) | 78.7 | 76.9 | 70.7 | 22.4 | 19.7 | 12.9 |
| Red spruce (Picea rubens) | 69.9 | 45.0 | 55.9 | 19.2 | 18.1 | 11.2 |
| Black spruce (Picea nigra) | 23.2 | 28.8 | 9.0 | 5.0 | 14.1 | 2.0 |
| White spruce (Picea glauca) | 34.0 | 33.8 | 21.3 | 7.6 | 6.2 | 2.3 |
| White pine (Pinus strobus) | 19.1 | 8.4 | 29.3 | 5.3 | 1.0 | 8.1 |
| Red maple (Acer rubrum) | 66.4 | 59.0 | 67.2 | 16.1 | 9.1 | 12.9 |
| Yellow birch (Betula alleghaniensis) | 26.6 | 27.4 | 41.1 | 5.3 | 3.2 | 5.5 |
| Sugar maple (Acer saccharum) | 13.5 | 21.1 | 25.7 | 4.4 | 5.7 | 5.8 |
| Paper birch (Betula papyrifera) | 36.6 | 46.6 | 49.0 | 3.6 | 5.2 | 5.5 |
| Quaking aspen (Populus tremuloides) | 10.4 | 23.8 | 20.6 | 1.5 | 6.8 | 2.6 |
Percent occurrence of common tree species by data source
| Species | Percent occurrence by number of plots | Percent occurrence by basal area per ha | ||||
|---|---|---|---|---|---|---|
| NS | NB | ME | NS | NB | ME | |
| Balsam fir (Abies balsamea) | 78.7 | 76.9 | 70.7 | 22.4 | 19.7 | 12.9 |
| Red spruce (Picea rubens) | 69.9 | 45.0 | 55.9 | 19.2 | 18.1 | 11.2 |
| Black spruce (Picea nigra) | 23.2 | 28.8 | 9.0 | 5.0 | 14.1 | 2.0 |
| White spruce (Picea glauca) | 34.0 | 33.8 | 21.3 | 7.6 | 6.2 | 2.3 |
| White pine (Pinus strobus) | 19.1 | 8.4 | 29.3 | 5.3 | 1.0 | 8.1 |
| Red maple (Acer rubrum) | 66.4 | 59.0 | 67.2 | 16.1 | 9.1 | 12.9 |
| Yellow birch (Betula alleghaniensis) | 26.6 | 27.4 | 41.1 | 5.3 | 3.2 | 5.5 |
| Sugar maple (Acer saccharum) | 13.5 | 21.1 | 25.7 | 4.4 | 5.7 | 5.8 |
| Paper birch (Betula papyrifera) | 36.6 | 46.6 | 49.0 | 3.6 | 5.2 | 5.5 |
| Quaking aspen (Populus tremuloides) | 10.4 | 23.8 | 20.6 | 1.5 | 6.8 | 2.6 |
| Species | Percent occurrence by number of plots | Percent occurrence by basal area per ha | ||||
|---|---|---|---|---|---|---|
| NS | NB | ME | NS | NB | ME | |
| Balsam fir (Abies balsamea) | 78.7 | 76.9 | 70.7 | 22.4 | 19.7 | 12.9 |
| Red spruce (Picea rubens) | 69.9 | 45.0 | 55.9 | 19.2 | 18.1 | 11.2 |
| Black spruce (Picea nigra) | 23.2 | 28.8 | 9.0 | 5.0 | 14.1 | 2.0 |
| White spruce (Picea glauca) | 34.0 | 33.8 | 21.3 | 7.6 | 6.2 | 2.3 |
| White pine (Pinus strobus) | 19.1 | 8.4 | 29.3 | 5.3 | 1.0 | 8.1 |
| Red maple (Acer rubrum) | 66.4 | 59.0 | 67.2 | 16.1 | 9.1 | 12.9 |
| Yellow birch (Betula alleghaniensis) | 26.6 | 27.4 | 41.1 | 5.3 | 3.2 | 5.5 |
| Sugar maple (Acer saccharum) | 13.5 | 21.1 | 25.7 | 4.4 | 5.7 | 5.8 |
| Paper birch (Betula papyrifera) | 36.6 | 46.6 | 49.0 | 3.6 | 5.2 | 5.5 |
| Quaking aspen (Populus tremuloides) | 10.4 | 23.8 | 20.6 | 1.5 | 6.8 | 2.6 |
Individual plot summaries
Plot summaries were calculated for each measured plot in NS and NB. In ME, the FIA clusters were averaged to obtain a single-cluster summary. For each remeasurement, the following periodic stand characteristics were calculated: quadratic mean diameter in centimeters (Dq), density in trees per hectare (N), ingrowth density in trees per hectare (I), mortality density in trees per hectare (M), and species-weighted additive relative density (SWARD). SWARD is a stand density measure developed by Ducey and Knapp (2010) that weights relative density of a multi-species stand using specific gravity. Because the focus of this study is on endogenous stand dynamics, rather than propensity of stands to exogenous natural disturbance or forest harvesting, plots with catastrophic mortality, defined as any plot having >30 percent loss of density over any remeasurement period and any plot experiencing human intervention were removed from the data set. The number of plots removed from NS, NB and ME data sets were 285, 110 and 625 respectively. Total number of plots (plot clusters for FIA data) used for NS, NB and ME were 2615, 1439 and 2475, respectively. Measurements in these plots span 39 years in NS, 23 years in NB and 9 years in ME.
Trajectories were calculated for the change in quadratic mean diameter and change in density for each plot. The change in quadratic mean diameter was the difference between quadratic mean diameter between measurement i and measurement i + 1. The change in density was expressed in terms of ingrowth and mortality. Periodic changes were calculated for each measurement period and annualized by linear interpolation (McDill and Amateis, 1993).
The overall behavior of individual trajectories through the self-thinning space was summarized using direction fields (Leary and Stephans, 1987). Trajectories were assigned to bins across the range of density and quadratic mean diameter, and average trajectory slope was calculated for each bin (Figure 1). The bins were equally spaced across the observed range of the self-thinning space (in log–log space) that created a 50 × 50 grid. Each grid space was a bin and all slopes within each bin were averaged. The average slopes provide a broad understanding of plot behavior across the self-thinning space and will be used as a benchmark for the random walk model developed subsequently.
Average observed trajectories—direction fields for observed data.
Random walk model
The stand model was developed in R (R Development Core Team, 2011). The structure of the stand model is shown in Figure 2. The basic structure of the informed random walk model is:
A starting point, defined by Dq, N and SWARD, is specified.
Nearest neighbor imputation is then used to develop a joint distribution of Ii, Mi and, ΔDqi Only three independent variables were used in the imputation process in this study: Dqi, Ni and SWARDi. The R package yaImpute (Crookston and Finley, 2008) was used for nearest neighbor imputation. This package requires a set of reference points and a set of target points. The reference points are a matrix of known stand characteristics [Dq, N, SWARD] and the target points are a matrix of the variables of interest [I, M, ΔDq] associated with the reference set. A weighted Euclidean distance function (Crookston and Finley, 2008) is used to rank the reference points from most similar to least similar, relative to the current trajectory state (Dqi, Ni, SWARDi). The k nearest neighbors are then selected from the reference points and the associated subset of attributes selected from the matrix of known stand attributes. The number of nearest neighbors selected (k) is set by the user, and is referred to as bandwidth in this study.
The subset of stand attributes [Ik, Mk, ΔDqk] selected in step 2 is then used to estimate the marginal kernel densities for Ii, Mi, ΔDqi. Kernel density estimates are obtained using the density function in R (Becker et al., 1988; Scott, 1992). Density estimates were calculated from the k nearest neighbors selected in the imputation step. An optional weighting function, based on the inverse of the Euclidean distance calculated in the imputation step, may be used in the density estimation.
The value for each step in the random walk is then determined by sampling m observations from the joint distribution of the marginal kernel density estimates using a normal copula (Genest and MacKay, 1986; Yan, 2007). The copula sampling approach used here was the same as described by Kershaw et al. (2010): (a) an m-by-3 matrix of random standard normal variates (µ = 0, σ = 1) were generated; (b) the columns of the random normal variate matrix were correlated using the Choleski decomposition (Andersen et al., 1999) of the partial correlation matrix estimated using Spearman's rank correlation and the k values of Ii, Mi, ΔDqi selected in step 2; (c) the normal marginals were then stripped off using the inverse (cumulative) normal distribution; (d) the m sample values for Ii, Mi, ΔDqi are then obtained using a custom quantile function (Kershaw, 2011), the cumulative probabilities and the marginal kernel density estimates. The values for Ii, Mi, ΔDqi defining the step are obtained by averaging the m random samples.
- The values for the step are then used to update Dq and N. SWARD is updated using a proportional change in SDI:where SDIi is stand density index (McCarter and Long, 1986) at time i and SWARDi is species-weighted average relative density at time i. Using a proportional change in SDI assumes that the species compositions used to calculate SWARD do not change.(4)
Steps 2–5 are repeated for the desired length of the informed random walk (years of projection).
Informed random walk model flow.
Determination of bandwidth and sample intensity
A range of grid points spanning the observed self-thinning space were used to assess potential impacts of bandwidth and sample size on model predictions across the self-thinning space (Table 2). For this stand-level imputation model, bandwidth refers to the number of neighbors sampled (k) to create the kernel distributions in step 3 described above. Sample intensity is the number of random samples drawn from the normal copula (m) described in step 4. Bandwidth and sample intensity were tested individually. To test bandwidth, the number of random samples (m) was held constant at 10 while 3 values of k were used: 10, 20 and 30. To test the effect of the number of neighbors sampled (k), the weighting described in step 3 above was not used as there was concern that any potential effect of k would be masked by the weighting process. Four levels of sample intensity (m) were tested: 1, 5, 10 and 20 while k remained at 20. These four levels were tested twice, once with weighting and once with no weighting. At each grid point × bandwidth × sample intensity level, 1000 random walks were initiated and run for 20 consecutive steps. In total, 11 model runs were initiated at each grid point – 3 levels of k (unweighted), 4 levels of m (weighted) and 4 levels of m (unweighted). An additional run to further test the influence of k was done on a subset of grid points. This additional run was a factorial analysis that used all k levels (10, 20, 30) and ran each k level for each sampling intensity (1, 5, 10, 20). Observed plots with initial density within ±50 trees per hectare and initial quadratic mean diameter within ±0.5 cm were compared with the random walk results. A maximum of 10 observed plots were examined at each grid point. Predictions of Ni, Dqi, Ii and Mi were compared with observed values in 5-year intervals beginning with year 5.
Sample space (grid cells) used for testing bandwidth and sample size for random walk model validation
| Quadratic mean diameter (cm) | Density (trees ha−1) | ||||||
|---|---|---|---|---|---|---|---|
| 250 | 500 | 1000 | 1500 | 2000 | 2500 | 3000 | |
| 32 | • | ||||||
| 30 | • | • | |||||
| 28 | • | • | • | ||||
| 26 | • | • | • | ||||
| 24 | • | • | • | ||||
| 22 | • | • | • | ||||
| 20 | • | • | • | • | |||
| 18 | • | • | • | • | • | • | |
| 16 | • | • | • | • | • | • | |
| 14 | • | • | • | • | • | • | • |
| 12 | • | • | • | • | • | • | • |
| Quadratic mean diameter (cm) | Density (trees ha−1) | ||||||
|---|---|---|---|---|---|---|---|
| 250 | 500 | 1000 | 1500 | 2000 | 2500 | 3000 | |
| 32 | • | ||||||
| 30 | • | • | |||||
| 28 | • | • | • | ||||
| 26 | • | • | • | ||||
| 24 | • | • | • | ||||
| 22 | • | • | • | ||||
| 20 | • | • | • | • | |||
| 18 | • | • | • | • | • | • | |
| 16 | • | • | • | • | • | • | |
| 14 | • | • | • | • | • | • | • |
| 12 | • | • | • | • | • | • | • |
Sample space (grid cells) used for testing bandwidth and sample size for random walk model validation
| Quadratic mean diameter (cm) | Density (trees ha−1) | ||||||
|---|---|---|---|---|---|---|---|
| 250 | 500 | 1000 | 1500 | 2000 | 2500 | 3000 | |
| 32 | • | ||||||
| 30 | • | • | |||||
| 28 | • | • | • | ||||
| 26 | • | • | • | ||||
| 24 | • | • | • | ||||
| 22 | • | • | • | ||||
| 20 | • | • | • | • | |||
| 18 | • | • | • | • | • | • | |
| 16 | • | • | • | • | • | • | |
| 14 | • | • | • | • | • | • | • |
| 12 | • | • | • | • | • | • | • |
| Quadratic mean diameter (cm) | Density (trees ha−1) | ||||||
|---|---|---|---|---|---|---|---|
| 250 | 500 | 1000 | 1500 | 2000 | 2500 | 3000 | |
| 32 | • | ||||||
| 30 | • | • | |||||
| 28 | • | • | • | ||||
| 26 | • | • | • | ||||
| 24 | • | • | • | ||||
| 22 | • | • | • | ||||
| 20 | • | • | • | • | |||
| 18 | • | • | • | • | • | • | |
| 16 | • | • | • | • | • | • | |
| 14 | • | • | • | • | • | • | • |
| 12 | • | • | • | • | • | • | • |
As time increased the number of observed plots over which RMSE was calculated decreased because some plots did not have a 20-year measurement period.
Results
General model behavior
Mapped error contours
RMSE for mortality (Figure 3) showed very little increase in error as time interval (i) increased. At 5 years, RMSE for mortality was <60trees ha−1. Mortality RMSE was higher when associated with lower quadratic mean diameters and increased towards the maximum size-density line (i.e. increased with increasing density when Dq is held constant, Figure 3). Time did not affect the prediction of ingrowth (Figure 4). The highest ingrowth RMSE was seen in the lower density–lower quadratic mean diameter region of the self-thinning space (Figure 4), and decreased as the maximum size-density line was approached and as Dq increased. The mapped error contour for Dq is not presented as no patterns were evident over the self-thinning space. RMSEs for Dq were generally small, and remained between 1 and 2 cm across the majority of the self-thinning space.
Contour map for mortality RMSE. Based on sample intensity (m) = 10 and bandwidth (k) = 20.
Contour map for ingrowth RMSE. Based on sample intensity (m) = 10 and bandwidth (k) = 20.
Trajectory behaviors
Examples of the 1000 random walk trajectories, each with 20 steps, are shown in Figure 5. Rather than showing each individual trajectory, the 95th percentile loess bounds for each 5-year period are shown (Galili, 2010). Note that the width of each of these bounds changed depending on where the starting point was in the self-thinning space (Figure 5). Near the maximum size-density line, the bounds of the trajectories were narrower than other regions in the self-thinning space. For each starting point, trajectory bound widths increased in size as time advanced (Figure 5). Generally errors associated with density were greater than errors associated with Dq.
Model trajectory 95 percentile bounds for 5 starting points on the self-thinning space. Each starting point used a sampling intensity (m) = 10 and a bandwidth (k) = 20. Trajectory bounds shown are five steps (years) per each band.
Bandwidth selection
Error patterns for all variables (Dq, I, and M) were similar across all bandwidth selection sizes (k). All results indicate that there was very little difference between the three k levels used here (10, 20, and 30 nearest neighbors). Figures 6 and 7 show the RMSE for ingrowth by both density (Figure 6) and quadratic mean diameter (Figure 7). In general, errors peaked at higher density levels. Errors were consistently larger for smaller quadratic mean diameters and decreased rapidly with increasing quadratic mean diameter. Bandwidth only seemed to make a consistent difference, though not significant (P > 0.05), at the higher densities (Figures 6 and 7): k = 10 had lower average errors as density increased. Even with the use of a sampling intensity of 1, no difference in RMSE is evident (Figure 8).
Root-mean-square error for ingrowth versus density for each time interval and each k value (sampling intensity = 10).
Root-mean-square error for ingrowth versus quadratic mean diameter for each time interval and each k value (sampling intensity = 10).
Root-mean-square error for ingrowth versus density for each time interval for each k value using a sampling intensity (m) = 1.
Sample intensity selection
Sample intensity results were very similar across variables (Dq, N, I and M), therefore, only mortality results are presented (Figures 9 and 10). Sampling only one random sample from the kernel density results in the highest RMSE (Figures 9 and 10). There was little difference between sampling intensities 5, 10 and 20 (Figures 9 and 10). Weighting did not have much influence on RMSE results (Figure 9 versus Figure 10). Mortality RMSE peaked at mid-range densities and decreased or plateaued above densities of 2000trees ha−1 (Figures 9 and 10). In general, errors did not increase over time (Figures 9 and 10).
Root-mean-square error for unweighted runs: mortality versus density for each time interval and each sampling intensity (m) value.
Root-mean-square error for weighted runs: mortality versus density for each time interval and each sampling intensity (m) value.
Discussion
General model behavior
Mapped error contours
An increase in RMSE over time was evident in the mortality contour error map (Figure 3), but the location of the highest errors remained the same over time. The pattern of errors across the self-thinning space in any time period is of interest here. Mortality RMSE increased as proximity to the maximum self-thinning line increased. The increase in mortality RMSE in proximity to the maximum self-thinning line is a reflection of the increase in mortality occurring in this region. Because our model averages over a number of samples drawn from the copula, the mortality rates are smoothed over time. Mortality in stands subject to high self-thinning pressure often occurs in pulses rather than continuously (Arisman et al., 2005). Other regions in the self-thinning space experienced less mortality, thus smaller prediction errors occur. In the region where mortality was occurring at increased levels, predictions could be improved by incorporating additional variables in the imputation. This model encompasses multi-cohort, multi-species stands and currently only SWARD is included in the imputation model to capture the multi-species effects. Species composition should affect mortality as each species competes for resources differently. For example, if shade-tolerant species are present, self-thinning may be slowed and mortality reduced. Site factors may also influence mortality; Weiskittel et al. (2009) found that the location of the self-thinning line was influenced by both site index and purity (measured as the proportion of basal area in the primary species). While we were not estimating the self-thinning line, it is reasonable to conclude that if the self-thinning line is sensitive to site and species composition, the amount of mortality occurring can also be sensitive to these same factors and RMSE for mortality may have decreased if these variables had been included in the informed random walk model. While SWARD may at least partially compensate for species composition in the location of the self-thinning line (Ducey and Knapp, 2010), species composition may still impact the kinetics of the approach to the line. Full exploration of this issue is outside the scope of this study but is worth further attention.
Ingrowth RMSE increased in the regions of the self-thinning space where one expects ingrowth to be occurring – in the lower density, smaller quadratic mean diameter regions. RMSE was lower in other regions of the self-thinning space simply because ingrowth was not occurring as much in these areas. Increasing the number of variables used in the model may reduce RMSE for ingrowth. Variables including species composition, tree size and site conditions were noted by Shifley et al. (1993) as potential factors influencing ingrowth. Many authors have found the basal area to be an important variable predicting ingrowth (e.g. Moser 1972; Ek, 1974). In their stand level ingrowth model, Li et al. (2011) note the importance of basal area and percent hardwood. Site may also be used as a predictor of ingrowth, although some authors have found it to be non-significant (e.g. Ek, 1974), others have found ingrowth to be higher on better sites (Tecle and Covington, 1991). Including any of these variables (basal area, species composition, or site) has the potential to better capture observed variability and reduce prediction errors for ingrowth.
No patterns emerged in RMSE for Dq across the self-thinning space, average RMSE remained between 1 and 2 cm across the majority of the self-thinning space. Much like ingrowth and mortality, the prediction of growth in quadratic mean diameter has the potential to be improved by the inclusion of additional variables in the informed random walk model. The same factors that affect mortality and ingrowth (e.g. species composition and site index) may influence growth in quadratic mean diameter.
Trajectory behaviors
Variation in random walk trajectories changed over the self-thinning space (Figure 5). The width of each 5 year bound represents the uncertainty in the predictions. Starting points closest to the maximum self-thinning line had smaller trajectory bounds than those starting points in the center or in regions with higher quadratic mean diameters and lower densities (Figure 5). Variation in the bounds of the random walk trajectories was linked to variation in observed trajectories. Stands closest to the maximum size-density line had increased error in mortality, indicative of the increased mortality occurring in observed trajectories, but error in ingrowth approached zero. The reduction in ingrowth error in this area of the self-thinning space resulted in more predictable trajectories and though mortality errors in this region were higher than in other regions of the self-thinning space, the mortality occurring in this region was more predictable relative to the variables used in the self-thinning space and was bounded by the maximum size–density relationship. When mortality occurred at lower densities, it was less predictable than when it occurred near the maximum size-density line. At lower densities, this less predictable mortality, and the increase in ingrowth occurring (and thus ingrowth error) resulted in more variation in predicted trajectories.
This model has the potential to be used to examine where predictions are most variable. Rather than a model that can examine errors only by time or covariate, this model offers a method to examine prediction errors throughout stand development. Expanding Leary and Stephans’ (1987) direction fields to build a distribution of trajectories to create this informed random walk model reveals more about observed stand dynamics and development than models based on individual plot points.
The width of the trajectory bounds increased as time advanced (Figure 5), though the RMSE in ingrowth and mortality did not increase as time advanced (Figures 6–10). The increase in the width of the trajectory bounds as time increased was not surprising as neighbors are sampled at each time step and the resulting step in each walk could veer left or right depending on the sample taken from the joint distribution built on those neighbors. Each time step then could veer further to the left or right than previous time steps even though the bandwidth has not changed. The lack of substantial increase in RMSE for either ingrowth or mortality as time advanced indicates that, on average, each step of informed random walk was still predicting the observed value well. Put simply, the veering to the right was canceled out by the veering to the left, and vice versa.
Variation in the trajectory bounds was always wider in the density axis than in the quadratic mean diameter axis (Figure 5). The increased variation in the density axis was a reflection of the growth pattern in the observed trajectories. Observed changes in density tended to be higher than observed changes in quadratic mean diameter.
For most of the self-thinning space, the trajectory bands reflected the average direction fields (Figure 5). In stands with low Dq and lower N, the trajectory bands appeared to move to the right (increasing N) while many of the direction fields appeared to suggest a vertical movement (increasing Dq with little change in N). The direction fields showed an average slope and did not convey any information about variability. In addition, the use of SWARD in our model shifted these direction fields as different subsets of plots were averaged. The combination of high variability in these low-density stands and the differential weighting introduced by the inclusion of SWARD most likely accounted for the apparent difference in the observed direction field and simulated trajectory bounds.
Bandwidth selection
The number of neighbors selected (k) to build the kernel density that is randomly sampled from did not significantly influence the RMSE for any of the variables (I, M, Dq, or N). Figures 6 and 7 show the ingrowth RMSE for each 5-year time period and k level. Even though weighting was not employed in these model runs, k still did not influence RMSE. The additional factorial runs also show no significant difference between k values when sampling intensity is at 1 (Figure 8). The visible difference between k values at higher densities was a result of the neighbors available to sample at these densities. At higher densities only one side of the distribution was available for sampling because there were few neighbors above the current density, thus when bandwidth increased neighbors selected were farther away. Within the range of bandwidths tested, the number of neighbors selected to build the joint kernel distribution was not as important as the number of random samples drawn from the distribution and averaged.
Sample intensity selection
Increasing sampling intensity above 1 resulted in lower RMSE (Figures 9 and 10). Improvements in RMSE above sampling intensity of five were consistent for both weighted and unweighted runs (Figures 9 and 10) and this pattern held for all variables examined. Weighting did not increase the benefit of increased sampling intensity (Figures 9 and 10).
Model applications
The informed random walk model is a unique approach to predicting ingrowth, mortality and quadratic mean diameter. The model bases the imputation on the trajectories of individual observations, and, as such, has the ability to capture stand dynamics in a new way. Bounds of multiple trajectories as shown in Figure 5 could be used to predict the range of potential stand growth over time. Examining the bounds of the random walk trajectories gives an indication of the certainty of the predictions and how much variation is seen in the observed trajectories. Bounds from these trajectories can provide a sense of what is biologically possible for stands in the region, and can be used to test new or existing empirical models of any scope. The trajectory through the self-thinning space can be calculated based on the predictions of any new or existing empirical model. If those trajectories fall outside the bounds of the trajectories seen in the informed random walk, then existing or new empirical models may not be capturing observed stand development.
Informed random walk trajectories could be used in a growthandyield model in several ways. One option is to use a predicted trajectory in a stand-level growth and yield model as a stochastic feature. Incorporating the informed random walk as a stochastic feature would better reflect the variability in the development of multi-cohort, multi-species stands. As a stochastic model, the trajectory of an individual stand would vary each time the model is run, and this is not always desirable in forest management. As such, this model could also be used as an empirical model if trajectories were averaged, or if a particle filter (Gove, 2009) was used to weight and smooth predicted trajectories.
The informed random walk model can also be used in conjunction with a tree-level model. While this model predicts the per hectare values of mortality and ingrowth, additional models could be used to predict which individual trees die in a mortality model (Monserud, 1976; Monserud and Sterba, 1999), or what ingrowth species occurs (Li et al., 2011). Several methods of linking a stand-level and tree-level model exist including disaggregation, constrained or combined approaches (Weiskittel et al., 2011). Disaggregation models use a top-down approach where stand-level models predict per unit area changes, and then these predictions are then allocated to individual trees (e.g. Ritchie and Hann, 1997; Cao, 2006; Qin and Cao, 2006; Zhang et al., 2011). In contrast, constrained models use a bottom-up approach where tree-level predictions are optimized at several levels (e.g. Cao, 2006). The combined approach uses a system where both individual tree-level and stand-level predictions are made and then modified through a feedback system (e.g. Yue et al., 2008; Zhang et al., 2011). Predictions from the informed random walk model could be used as the stand-level model in both the disaggregation and combined methods of using stand-level models and tree-level models in conjunction.
Conclusions
The informed random walk model is a unique stand-level model capable of predicting three main components of stand dynamics: net growth (ΔDq), ingrowth (I) and mortality (M). This model can be used on its own to examine the stand dynamics, or in conjunction with a tree-level model to improve overall predictions. By simulating many trajectories from the same starting point, as shown in Figure 5, the model provides its own measure of certainty – narrow bounds indicate increased certainty. The informed random walk has the benefit of not needing to estimate the size-density line to produce the typical self-thinning behavior. This is a major benefit as there is no need to subsample plots experiencing self-thinning to estimate the maximum size-density line. Finally, this model is not restricted to even-aged or single species stands, as long as sufficient observed data are available, this model is applicable to all stand structures and species compositions.
