Tree Height Increment Models for National Forest Inventory Data in the Pacific Northwest, USA

The United States national inventory program measures a subset of tree heights in each plot in the Pacific Northwest. Unmeasured tree heights are predicted by adding the difference between modeled tree heights at two measurements to the height observed at the first measurement. This study compared different approaches for directly modeling 10-year height increment of red alder (RA) and ponderosa pine (PP) in Washington and Oregon using national inventory data from 2001–2015. In addition to the current approach, five models were implemented: nonlinear exponential, log-transformed linear, gamma, quasi-Poisson, and zero-inflated Poisson models using both tree-level (e.g., height, diameter at breast height, and compacted crown ratio) and plot-level (e.g., basal area, elevation, and slope) measurements as predictor variables. To account for negative height increment observations in the modeling process, a constant was added to shift all response values to greater than zero (log-transformed linear and gamma models), the negative increment was set to zero (quasi-Poisson and zero-inflated Poisson models), or a nonlinear model, which allows negative observations, was used. Random plot effects were included to account for the hierarchical data structure of the inventory data. Predictive model performance was examined through cross-validation. Among the implemented models, the gamma model performed best for both species, showing the smallest root mean square error (RSME) of 2.61 and 1.33 m for RA and PP, respectively (current method: RA—3.33 m, PP—1.40 m). Among the models that did not add the constant to the response, the quasi-Poisson model exhibited the smallest RMSE of 2.74 and 1.38 m for RA and PP, respectively. Our study showed that the prediction of tree height increment in Oregon and Washington can be improved by accounting for the negative and zero height increment values that are present in inventory data, and by including random plot effects in the models.


