Influences of mean top height definition and sampling method on errors of estimates in New Zealand ’ s forest plantations

Background: A study was undertaken of 51 stand inventories to compare two alternative mean top height (MTH) calculation methods prevalent in New Zealand, and to evaluate the consequences of creating height versus diameter at breast height (H-D) curves at a stand-level during inventories as opposed to fitting H-D curves at a plot-level. Methods: The dataset was separated into two groups; one with plots having less than 6 heights measured and one with more than 5 heights measured. MTH was calculated using all possible combinations of the two calculation methods and with H-D curves either at a stand-level or a plot-level. Graphs were prepared to compare the 4 alternative MTH estimation techniques for all plots. In addition standard deviations of MTH between plots were calculated within stands, and then these were compared for different MTH calculation methods using interleaved histograms and with a mixed effects analysis


Introduction
Mean top height (MTH) is an important variable estimated in inventories and permanent sample plots (PSPs) in New Zealand (Goulding 2005), and it is a critical variable for the estimation of site productivity through site index.The relationship between height growth and productivity has long been recognised (Bauer 1881), and MTH at age 20 is commonly used as an index of site productivity of Pinus radiata D.Don. in New Zealand (Goulding 2005).
Site index has been used as a final test calibration for eco-physiological models of site productivity in Sweden (Mason et al. 2017), and New Zealand plantation owners are keen to access similar eco-physiological estimates for their forests calibrated with site index.Calculations of site index using MTH and age estimates from PSPs and forest inventory plots have been proposed for calibrating site productivity models at high resolution across their estates, but inventory estimates and PSP estimates differ in two potentially important respects.
Firstly, MTHs in inventories and PSPs are typically calculated in different ways; the Forest Research Institute (FRI) method (Goulding 2005), and the Carter Holt Harvey (CHH) method (Woollons 2003).
The FRI method involves fitting a height-dbhob (H-D) curve using a Näslund equation (Equation 1) (Näslund 1936) to all trees with heights and diameters at breast height outside bark (dbhobs, generally measured at 1.4 m above ground level in New Zealand), and then dbhobs of trees representing the largest 100 dbhobs/ha are averaged to determine "mean top diameter" (MTD).Finally, the MTD is inserted into the Näslund equation to determine MTH.Several New Zealand references refer to the Näslund equation as the "Petterson equation", but although Petterson (1955) reported on it, Näslund, who worked with Petterson, invented it and so it should always bear his name.
The CHH method also involves fitting a Näslund equation to all trees with heights and dbhobs measured, but then the equation is used to estimate heights of unmeasured trees among those representing the largest 100 dbhobs/ha, and the heights, both measured and estimated, of trees representing the largest 100 dbhobs/ ha are averaged to determine MTH.H=1.4+ (b+ a/dbhob) (-2.5)  (1) where H = estimated height, and dbhob = diameter at breast height outside bark.Coefficients a and b can be estimated using simple linear regression techniques following a transformation (Equation 2).

Y=a+b* dbhob (2)
where Y = dbhob/(H−1.4) 0.4 and other variables are as described previously.Woollons (2003) compared the two MTH estimation techniques in experimental plots with all trees' heights and diameters measured and found that under those circumstances they were almost identical.However, if only a few heights were measured in each plot and H-D curves were fitted at a stand-level rather than a plotlevel, the two techniques might yield different estimates of MTH, and this has not been tested previously.
Secondly, MTHs of plots estimated during inventories can lack independence, while those in PSPs do not.Inventory estimates in New Zealand generally lack independence because only a very few heights (often 2-3 and in some plots none) are measured in each plot, and so H-D curves are fitted at a stand-level rather than one curve per plot.PSPs, on the other hand, generally contain a minimum of 12 trees with heights measured, chosen to cover the range of dbhobs in a plot (Ellis & Hayes 1997), and so H-D curves can be fitted independently to each plot.
Preliminary investigations in New Zealand have suggested that PSP estimates of MTH might yield more diverse and somewhat different calibrations of ecophysiological estimates of productivity than calibrations with inventory MTH estimates.In addition, when inventory MTH estimates were plotted on maps, there often appeared to be very little variation within stands but large variation between adjacent stands (Figure 1).
Inventory plots are generally much more numerous than PSPs, however, and so it would be good to be able to use inventory estimates for calibrations.This led to the question of how much these two sources of MTH estimates might differ, and also to whether or not forest managers might be under-estimating variation between plots in their inventory samples.
Rayonier (NZ) Ltd. conducted 51 inventories in 2017 with more heights measured per plot than usual, and the large inventory dataset enabled a study to answer the following three questions: 1. Is there a substantial difference between inventory estimates of MTH when H-D curves are fitted within plots compared to when they are fitted within stands?
2. When only a few heights are measured in each plot and H-D curves are fitted at a stand level, does the FRI MTH calculation technique differ substantially from the CHH technique?
3. Does the lack of independence in inventory estimates of MTH result in an underestimate of variation between plots within stands?