Introduction
Forest growth models are used to quantitatively generalize forest stand development [1]. A variety of different data sources have been used to develop growth models: permanent plots measured over a full rotation, interval plots measured over one or more growth periods, and temporary plots [2]. Permanent and interval plots may come from inventory or monitoring programs or growth and yield studies that assess silvicultural treatment effects (e.g., thinning, fertilization). Köhl et al. [3] compared permanent inventory and growth and yield study plots for the purpose of growth model development and identified the fact that inventory plots are representative of the total population as major strength for using permanent inventory plots for model development. Inventory plot data have been used to 2 of 16 develop and calibrate equations for growth models [4][5][6][7]. Growth model development is challenging to validate if the data used for modeling are not from a well-controlled growth and yield experiment [1], for which detailed information about site conditions and management history are available. Yet, these experimental data sets are not representative of the total population [3] across a large geographical area and gradients of environmental conditions. Tree-level growth models typically contain equations for diameter increment, height increment, crown recession, and mortality [1]. Tree height and tree volume are closely related [8], and therefore height increment equations are frequently used to predict future tree height in stand development simulations [9][10][11] or to update previously measured tree height in inventories with the ultimate goal to predict and estimate tree volume. Two approaches are commonly used to model height increment [1]: (1) a potential multiplied by modifier approach, where the potential height increment is adjusted by a modifier that accounts for size, site, and competition [12,13] or (2) a linear or nonlinear equation that models realized height increment [14,15]. The use of random effects to account for and take advantage of hierarchical data structure has emerged in the late 1990s [1]. Mixed-effects models allow incorporating locality effects [16,17] and the hierarchical data structure of plot-based forest inventory data into height increment models [18,19].
The Forest Inventory and Analysis (FIA) Program of the United States Forest Service monitors forest resources with a nationwide systematic sample of all forested lands [20]. Due to budgetary constraints, the data collection protocol in the Pacific Coast states was changed in 2015 to only sample a subset of trees for height measurement within each plot [21], which leads to the need of predicting the height of the remaining trees. Currently, the Pacific Northwest (PNW) research station employs a multi-step approach for tree height estimation at a particular time: (1) obtain diameter at both times either by measurement or modeling; (2) predict the height at both times from the diameter using height-diameter equations [22]; (3) estimate the height increment as the difference in the predicted heights at time 1 and time 2; (4) add the estimated height increment to the height at time 1 to estimate the height at time 2. This approach ensures the height at time 2 is always greater than that at time 1. The same method is used in other FIA work units (e.g., Northern and Southern research stations), based on regional calibration of diameter increment models [6,23]. Regional height-diameter equations are often built using plantation-grown trees [24][25][26] that do not represent the total population, and height increment is indirectly obtained as the difference between two heights estimated from these regional height-diameter equations [27]. This approach requires multiple modeling steps that are each associated with model uncertainties [28], and demand future measurements or predictions of tree diameters to obtain height increment [25]. Moreover, this approach does not incorporate plot-or stand-level variation. Direct modeling of tree height increment and a mixed effect modeling approach reduce the steps for future height increment prediction and localize prediction, respectively [27].
Even though height increment models developed from inventory and monitoring data can account for regional and species-specific characteristics by incorporating random effects in the model [18], they may still encounter difficulties in application pertaining to the availability of predictive variables and measurements. Age, one of the tree attributes that is frequently used in height increment modeling, is hard to measure in natural stands [29][30][31] and is not necessarily available from inventory and monitoring data. Forest inventory field measurements often include many unusual data points of tree heights (e.g., negative or zero height increments) due to mensuration errors [32] or any external events such as disease, fire, weather, or mechanical damage [33,34]. However, dropping these unusual measurements from the dataset used for modeling may result in the equations not being representative of the whole population of trees, which is one of the major strengths of using inventory data for growth modeling [3]. Dropping negative or zero height increments would tend to bias estimates of height increment upwards. Therefore, models which allow height increment predictions to all tree records available in the inventory and monitoring database-including trees with atypical height increment observations-are needed. Height increment is defined as the rate of change in height over a specified period of time, which is not the same as height growth, the increase in height of each individual tree Forests 2020, 11, 2 3 of 16 through time [1]. Therefore, height growth is always positive, while height increment could also be negative or zero [35].
The purpose of this study was to develop a tree height increment model that is compatible with continuous regional and national forest inventory and monitoring programs and does not use stand age as predictor variable. The model should be applicable to the inventory data from a broad population of trees with different growth conditions and damage records. For this purpose, we utilized a large sample of FIA inventory plots with 10-year remeasurement interval for ponderosa pine and red alder in forests of Oregon and Washington, USA. The main objectives were to (1) predict realized tree height increment in a single step, thus simplifying the current multi-step approach used; (2) include observations of negative or zero height increment in the modeling process; and (3) achieve more accurate prediction than the current procedure which applies regional height-diameter equations to obtain local height increment predictions. To meet these objectives, we examined six different models that (1) use predictor variables that are commonly measured in regional and national inventories; (2) account for unusual observations (i.e., negative or zero height increment) in the model form; and (3) incorporate random plot effects to improve prediction accuracy by localizing the predictions. We compared model performance through cross-validation.

Data
Forest Inventory and Analysis (FIA) data from Oregon (OR) and Washington (WA) measured between 2001 and 2015 for two tree species-ponderosa pine (PP, Pinus ponderosa Dougl. ex P. Laws. & C. Laws.) and red alder (RA, Alnus rubra Bong.)-were used in this analysis. FIA data contain a spatially balanced sample of inventory plots every 24 km 2 [20]. All live trees over 2.54 cm of diameter at breast height (DBH) were measured within accessible forested plots (i.e., at least 10 percent of potential canopy cover, excluding agricultural and urban areas). We retained plots with five or more trees for modeling and validation purposes. This resulted in 2172 RA trees from 193 plots and 6977 PP trees from 573 plots. Each tree had repeated observations (Time 1 and Time 2) with a 10 year remeasurement interval, except for a few remeasurements with a 9 year interval (i.e., 46 for PP and 26 for RA).
For each tree, the field crews noted if the top was broken and missing, and then (1) recorded the actual height of the tree, and (2) if the top was broken, estimated the height that the tree would have had if the crown was intact [36] (p. 92). Thus, trees with a broken top have a different actual height and field-estimated total height, and were excluded from the analysis. A previously broken top was considered recovered and, therefore, unbroken, if the diameter of the new leader was 1/3 the diameter at the point of breakage. When measuring heights, forked trees were treated in the same way as unforked trees [36]. In our analysis, only the trees with an unbroken top were included in the data. However, the analysis included trees with previously broken and recovered tops, even though the recovered height may be less than that before breakage.
The response variable-∆HT, 10 year height increment in meters-was computed as the difference in heights measured at Time 2 and Time 1 for trees with unbroken tops, as defined in the previous paragraph (i.e., ∆HT = HT 2 − HT 1 ). The height increment ranged from −6.70 to 15.24 m for PP, and from −7.62 to 17.68 m for RA ( Figure 1). Most height increment equations allow positive values only, therefore it is impossible to model negative or zero values of ∆HT with those equations. In our analysis, however, the negative or zero ∆HT accounted for noticeable percentages of the tree data (e.g., 13.21% for PP and 15.75% for RA). We tried to exclude negative observations caused by physical damage or measurement error based on information collected by field crews as damage codes and optional notes in the tally sheet. The attempt of excluding negative height increment values this way, however, only accounted for a small portion of all negative values. When we filtered the data by damage codes that addressed damage to tree tops (e.g., broken or suppressed due to weather or any mechanical injury), 230 out of 922 and 27 out of 342 negative height increment observations were removed for PP and RA, respectively. When we examined the optional field crew notes including the words referencing broken tops or measurement error, additional 42 and 40 negative height increment values for PP and RA were removed, respectively. However, both filtering methods eliminated observations with positive growth as well (e.g., 1196 using the damage code and 822 using the optional field crew notes out of all 7885 positive height increment observations). Because our approach of selectively removing observations based on the available information could not successively withdraw all the observed negative values, we decided to include all data as is and deal with the negative values during the modeling stage.
Auxiliary tree-level variables such as height, DBH, compacted crown ratio, and crown class were selected as potential predictor variables based on existing literature [37][38][39]. Plot-level variables such as site productivity, basal area of all species, elevation, slope, and aspect were also considered to be important predictors [5,15,37]. Full description of the FIA-measured variables can be found in [21]. All tree-level predictor variables and basal area per ha were from Time 1 measurements. Summary statistics of the response and the auxiliary variables obtained from the FIA database are displayed in Table 1.  1 Compacted crown ratio: percent of tree bole with live foliage. 2 Site productivity: capacity to grow crops of industrial wood (potential growth based on the culmination of mean annual increment of fully stocked natural stands). The FIA database classifies site productivity into seven categories based on the potential volume growth [21], and the categories were reclassified in this analysis into three factors.

Models
The current process for estimating tree height increment in the FIA program of the Pacific Northwest (PNW) research station requires DBH measurements at both Time 1 and Time 2, because the tree heights are estimated at both times using the Chapman-Richards height-diameter equation [40] and the difference between the estimated tree heights is added to the height at Time 1. The current method was considered as our baseline model estimating the height increment ( ∆HT base ) over time (t) as follows; where α = 28.26, β = 0.0590, and γ = 1.3040 for red alder (RA), and α = 58.81, β = 0.0123, and γ = 1.1766 for ponderosa pine (PP), respectively [22]. In this approach, the predicted height increment is always positive ( Table 2). The models that we considered as alternatives to this baseline approach needed to accommodate the negative height increment observations in the dataset. Three strategies for dealing with negative response values were taken into consideration: (i) adding a constant to the height increment values to shift the entire distribution of the response to be positive; (ii) setting all negative values to zero; and (iii) using a nonlinear model, which allows negative height increment observations but only predicts positive values. These three strategies are not uncommon in regression analyses that involve non-negative variable transformation or mixed distributions [41,42]. The first strategy requires prior knowledge about the distribution of the response variable, and the second resembles the replacement of missing data strategy [43]. Two models for each of the three options of treating the observed negative height increment values were explored ( Table 2). In order to have positive height increment for all observations, we chose the rounded absolute of the maximum negative value +1 as the constant c (i.e., c = 8 m for PP and c = 9 m for RA) for the log-transformed linear and gamma models. The negative increment observations were also dealt with using a nonlinear exponential model (i.e., ∆ HT = exp(Xβ + Zu)), where X and Z are the fixed and random-effects matrices and β and u are the vectors parameters to be estimated. This resulted in a total of six height increment models that were examined in this study (Table 2). Depending on the model type, different procedures (Table 2) in SAS 9.4 (SAS Institute, Cary, NC, USA) were chosen to fit the models. Except for the baseline model, the models included fixed effects of the tree-and plot-level variables (Table 1) as well as a random plot effect. After fitting the models by species using the full data (i.e., n = 6997 for PP and n = 2172 for RA), variables that were not significant at α = 0.05 were excluded from the model. Significance of the variables was tested using likelihood ratio tests.