Methods
Fifty-one inventories were conducted in stands of Pinus radiata D.Don in southern New Zealand in a forest estate of slightly less than 44,000 ha planted with approximately 20% open-pollinated and 80% controlled-pollinated nursery stock.In total there were 865 plots distributed on randomly located grids throughout the stands.This arrangement of plots has been recommended as one that provides unbiased estimators of crop condition (Gordon & Pont 2015).Features of the plots are shown in Table 1.

TABLE 1: Features of plots in the inventories
In each plot all the dbhobs were measured with a diameter tape, and a subset of trees' heights (as indicated in the table) were measured with a Vertex hypsometer.The distribution of numbers of height measurements per plot was bimodal, with peaks at 3 and 7 heights per plot (Figure 2).The dataset was separated into two groups; one with plots having less than 6 heights measured and one with more than 5 heights measured.MTH was calculated in the following four ways (where numbers of height measurements allowed): 0. 1.
H-D curves fitted at a stand level and MTH calculated using the CHH method 2.
H-D curves fitted at a stand level and MTH calculated using the FRI method 3.
H-D curves fitted at a plot level and MTH calculated using the CHH method 4.
H-D curves fitted at a plot level and MTH calculated using the FRI method Graphs were prepared to compare the four alternative MTH estimation techniques for all plots with more than 5 heights measured per plot.In some cases plots were also examined using plots with <6 and >2 heights measured.In 15 cases, where heights measured within stands were few, MTH calculations at a plot level were unreliable owing to illogical curve coefficients, and so these estimates were excluded from the study.
Standard deviations of MTH were calculated within stands, and then these were compared for different MTH calculation methods using interleaved histograms.Comparing standard deviations between methods using within-plot H-D curves and within-stand H-D curves involved comparisons using H-D curves from differently sized samples.A further calculation of MTH and standard deviation was conducted using height and dbhob measurements from within stands that had a similar median and distribution to those within plots where numbers of heights measured was greater than five.This latter step was accomplished by firstly randomly selecting only one measurement per plot and then randomly adjusting the numbers in each stand to be within an appropriate range.The median in each case was 7 heights measured.Comparing these standard deviations with those from MTH calculation methods involving H-D curves would therefore indicate the extent to which using a common H-D curve for all plots on stand would influence variation between MTH estimates between plots in a stand.Standard deviations were transformed (λ=0.53) using a Box-Cox (Box & Cox 1964) method (Equation 3) to make their frequency distribution as normal as possible.A mixed effects model was then fitted to examine the joint effects of calculation method and level of the H-D curve (Equation 4).The level of H-D curve had three categories: plot-level, standlevel using all heights, and stand-level using a random selection of heights with a median of seven per stand.Stand was a random effect in this analysis.
Where sd = transformed standard deviation of MTH between plots in stand, P = level of H-D calculation, M = method of MTH calculation, S = random effect of stand, and ϵ = random error.
In order to further explore reasons why variation of MTH estimates within stands might be influenced by a common H-D curve, MTH and stems/ha estimated from each plot were divided by stand average MTH and stems/ha respectively calculated from all plots, to