Evaluation of Predictive Performance of Models
We performed leave-one-out cross-validation for model evaluation [48], similar to Temesgen et al. [18]. To account for the random plot effect during the process, we conducted the cross-validation as described in Figure 2. The cross-validation process produced one validation dataset accommodating observations from every plot during steps 1-6 ( Figure 2) for each species, and the height increments of the validation dataset were predicted by the six models (Table 2). With 200 iterations of random subsampling and variation in the number of subsamples (m = 1, 2, 3, 4), 800 validation datasets by species were created and used for model performance evaluation.  With the predicted height increment for the validation dataset, we calculated the bias (average difference) and root mean square error (RMSE) at iteration j ( j = 1, . . . , 200) of model k (k = 1, . . . , 6) as follows: where y ijk andŷ ijk are the observed and predicted height increment of the ith tree for model k at iteration j, respectively. The number of plots (l) in the data were 193 and 573 for RA and PP, respectively. Consequently, each validation dataset with the subsample size m had n − lm observations. For the log-transformed linear and the gamma models, y ijk was corrected with the constant added, while the other models used the actual height increment including negative values. Log-bias in the log-transformed linear model was corrected following Duan's smearing estimator [49]. All cross-validation and following analysis were conducted using the SAS statistical software 9.4 (SAS Institute, Cary, NC, USA).

Results
There were 922 observations showing negative or zero height increment for PP and 342 for RA, accounting for approximately 13.21% and 15.75% of the PP and RA data sets, respectively. The observed 10-year mean height increment when negative or zero values were included was 3.3 m for RA, and 1.8 m for PP. Excluding the negative or zero values, the observed 10-year mean height increment was 4.2 m for RA, and 2.2 m for PP. When the negative values were replaced with zero, the mean increment was 3.5 and 1.9 m for RA and PP, respectively.

Model Fitting
Among the eight tree-level and plot-level predictor variables, DBH at Time 1, height at Time 1, crown class, site productivity, and slope were found to be significant in predicting height increment of RA in log-linear, gamma, and quasi-Poisson models (Table 3). Height increment of PP had the same predictor variables, except slope was replaced by elevation (Table 3). Compacted crown ratio also exhibited a strong relationship (e.g., p < 0.0001) with the height increment of PP. In the zero-inflated Poisson model, the production of negative or zero height increment (i.e., zero model) was associated with DBH and height for RA, and DBH, height, compacted crown ratio, and elevation for PP (Table 3). Basal area per ha was an additional predictor variable for RA height increment. In the zero-inflated models, the probability of having negative or zero tree increment was larger for taller trees in both species. DBH of RA, and compacted crown ratio and elevation of PP were positively related to height increment, while they were negatively associated in the zero model. The plot-level random effect was significant for all models (i.e., p < 0.0001).