Feature
Min. Median Mean Max.
Years The dataset used for the study contained a range of plot sizes (Table 1), and this can cause a bias in estimation of mean top height (Garcia 1998;Garcia & Batho 2005;Magnussen 1999;Ochal et al. 2017).Garcia (1998) suggested that if MTH was defined as the mean height of the largest dbh tree in each 0.01 ha of a stand as advocated by Rennolls (1978) then in plots where tree positions have not been defined a U-statistic estimator can be employed to provide a nearly unbiased estimate of mean top height defined in this way.The U-statistic estimator is the height of all trees weighted by the number of times they appear as the largest tree in all possible mean top-sized groups of trees in a plot.This procedure has since been employed successfully (Garcia & Batho 2005;Ochal et al. 2017).To check whether or not the conclusions of the study reported here were affected by bias due to varying plot size, two alternative definitions were developed that were U-statistic estimator analogues of the FRI and CHH definitions used in New Zealand, and the analysis was repeated using these two definitions.The definitions used either all Näslund curve estimates of height (the FRI analogue), or measured estimates when possible and otherwise curve estimates (the CHH analogue).Neither of these definitions is currently employed in New Zealand, and so they were used simply to check on the validity of the study's conclusions.
All calculations, analyses, and geographical Information system operations were done in R software (R Development Core Team 2004).

Results
Comparing the CHH and FRI calculations with H-D curves at a plot level, in the same manner as Woollons (2003), replicated Woollons' result showing that two methods appeared to be equivalent for practical purposes (Figure 3).However, the two methods gave different results when compared with H-D curves estimated at a stand level, and the range of estimates MTHs was noticeably smaller (Figure 4).There was also a hint that estimates for at least one of the methods might be biased for smaller MTHs.When fewer than six heights were measured in each plot, the differences between the two methods became even more pronounced (Figure 5).For both calculation methods, the stand-level and plot-level estimates differed substantially (Figures 6 and  7), but the divergence was greater for the FRI method.In addition, there was evidence of bias at the extremes of these plots.FIGURE 6: A comparison of MTH calculation methods using stand-level and plot-level H-D curves in plots where at least six heights were measured per plot, using the FRI calculation method.The line shows where two estimates would be exactly equivalent.Histograms of standard deviations by MTH calculation method and level of H-D equation showed that there was a tendency for stand-level equations to have smaller standard deviations than plot-level equations and that MTH calculation method influenced standard deviation when stand-level methods were applied (Figure 8).
The mixed effects analysis of standard deviations with stand as a random effect showed clearly that the interaction between H-D curve level and MTH calculation method was statistically significant (P<0.0001).Furthermore, standard deviations of stand level curves calculated with only 7 heights did not differ significantly from those calculated using all heights in a stand (Figure 8).
Relative MTH tended to show a negative slope (P=2.2e-16) when plotted against relative stems/ha using within-stand H-D curves, and this trend more evident when using the FRI calculation method (Figure 9) than when using the CHH one.By contrast, using plot-level H-D curves exhibited a slightly positive trend (P=2.358e-06)(Figure 10).Repeating the analysis with contrasting definitions of mean top height employing the definition of "mean top" suggested by Rennolls (1978) and U-statistic estimators developed by Garcia (1998) yielded very similar results to those obtained with the New Zealand definitions.The magnitudes of differences in standard deviation shown in Figure 8 were slightly less when using the U-statistic estimator MTH definitions.For instance instead of the 69% reduction in standard deviation observed between plot-level and stand level Näslund curves for the FRI technique, the reduction was 60%.All patterns and trends shown in the Figures for the New Zealand definitions of MTH were the same when using the two U-statistic estimator definitions.