Comparison of Model Performance
The baseline model substantially overestimated the height increment for both species (e.g., 1.4 m for RA and 0.5 m for PP on average) ( Table 4). The magnitude of mean bias relative to the observed mean height increment was larger for RA (42.81%) than for PP (27.62%) in the baseline model. The mean bias was lowest for the log-transformed linear model for both species across all subsample sizes m, except for the gamma model with m = 1 and m = 2. The average mean bias across subsample sizes was the largest for the quasi-Poisson model (0.20 m for RA and 0.08 m for PP). All suggested models reduced the prediction bias at least by 83% compared to the baseline model (Table 4).  For RA, we obtained lower mean RMSE values than for the baseline model for all suggested models (Figure 3). The magnitude of RMSE reduction was different by model and subsample size. There was a clear tendency for RMSE values to decrease for all models with random plot effects as the number of subsamples included in the modeling set increased (Table 4; Figure 3). When m = 1, RMSE was lowest for the gamma model followed by the zero-inflated Poisson, quasi-Poisson, and then the log-transformed linear model. The gamma model consistently showed the lowest RMSE value across all models, while the log-transformed linear model exhibited a rapid decrease of RMSE as m increased (Figure 3). When m = 4, the quasi-Poisson and log-transformed linear models outperformed the zero-inflated Poisson model. The largest difference from the mean RMSE of the baseline model  For RA, we obtained lower mean RMSE values than for the baseline model for all suggested models (Figure 3). The magnitude of RMSE reduction was different by model and subsample size. There was a clear tendency for RMSE values to decrease for all models with random plot effects as the number of subsamples included in the modeling set increased (Table 4; Figure 3). When m = 1, RMSE was lowest for the gamma model followed by the zero-inflated Poisson, quasi-Poisson, and then the log-transformed linear model. The gamma model consistently showed the lowest RMSE value across all models, while the log-transformed linear model exhibited a rapid decrease of RMSE as m increased (Figure 3). When m = 4, the quasi-Poisson and log-transformed linear models outperformed the zero-inflated Poisson model. The largest difference from the mean RMSE of the baseline model (3.3 m) for RA occurred in the gamma model at m = 4 (2.5 m) with 24.32% reduction in RMSE (Table 4).  For PP, on the other hand, the baseline model had lower RMSE at m = 1 than the other models except for the gamma model. There is a similar pattern of decreasing RMSE with increasing subsample size m as was observed for RA ( Figure 3). Nonetheless, the log-transformed linear model performed particularly poorly with small m values (e.g., m = 1, 2) and only achieved a smaller RMSE than the baseline model for m = 4 ( Figure 3). This was due to extreme observations (i.e., height increment of −6.7 and −6.4 m) that were included in the modeling or validation dataset, but the impact of such unusual observation drastically subsided as m reached 4. All the other models managed to outperform the baseline model with subsample sizes of m ≥ 2. The gamma model performed best in terms of the smallest RMSE for PP as well as for RA (Table 4), and the amount of reduction in RMSE was about 8% of the currently used baseline model (mean RMSE of the baseline model = 1.4 m).
There was large variation in RMSE across the 200 simulation sets for the log-transformed linear and quasi-Poisson models especially for m = 1 for both species, while the gamma model showed less variability in RMSE regardless of m (Figure 4). The large variation of RMSE subsided as the modeling dataset included more subsamples m per plot (i.e., m = 3, 4). For PP, the variability in RMSE was substantial for both m = 1 and m = 2, but leveled-off for m ≥ 3. The variability in RMSE of the baseline model slightly increased as the number of trees in the validation dataset decreased along with increases in m. and quasi-Poisson models especially for m = 1 for both species, while the gamma model showed less variability in RMSE regardless of m (Figure 4). The large variation of RMSE subsided as the modeling dataset included more subsamples m per plot (i.e., m = 3, 4). For PP, the variability in RMSE was substantial for both m = 1 and m = 2, but leveled-off for m ≥ 3. The variability in RMSE of the baseline model slightly increased as the number of trees in the validation dataset decreased along with increases in m.

Discussion
The negative height increment values may be due to measurement errors. When height increments are recorded based on remeasurements on the same tree, the measurement errors can be substantial, especially for slow-growing trees and dense stands with steep terrain [27]. Observed negative or zero height increment values may further be due to poor tree shape or physical and mechanical damages such as topped and partially recovered trees [36]. The proportions of negative height increment values (i.e., 13.21% for PP and 15.75% for RA) were too large to be ignored and to be dropped in the modeling process without substantially biasing the resulting estimators. The mean increment values for both species also differed substantially when calculated with and without the negative or zero values, which further indicated that excluding the negative or zero values from the modeling process may impact modeling results as truncating the data from below will result in overpredictions. This will produce models that are not representative of the whole population of trees [3]. Thus, we decided to explicitly include all the negative or zero observations in the modeling process.

Interpretation of Model Fit and Variables
The general relationship between the height increment and the tree-and plot-level variables was consistent among the species and models. Shorter, fatter, and dominant trees tended to grow more in less steep and lower elevation conditions. The tendencies mostly agree with previous literature characterizing tree height increment in relation with individual tree attributes and competition measures [50][51][52]. Initial height and DBH were among the most important covariates for predicting height increment in previous research [53,54]. Similar to our findings plot or stand-level variables, such as basal area per hectare and site factors (e.g., slope, aspect, elevation, and forest type) have also been found to be highly associated with tree height increment [5,15,37,55].
During data exploration, we found that negative values tended to have higher prevalence in tall trees than in short trees. These tendencies appeared in the estimated coefficients as well. This may be explained by the fact that when measuring standing trees in the field, the trigonometry-based instruments such as hypsometers are more sensitive to errors as the top height is higher, thus there may be more chance of obtaining erroneous height measurement for taller trees [56,57]. Another explanation could be that taller trees have slower height growth than shorter trees [58]. The opposite trends of variables in the fixed effects and in the zero models for the zero-inflated Poisson model suggest that the trees with negative or zero growth may possess a different data-generating process than the positive increment values. Thus, zero-inflated or hurdle models [59], which allow modeling the negative or zero increment observations separately from the positive increment observations may be useful in modeling tree height increment when the data include negative or zero increment [60].
The significant plot-level random effects across all models suggest that the inclusion of a random plot effect contributed substantially to the reduction of random variance [18]. This was also confirmed by the reduction in RMSE as subsample size m increased, because a larger number of subsamples m included in the modeling dataset improved the estimation of the random plot effect. Previously shown gains in precision of height prediction by inclusion of random effects [18] support our decision of incorporating the plot-level random effect to account for the hierarchical data structure of the inventory data, thus allowing localized height increment predictions. As subsamples of height are typically measured on regional and national inventory and monitoring plots, random effects can be estimated from existing databases during future model applications when height increment is predicted for trees without height measurements.

Comparison of Model Performance
Differences in the relative mean bias between RA and PP (e.g., 42.81% vs. 27.62%) imply that the current height diameter equation [22] may be less accurate for RA than for PP. It is generally considered that height measurements tend to be challenging for irregular or round-crowned trees, hence height increment models are more difficult to develop for hardwood species than for conifers [1]. Average absolute gains in RMSE were smaller for PP (0.04-0.10 m) than for RA (0.60-0.83 m) for the gamma model that performed best among any of the suggested models. At least for RA, this may be a considerable improvement over the baseline model, as RA trees are much shorter than PP trees (e.g., Hanus et al. [61]; average height of 25.6 m for PP vs. 9.2 m for RA). Differences in performance among the suggested models themselves were quite small, especially for larger subsample sizes (m ≥ 3; approximately 0.15 m).
Based on the previously mentioned small differences in performance among the suggested models, it may not matter too much which model form is chosen as long as the model accounts for negative or zero observations in some way and is not too sensitive to extreme observations in the data set. Especially, the log-transformed linear model was shown to be sensitive to extreme values when small subsample sizes were used. Changyong et al. [62] demonstrated that log transformation can result in more variable and more skewed data through simulation. Another model performance characteristic to consider when selecting the model may be the variability of the model across subsample sizes m. The gamma model, which had the lowest RMSE for all subsample sizes, demonstrated very low variability in performance across subsample sizes, which provides another advantage of this modeling approach. Similarly, the nonlinear model showed little variability across the subsample sizes. On the contrary, the zero-inflated and quasi-Poisson models exhibited high variability, especially for smaller subsample sizes. This may suggest sensitivity of these models to extreme values when only small subsample sizes are available [63], and makes these two models less desirable to use for height increment predictions.
Zero-inflated and quasi-Poisson models adopted the second strategy of negative value treatment (i.e., set negative values to zero), inevitably showed larger mean bias, because the response values were modified during the modeling process (Table 2). However, the differences in bias across the models were minor, and all the suggested models managed to reduce the prediction bias at least by 83% compared to the baseline model. Although this result is possibly also due to the fact that the suggested models were built from the same dataset used for the cross-validation while the baseline model was not, the improvement in the accuracy of prediction is noteworthy. Acquisition of a new dataset to validate the model will confirm the true performance of our suggested models contrasted with the current method of height increment prediction.
The quasi-Poisson model exhibited the largest mean bias but smaller RMSE, showing typical variance-bias trade-off in model prediction [64]. It is noteworthy that for the shifted models with the constant-added response (e.g., log-transformed linear model and gamma model), the constants should be determined before the modeling process. We chose the constants from the observed maximum absolute value of the negative height increment, c = 8 m for PP and c = 9 m for RA, but these values may vary with different datasets. Considering the need of prior information to find the constant to be added as a drawback of the gamma model, we recommend the use of the quasi-Poisson model for future height increment predictions based on its second-lowest RMSE value across the two species examined in this study.
The zero-inflated Poisson model provides the ability to explicitly model negative and zero height increment values separately from the positive values [65]. Therefore, this model may be able to capture the mechanism that creates the negative or zero height increment values, provided any of the available predictor variables have the explanatory power to do so. We also implemented a zero-inflated negative binomial model during the model fitting and validation process. A zero-inflated negative binomial model would similarly to the zero-inflated Poisson model allow capturing the mechanism that creates the negative or zero values [65]. However, the results were not reported in this paper because the differences in performance between the zero-inflated negative binomial and the zero-inflated Poisson were minimal. However, the zero-inflated negative binomial model required excessive time for model fitting with low convergence rate compared to the zero-inflated Poisson model.
The main improvement for all models appeared to be for larger subsample sizes (m ≥ 3). This may suggest that the inclusion of the random effect, which improves with larger subsample sizes [18], may have provided the most improvement over the baseline model, which does not include a random effect and is therefore unable to localize height increment predictions. The inclusion of a random effect to allow localized predictions may be the most important strategy in improving upon the currently used baseline model. It may be possible to include random plot effects in the currently used height prediction models [22] that are used in the baseline model. This could lead to similar results for height increment predictions as were achieved by the gamma and the other models that allowed localized height increment predictions.