Discussion and Conclusions
Independence of sampling units is a fundamental assumption for the use of statistics such as t to calculate confidence limits around sample estimates ('Student' 1908;Fisher & MacKenzie 1923)."Independence of sampling units" means that sampling units are not related; that the value obtained from one sampling unit is unaffected by the value obtained from another.Using heights and dbhobs from other sampling units to formulate an H-D relation for MTH calculation clearly violates this assumption.If sampling units lack independence, then standard errors may be underestimated, leading to overestimates of precisions of sample estimates.In the study reported here, using stand-level H-D curves reduced the width of estimated confidence intervals by 69% using the FRI MTH calculation method, and by 29% using the CHH method, with a median of seven height measurements per plot.The reason for this is clear, H-D curves obtained independently from plots have different levels on an H-D graph, and using only one "average" H-D curve across a stand constrains the range of possible estimates of MTH (Figure 11).The CHH method gives greater weighting to local height measurements, thereby reducing the extent of lack of independence compared to the FRI method.Measuring a median of only three height measurements/plot, as is common in New Zealand's plantation inventories, would cause the CHH method to underestimate precision of estimated MTH to a greater extent than outlined in the study reported here.300 index calculations (Kimberley et al. 2005) would also be affected by lack of independence, because site index is estimated before solving a basal area/ha function to establish a 300 index estimate.An alternative explanation for greater variation among plot-level H-D estimates was that the larger number of measured heights and dbhobs used to make a stand-level H-D curve might provide a more accurate H-D equation, but this possibility was shown to be of little consequence by the comparison between standlevel curves using a median of 7 heights and dbhobs and plot-level H-D curves constructed with the same number of measured heights and dbhobs (Figure 8).However, stochastic simulations of samples from a large population with linear correlations similar to those between height and dbhob suggest that a median of seven heights is barely adequate (results not shown), and that Ellis & Hayes (1997) are right in recommending at least 12 height measurements in PSPs to get more precise estimates of MTH.More height measurements per plot is likely to reduce variation in MTH between plots, but not to the extent that the variation is as small as variation between MTHs estimated from a stand-level H-D curve (Figure 11).A future study should examine the extent of variation in stands with all heights measured, perhaps by LiDAR (Saremi, Kumar, et al. 2014;Saremi, Lalit, et al. 2014) and compare alternative sampling strategies when all heights are known.
Stems/ha varied enormously within stands (Figure 9), to the point where locally high stocking led to lower estimates of MTH with stand-level H-D curves compared to estimates of height in areas with locally low stocking.A stand-level H-D curve would tend to wrongly indicate that trees should be shorter where locally high stocking led to smaller diameters, but when plot-level H-D curves were used the reverse was true, albeit to a lesser extent (Figure 10).The latter result may reflect a small increase in mean and mean top heights with stems/ha, as has been reported in previous studies (Maclaren et al. 1995;Mason 1992).In addition, having more trees to sample from may cause a small increase in the heights of the largest 100 stems/ha and hence MTH when a local H-D curve is employed.Stems/ha estimates from each plot may be useful as a covariate in estimating MTH if this effect is characterised across a wide range of stands.
The effect of stand-level H-D curves in reducing estimates of MTH variability very likely explains the bias observed at extremes in Figures 6 and 7.This challenges our definition of the word "stand" as an area of land with more or less homogeneous site quality, forest species composition, stand structure and age structure, unless we give a great deal of latitude to the phrase "more or less".
The results of this study confirmed the findings of Woollons (2003) that when H-D equations were estimated at a plot level, the CHH and FRI MTH calculation methods differed little in their respective estimates (Figure 3).However clearly the two methods differed substantially when H-D curves were formed at a stand level (Figure 4).This latter finding can be explained by the fact that using a stand-level H-D curve causes estimates of MTH to be less independent with the FRI calculation method than with the CHH method.
Implications of these findings for forest managers are that: Estimates of stand-level MTH and volume will be much less precise than statistical calculations indicate if H-D equations are formed at a stand-level.When using the CHH MTH calculation method this effect will become larger as numbers of heights measured per plot is reduced.
Studies designed to evaluate site productivity within stands for the purposes of precision forestry will be misled by site indices derived from MTHs where H-D curves were formed using stand-level data, especially if a median of only 3 height measurements are obtained per plot.
Consideration should be given to increasing the numbers of heights measured in each inventory plot in order to allow plot-level H-D curves to be calculated.In New Zealand PSPs are established with at least 12 height measurements/plot (Ellis & Hayes 1997), which may seem too great an expense for inventories, except that measuring larger plots may reduce overall inventory costs because fewer plots may be required for any given level of inventory precision (when precision is calculated properly, from independent sampling units).Ellis & Hayes recommend that PSPs should be large enough to accommodate at least 20 trees that will be in the final crop, and median sizes of PSPs are generally around 0.1 ha, compared to a median of 0.022 ha for the inventory plots used for this study.Optimising inventory design should be a priority, given the new findings described here.Trees for height measurement should be chosen across the range of dbhobs in a plot, including the largest and smallest dbhobs, but with more weighting to larger trees if MTH estimates are the objective.
Forest managers should consider creating a consistent definition of MTH.As it is more robust when sampling units lack independence, the CHH method is recommended.However, sampling with truly independent sampling units should be the preferred option.Previous studies have found that the "mean top" criterion is best established as an absolute number of trees per hectare rather than a proportion (Ritchie et al. 2012), and in order to reduce dependence of MTH estimates on plot size, the British Columbian definition based on the mean of largest trees in 0.01 ha plots (Garcia & Batho 2005), estimated with U-statistic estimators in larger plots (Garcia 1998), should also be considered.The practicality of these alternatives in New Zealand conditions will be the subject of future studies.

FIGURE 1 :
FIGURE 1: Site indices (MTH at age 20) for Pinus radiata D.Don estimated by standard inventory procedures across part of a forest estate, coloured by site index.Each point represents a plot location, and those within the same grid and at the same spacing are in the same inventory.The lines are forest boundaries.Typically 1-3 heights were measured per plot, and all H-D curves were estimated at a stand level.Note how different inventories show markedly different site indices for almost the same points, and how little variation in MTH estimates exist within inventories.

FIGURE 2 :
FIGURE 2: Frequencies of height measurements per plot.

FIGURE 3 :
FIGURE 3: A comparison of the FRI and CHH MTH calculation methods using plots where at least six heights were measured per plot with plot-level H-D curves.The line shows where two estimates would be exactly equivalent.

FIGURE 4 :
FIGURE 4: A comparison of the FRI and CHH MTH calculation methods using plots where at least six heights were measured per plot with stand-level H-D curves.The line shows where two estimates would be exactly equivalent.

FIGURE 5 :
FIGURE 5: A comparison of the FRI and CHH MTH calculation methods using plots where less than six heights were measured per plot with stand-level H-D curves.The line shows where two estimates would be exactly equivalent.

FIGURE 7 :
FIGURE 7: A comparison of MTH calculation methods using stand-level and plot-level H-D curves in plots where at least six heights were measured per plot, using the CHH calculation method.The line shows where two estimates would be exactly equivalent.

FIGURE 8 :
FIGURE 8: Mean standard deviations of plot MTH estimates within stands versus level of H-D curve ("Plot" = H-D at a plot level, "Stand" = using all heights in a stand, and "Stand7" = using a median of 7 heights/stand), and MTH calculation method (CHH or FRI).Bars with the same letter are not significantly different.

FIGURE 9 :
FIGURE 9: Relative MTH versus relative stems/ha, using stand-level H-D curves in plots where at least six heights were measured per plot, using the FRI calculation method.The line shows a fitted linear regression (P=2.2e-16).

FIGURE 10 :
FIGURE 10: Relative MTH versus relative stems/ha, using plot-level H-D curves in plots where at least six heights were measured per plot, using the CHH calculation method.The line shows a fitted linear regression (P=2.358e-06).
FIGURE 11: H-D curves at a plot-level (blue), and stand-level (red) for one example stand.Points show MTH estimates, using the FRI calculation method.