Height Increment vs. Height Growth Models
From an inventory and monitoring point of view, it is necessary to predict height increment for trees for which heights were not measured. The main interest for inventory and monitoring applications is therefore to estimate changes in measured height in the population of all trees, under the same measurement conditions as those of the measured trees. Therefore, the focus is strictly on height increment prediction, where increment is defined as the rate of change in height in a specified period of time [1]. On the other hand, height growth refers to the actual increase in height of each individual tree through time [1]. Height growth can be considered a latent process [66], which contributes to the observed differences in height increment among other possible factors (e.g., diameter increment, crown recession, and mortality). Using height increment models to infer height growth assumes that a population-level change is the same as growth at the individual tree level. This is especially a big assumption when negative or zero height increment values are included in the modeling process as suggested in this study for predicting height increment in the framework of inventory and monitoring databases.

Conclusions
Height increment models are major components of growth and yield models and are also needed to update inventory and monitoring databases when tree heights are only subsampled. Inventory data contain many negative or zero height increment observations that cannot be explained by damage records. Therefore, these negative or zero height increment values are part of the total population of trees, as measured using the inventory field protocols. Yet, when developing height increment equations from inventory and monitoring databases, negative or zero height increment observations are often dropped from the modeling data set, which results in models that may not be fully representative of the total population of tree measurements. Therefore, negative or zero height increment observations should not be ignored in model development.
In this study, we demonstrated that the currently used baseline method of the FIA inventory in the Pacific Northwest region-predicting height increment as difference in modeled heights-can be improved upon by using random effects models that explicitly account for negative or zero height increment observations in the modeling process and include random plot effects to allow localized predictions. We employed a nonlinear model, two models that added constants to the negative observations to shift data into the positive range, and two models that set all negative values to zero.
Models for red alder (RA) improved the baseline model substantially, while the reduction in RMSE was not that large for ponderosa pine (PP) models. The gamma model performed best for both species across all subsample sizes (m = 1, 2, 3, 4), followed by the quasi-Poisson and non-linear models for larger subsample sizes (m ≥ 3). Gains in model performance over the baseline model were more pronounced for larger subsample sizes.
The somewhat moderate observed gains in model performance, especially for PP, may not justify the need to fit more complicated models. As observed differences between the baseline model and models that account for negative or zero observations and include a random effect were not that large, we suggest that models are used that are most convenient for a specific application. The inclusion of random effects to localize predictions are suggested for any model and could improve the currently used baseline model as well.