Wood density estimates of standing trees by micro- drilling and other non-destructive measures Christine L. Todoroki1*, Eini C. Lowell2 and Cosmin N. Filipescu3 1 Scion, Titokorangi Drive (formerly, Longmile Road), Whakarewarewa, Rotorua 3010, New Zealand 2 USDA Forest Service, 620 SW Main St. Suite 502, Portland, OR 97205, USA 3 Canadian Wood Fibre Centre, 506 Burnside Road West, Victoria, BC, V8Z 1M5, Canada *Corresponding author: christine.todoroki@scionresearch.com (Received for publication 11 October 2019; accepted in revised form 7 April 2021) Abstract Background: Accurate estimates of wood density are needed by the forest sector to increase value along the tree-to-product value-chain. Amongst tools supporting in-situ assessments, micro-drills and acoustic hammers have become increasingly popular. Our objective was to use these tools, and other easily-obtained measures, to develop predictive wood density models for in-situ assessments of Douglas-fir (Pseudotsuga menziesii (Mirb.) Franco) trees in western North America. Methods: Wood density estimates of 133 trees, 60–75 years-old, were benchmarked against X-ray densitometry data using linear mixed-effects models. Mean resistograph amplitude (unadjusted, adjusted, and standardised variants), and combinations of acoustic velocity, tree diameter, stand age, and site index were considered as fixed effects. Plots, comprising differing treatments, and sites were considered as random effects. Candidate models were selected based on fit statistics, and further evaluated with an independent external dataset comprising 37 Douglas-fir trees. Results: The optimal model comprised amplitude (adjusted), site index (transformed), and the quotient of velocity and age. It had a mean absolute percentage error, MAPE, of 4.1%, mean absolute error, MAE, of 19.4 kg.m-3, a root-mean-squared- error, RMSE of 25.0 kg.m-3, and marginal R2 for fixed effects, R2marg of 0.60. With external data, MAPE was 8.7%, MAE 52.4 kg.m-3 and RMSE 59.5 kg.m-3. Fit statistics for a simpler two-variable model (standardised amplitude and transformed site index) were: MAPE 4.9%, MAE 23.2 kg.m-3, RMSE 28.0 kg.m-3, and R2marg, 0.48, and with external data MAPE was 8.5%, MAE 51.6 kg.m-3 and RMSE 59.3 kg.m-3. Thus, with external data, the simpler model produced greater accuracy than the optimal model. Amplitude, and all other single-variable models, recorded poorer levels of accuracy. Conclusions: Micro-drilling alone, though highly significant as a predictor, is insufficient for providing accurate wood density estimates of individual trees. Site effects need to be considered too. Standardisation of mean amplitudes to z-scores makes models highly portable across a range of resistance tools and operating speeds, and therefore more practical. As noted in the literature, optimal models are not necessarily best for predicting outcomes with other datasets, therefore model evaluation with external data is critical to determining how well a model will perform in practice. New Zealand Journal of Forestry Science Todoroki et al. New Zealand Journal of Forestry Science (2021) 51:6 https://doi.org/10.33494/nzjfs512021x74x E-ISSN: 1179-5395 published on-line: 22/06/2021 © The Author(s). 2021 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Research Article Open Access introduced tree species in New Zealand (Wang et al. 2001; Maclaren 2009), France, Germany, the UK, Spain, Belgium, and the Czech Republic (Zeidler et al. 2017; Spiecker et al. 2019). Its popularity arises from its relatively fast growth, high productivity, and desirable wood properties. Wood density is one such desirable property and an important indicator of the quality of solid wood products Introduction Douglas-fir (Pseudotsuga menziesii (Mirb.) Franco) is widely recognised as being the most commercially significant species for structural applications in North America, particularly the Pacific Northwest of the United States, and the west coast of Canada. Outside of North America, Douglas-fir is amongst the most commonly Keywords: Resistograph; acoustic velocity; Douglas-fir; mixed-effects models http://creativecommons.org/licenses/by/4.0/), Todoroki et al. New Zealand Journal of Forestry Science (2021) 51:6 Page 2 (Zobel & Van Buijtenen 2012). Wood density directly influences and is correlated with stiffness and strength (Zeidler et al. 2017), thus is used as a selection criterion in tree breeding programs for improving yield of high quality structural lumber (Howe et al. 2006). Because Douglas-fir wood density varies considerably from one tree to another, from one site to another, and from one region to another (Filipescu et al. 2014; Kimberley et al. 2017), forest managers and wood processors need to be able to assess this wood property to best match the raw material to the final product. When assessed in-situ, there are greater opportunities to increase revenue to both land and mill owners, through improved selection, sorting, allocation to the correct processing pathway, and through reduced wastage and manufacturing costs. Therefore, there is a need for rapid, accurate, and non- destructive methods for assessing wood density in standing trees. A range of non-destructive tools and methods have been developed for in-situ assessments of wood density (Gao et al. 2017; Schimleck et al. 2019). Examples include wood increment borers, the Pilodyn Wood Tester, torsiometers, and micro-drills (Gao et al. 2012; Wessels et al. 2011). Acoustic tools, although developed primarily for assessing wood stiffness, have also been used in wood density studies (e.g. Chauhan & Walker 2006; El-Kassaby et al. 2011; Newton 2017). Increment cores, when analysed using X-ray densitometry techniques (Walker & Dodd 1988; Eberhardt & Samuelson 2015), enable very accurate measurements of wood density. However, X-ray densitometry is expensive and time-consuming (Chantre & Rozenberg 1997), requires rigorous methods for preparing and processing cores, and therefore is neither a rapid technique nor can it be applied in-situ on standing trees. Other tools such as the Pilodyn and torsiometer, which penetrate only a small way through the bark, are less destructive than increment borers, but, according to Gao et al. (2012) they cannot be considered as substitutes for the increment borer as they have had only limited success. In contrast, micro-drill resistance tools have demonstrated greater potential (Isik & Li 2003; Gao et al. 2012), are less invasive than increment borers due to the smaller holes (3 mm) left following drilling (Rinn 1988) and are used extensively in assessing progeny trials (Bouffier et al. 2008; Gwaze & Stevenson 2008; El- Kassaby et al. 2011; Desponts et al. 2017). Micro-drill resistance tools record the amplitude of the resistance to turning (torque) experienced by a fine drill when driven through wood at a given forward speed (cm per min) and rotational frequency (rpm) (Rinn 1988). As the bit progresses through the stem, resistance due to friction generally increases, thus creating an increasing trend. Therefore, adjustments need to be made to the data prior to analyses, to remove any potential sources of bias. Methods developed to adjust, or detrend, resistance profiles include trigonometric approaches (Gantz 2002; Fundova et al. 2018), smoothing functions, and translation functions that shift the baseline to correct for bias (Isik & Li 2003; Eckard et al. 2010; Fundova et al. 2018). Due to the nature of these detrending methods, the full bark-to-bark profile is required. Shorter profiles (e.g. the first 5 cm of the inside-bark profile) have also been evaluated (Bouffier et al. 2008) and moderate relationships (i.e. R2 ≈ 0.41 and 0.48, based on correlations of 0.64, 0.69) reported between mean (adjusted) amplitude and mean wood density assessed in progeny trials of maritime pine (Pinus pinaster Ait.) at two sites. In general, only weak to moderate relationships have been found for individual trees. Gwaze & Stevenson (2008) reported an R2 of 0.23 for 25-year-old shortleaf pine (Pinus echinata Mill.) in Missouri, USA, while Walker et al. (2019) reported an R2 of 0.47 for 6-9-year-old loblolly pine growing in southeastern US and noted improved correlations with inclusion of site effects in their predictive models. Isik and Li (2003) obtained R2 values of 0.21, 0.24, 0.31, 0.44 for 11-year-old loblolly pine (Pinus taeda L.) at four sites in South Carolina, USA, and when all sites were combined, the phenotypic correlation was just 0.12. Phenotypic correlations, often estimated using product-moment correlation statistics (also called Pearson’s r) should not be confused with genetic correlations that, based on correlations between family means, will always be stronger than their phenotypic counterparts. It should also be noted that the coefficient of determination, R2, equal to the square of r, describes the explanatory power of a model with the dataset at hand, not the predictive ability or usefulness of the model to new data obtained from other settings. Therefore, models that fit well to in-sample data (i.e. with a high R2 value), may not necessarily provide accurate predictions when applied to new data (Mendenhall & Sincich 2012). The same holds true for other variates of R2 including marginal and conditional R2 values, R2marg and R 2 cond respectively, (Nakagawa & Schielzeth 2013) where R2marg represents the proportion of variance explained by fixed effects alone (akin to R2 in a fixed-effects model) while R2cond represents the proportion of variance explained by both fixed and random effects. A low R2marg and high R2cond implies that the fixed effects (measured variables) explain little of the variance while random effects (unmeasured) explain a far greater proportion of the variance, indicating that there are other factors that the model may have failed to capture. In addition to examining goodness-of-fit statistics, the use of an independent dataset is considered the “gold standard” for assessing the predictive power of models (Alexander et al. 2015). Moreover, this new dataset adds a further level of scrutiny to the model construction process (Snee 1977). When it is not possible to collect new data, techniques such as data-splitting, cross- validating, or bootstrapping can be applied (Snee 1977; Dankers et al. 2019). Despite the results of previous studies demonstrating the potential of micro-drilling to estimate wood density, a factor that has limited its practical use is the availability of accurate, reliable models that make predictions with new data from other settings. The objective, and key challenge of our research, was to develop individual tree-based predictive models and to access their accuracy with both in-sample and external datasets. Using our predictive models, we wished to resolve three key questions: 1) Can micro-drill resistance tools alone provide accurate assessments of wood density for a diverse set of trees? 2) When used in combination with acoustic velocity tools, and other easily measured variables (diameter, stand age, site index), what level of accuracy can be achieved? 3) How portable are the predictive models to new data? Overall, we wanted to develop robust, portable models for rapid in-situ assessments of wood density in individual standing trees. Methods Study sites The Douglas-fir trees of this study were located in coastal western North America. They ranged in age from 60 to 75 years old and were selected from six experimental installations, four of which were planted, and the remaining two of natural origin, both of which had regenerated naturally after wildfire. Site index (at 50 years) ranged from 27 to 41 m, Table 1. The installations were established between 1963 and 1970, on sites with elevations ranging from 274 to 823 m above sea level. Measurements from five of the six stands were used to develop models for predicting wood density (known hereafter as the in-sample dataset). Following model development, and candidate model selection, measurements from the sixth stand were obtained to test model portability and predictive ability with new observations (the external dataset). The external dataset was collected one year after the in-sample dataset. The five stands of the in-sample dataset were part of the Levels-Of-Growing-Stock (L.O.G.S) Cooperative Study in Douglas-fir (Williamson & Staebler 1971; Marshall et al. 1992). Each of the L.O.G.S stands comprised plots with three differing treatments; a control (i.e. no treatment and plot density greater than 2450 stems per hectare), a light thinning treatment (for which 70% of basal area was retained), and a heavy thinning treatment (for which 30% of basal area was retained). Plot densities were maintained over time, relative to the controls, through repeated thinning treatments. The sixth stand, though adjacent to one of the L.O.G.S sites, was in no way related to L.O.G.S. This completely independent stand was managed differently. It comprised four combinations of thinning and fertilisation treatments; 1) no thinning or fertilisation; 2) no thinning, but fertilisation with 448 kg N ha-1; 3) thinning with 1/3 of the basal area retained, but no fertilisation; and 4) thinning with 1/3 of the basal area retained and fertilisation with 448 kg N.ha-1. The sixth stand, with a different study design and evaluated with a different instrument, using different settings, provided an external dataset (i.e. completely separate data) and an opportunity to evaluate model portability and accuracy. Tree measurements Trees were randomly chosen from plots by stratification into three diameter at breast height (DBH) classes which differed with site and treatment. The three classes approximately represented boxplot statistics (i.e. minimum DBH to lower quartile, lower quartile to upper quartile, and upper quartile to maximum DBH) thus ensuring that DBH distributions were approximately normal. In total, 172 trees were sampled with 133 trees used for model development (five locations x three treatments x three replicates x three trees = 135, minus two trees for which data were missing due to issues with sample preparation and collection of field micro-drill data) and 37 trees for evaluation of model portability (one location x four treatments x two replicates x five trees = 40, minus three trees with unreliable micro- drill data). The two datasets were collected during late summer / early autumn at the end of growing season. There was a one-year difference between collection of the two datasets, however weather conditions were very similar (dry and warm) and typical for the season in the Pacific Northwest. Each tree was measured at breast height for diameter, by micro-drill resistance methods for amplitude, and by time-of-flight tools for acoustic velocity. A 5-mm core sample, also taken at breast height, was extracted from each tree for wood density assessment using X-ray Todoroki et al. New Zealand Journal of Forestry Science (2021) 51:6 Page 3 TABLE 1: Description of the study sites Site Latitude Longitude Elevation Origin Site Index Established* Age in study (m, a.s.l.) (m, age 50) (years) Hoskins 44°41’ 123°30’ 305 Natural 41 1963 66 Iron Creek 46°24’ 121°59’ 762 Planted 40 1966 60 Sayward 50°04’ 125°35’ 274 Planted 34 1969 62 Shawnigan 48°38’ 123°43’ 335 Planted 29 1970 64 Stampede 42°53’ 122°49’ 823 Natural 34 1968 75 Shawnigan 48°38’ 123°43’ 335 Planted 27 1971 64 *Year of establishment indicates when experimental sites were established, not planted. TABLE 1: Location and description of study sites comprising in-sample and external datasets. densitometry. This approach was used to provide an accurate measure of wood quality, and as a benchmark against which the non-destructive measures could be evaluated. Cores were initially frozen to prevent mold and stain, then dried at 50℃ for 24 hours. Following conditioning (for at least 48 hours) to attain a uniform moisture content of 8%, X-ray densitometry profiles were generated at a resolution of 0.06 mm. Profiles included ring width, earlywood and latewood densities, and earlywood/latewood proportions. Ring density was calculated using the weighted areas of earlywood and latewood. Mean wood density of each tree was calculated using the weighted average of ring density for the length of the core from inside bark to pith. Micro-drill resistance methods were applied from bark to pith in close proximity to the location of the core samples. In-sample dataset trees were drilled using an IML RESI F400-S tool at the maximum constant forward speed of 150 cm/min, while external dataset trees were drilled using a newer Resi PD500 at a constant forward speed of 25 cm/min. The rotational frequency used for the latter tool was 1500 rpm while that for the former, while not recorded, was known to be at a value between 400-1200 rpm. Together, the selection of rotational frequency and forward speed settings is important for preventing overloading the motor, while drilling speed selection can influence drilling resistance measurements and subsequent prediction of wood density (Sharapov et al. 2019a). In general, resistance amplitudes tend to be lower at lower forward speeds, but the difference may be equal to several orders of magnitude (Mattheck et al. 1997, Rinn 2015). Both instruments converted variations in torque into graphical and digital outputs of path length of the drill (measured at a resolution of 0.1 mm) and the relative resistance, given as an amplitude percentage, that the drilling bit encountered. Before average amplitude was calculated, the initial portion of the profile through the bark was removed as was any portion that extended beyond the pith; identified by a bowl-like shape (Rinn 2012; Fundova et al. 2018). An average of the amplitude profile for each tree was calculated using six approaches. The first approach, the simplest, determined the arithmetic mean of the profile. Since mean values are sensitive to outliers, with further bias introduced through increasing drill resistance, the second approach overcame these shortcomings and detrended the profiles through construction of trendlines that smoothed over fluctuations, and through correction of the baseline. The third approach was similar to the second, but after adjusting/detrending the profile, a value determined within the first 10 cm of the inside-bark profile was added to the mean. The third approach was motivated by the study of Bouffier et al. (2008). However, rather than using 5 cm, as in their approach, we chose 10 cm due to our trees being considerably older, and hence larger. The remaining three approaches scaled the mean amplitudes to their respective z-scores to enable comparison of values from different samples (which may have different means and standard deviations, as can be the case when drills are Todoroki et al. New Zealand Journal of Forestry Science (2021) 51:6 Page 4 operated at different speeds, or different brands of drills are used). The mean amplitude values, X, were scaled using the usual method of calculation (Draper & Smith 1966) i.e. subtraction of the mean, μ, and dividing the result by the standard deviation, σ, i.e. Z = (X – μ) / σ. Details of the six approaches follow: A0: The arithmetic mean amplitude of unadjusted data. A1:The average difference between centered moving means and centered moving minimums (window widths of 10 and 100 respectively), plus the average centered moving minimum within 100 mm of the inside-bark profile (following Isik & Li 2003). A2: The average difference between centered moving means and centered moving minimums (window widths of 10 and 100 respectively), plus the average amplitude within 100 mm of the inside-bark profile. Z0: The standardised equivalent of A0. Z1: The standardised equivalent of A1. Z2: The standardised equivalent of A2. With A0 being the unadjusted mean, A1 will be less than or equal to A0, and A2 greater than A0 (Figure 1). Data distributions for A0, A1, and A2 are preserved after standardisation to Z0, Z1, and Z2 (Figure 2), and therefore performance metrics arising from linear models with either standardised (Z0, Z1, Z2) or non- standardised data (A0, A1, A2) will be identical for in- sample datasets. Trees were also acoustically assessed for time- of-flight using a Hitman ST300 (Paradis et al. 2013). Sensors were centered at breast height and placed approximately one metre apart. Measurements were taken on opposite sides and oriented perpendicular to the slope orientation to avoid reaction wood. Four measurements were taken per tree, and the mean of the measurements for each tree (automatically calculated by the Hitman) used in the analyses. A summary of all data collected, for both in-sample and external datasets, is given in Table 2. Analyses To account for experimental variability and structure inherent in our data, linear mixed-effects models that allow for random group effects, were developed to estimate wood density. The models were developed using data from the 133 in-sample trees, and were formulated as: y = Xβ + Zu + ε where y is the response vector (mean wood density obtained by X-ray densitometry for each of the 133 trees), X and Z are matrices of explanatory variables corresponding to fixed (observed/measured variables) Todoroki et al. New Zealand Journal of Forestry Science (2021) 51:6 Page 5 FIGURE 1: Example of amplitude profiles using an IML RES F400-S (a) and a Resistograph PD500 (b), before (red) and after (blue) adjustment. Mean amplitudes for the unadjusted profile, A0, and adjusted variants, A1, A2, are indicated by horizontal solid red, and dashed blue and green lines respectively. FIGURE 2: Kernel density profiles of non-standardised mean resistograph amplitudes (A0, A1, A2), and standardised equivalents (Z0, Z1, Z2) for in-sample and external datasets. Sample size is indicated by N, and bandwidth by bw. and random effects (unobserved/unmeasured variables) respectively, β and u are the corresponding vectors of parameters for the respective fixed and random effects, and ε is a vector of random errors. Explanatory variables included acoustic velocity, V (km.s-1), breast-height diameter, DBH (cm), stand age (years), site index, SI (m, age 50), and micro-drill amplitude (% for A0, A1, A2, unitless otherwise). Each of the six amplitude variants was sequentially evaluated. The quotient of V and Age, V/Age, was also evaluated as an explanatory variable. Inverse and natural logarithmic transformations of variables were examined and selected by inspecting plots for linearity and constant variance with increasing mean values of the dependent variable. Models developed with combinations of explanatory variables, including single variable models, were explored. For all models, site and plot (plot nested within site) were formulated as random effects. The models were fitted using the Restricted Maximum Likelihood method (REML, Searle et al. 1992) and developed using the linear and nonlinear mixed-effects models package, “nlme”, (Pinheiro et al. 2019) within the R environment (R Core Team 2018). The significance of explanatory variables was evaluated with α = 0.05 using conditional t-tests and F-tests (Pinheiro & Bates 2006). Moving (rolling, running) means were facilitated through the “caTools” package (Tuszynski 2019) and data frame manipulations through the “plyr” package (Wickham 2011). Performance of models that satisfied the conditional t and F tests, were evaluated using: AIC (Akaike 1974), R2marg and R 2 cond following Johnson (2014), mean absolute error, MAE (Equation 1), mean absolute percentage error, MAPE (Equation 2), and root-mean-squared error, RMSE (Equation 3). Multiple metrics were applied because the use of a single metric could lead to an incorrect interpretation. Though we report both R2marg and R 2 cond, our focus is primarily on the former metric, similar to the coefficient of determination, R2. Diagnostics of the models included plots of residuals against fitted values, and plots of observed values versus fitted values. All predictions and subsequent calculations of MAE, MAPE, and RMSE were made at the population level, because, in practice, contributions due to random effects are unknown. (1) (2) (3) where yi and ŷi are observed and predicted values, and N the sample size. Candidate models were selected using all but the R2cond performance metric. The models were first grouped by number of variables (as this influences performance metrics). Within each model subgroup, any model having at least one performance metric in the top two of that metric was selected as a candidate model. This process essentially filtered out poorer-performing models. Candidate models were then evaluated using external data and the two best candidate models selected based on MAPE, MAE, RMSE metrics and their model parameters presented. The simpler of the two models (since simple models are generally better in practice) was then selected for further analyses with overall performance demonstrated by the percentages of all trees within 25 kg.m-3, 50 kg.m-3, 5% and 10% of their true values. Summary statistics are reported as means with variability indicated by standard errors (SE). Todoroki et al. New Zealand Journal of Forestry Science (2021) 51:6 Page 6 Site n Density DBH V A0 A1 A2 Z0 Z1 Z2 (kg.m-3) (cm) (km.s-1) (%) (%) (%) In-sample dataset Hoskins 25 465±7 46±2 4.4±0.1 45±3 40±3 54±3 0.5±0.3 0.3±0.3 0.4±0.3 Iron Creek 25 450±7 42±2 4.5±0.1 35±2 32±2 45±2 -0.5±0.2 -0.5±0.2 -0.5±0.2 Sayward 27 500±5 33±2 4.7±0.1 40±2 39±2 53±2 0.1±0.2 0.2±0.2 0.3±0.2 Shawnigan 27 522±6 27±2 4.6±0.1 40±2 39±2 50±2 0.0±0.2 0.2±0.2 0.0±0.2 Stampede 27 474±6 39±2 4.6±0.1 38±1 35±1 47±1 -0.1±0.1 -0.2±0.1 -0.3±0.1 All sites 133 483±4 37±1 4.6±0.1 40±1 37±1 50±1 0.0±0.1 0.0±0.1 0.0±0.1 External dataset Shawnigan 37 590±6 24±1 4.6±0.1 6±0 6±0 8±0 0.0±0.2 0.0±0.2 0.0±0.2 TABLE 2: Summary data (mean ± se) of the in-sample and external datasets. DBH = diameter; V = acoustic velocity; A0 = mean amplitude; A1 = mean adjusted amplitude shifted by the mean rolling minimum amplitude within 100 mm of bark; A2 = mean adjusted amplitude shifted by the mean amplitude within 100 mm of bark; Z0, Z1, Z2 are standardised equivalents of A0, A1, A2 respectively. Results Explanatory variables All variants of mean resistograph amplitude were highly significant as explanatory variables for predicting wood density, as was site index (transformed), acoustic velocity, V, and the quotient V/age. However, age was not significant as an explanatory variable. Diameter was significant as a predictor of wood density only when in combination with V (Table 3). Model performance with in-sample data and candidate models The optimal model, i.e. the model with the best performance metrics (Table 3), comprised three explanatory variables: 1/SI, V/Age, and A2 (or Z2). For this model AIC, MAE, MAPE, and RMSE were the lowest amongst all models, and R2marg (0.60) the highest. R2cond was 0.61, which indicated that the fixed effects accounted for the majority of the variance. MAE was 19.4 kg.m-3, MAPE 4.1%, and RMSE 25.0 kg.m-3. Another three-variable model, with fixed effects X = [1/SI, V, A2] (or similarly X = [1/SI, V, Z2]), also recorded very good performance statistics: R2marg 0.58, MAE 19.9 kg.m -3, MAPE 4.2%, and RMSE 25.4 kg.m-3. These two models were selected as candidate models for further evaluation with the external dataset. R2marg for all remaining three- variable models ranged 0.39 to 0.55, MAE from 12.0 to 24.6 kg.m-3, MAPE from 4.5 to 5.2%, and RMSE from 26.5 to 30.8 kg.m-3 (Table 3). The poorest performance metrics were all attributed to one model with X = [1/ SI, V, DBH]; the only model within this group without amplitude as an explanatory variable. Amongst two-variable models, all models with transformed site index and amplitude (all variants) were selected as candidate models, as they all had at least one Todoroki et al. New Zealand Journal of Forestry Science (2021) 51:6 Page 7 Explanatory variables 1, X R2marg R 2 cond AIC MAE (kg.m-3) MAPE (%) RMSE (kg.m-3) 1/SI 0.33 0.42 1307 26.0 5.5 32.2 V 0.03 0.48 1309 30.6 6.5 38.7 V/Age 0.05 0.50 1308 30.6 6.5 38.9 A0 (Z0) 0.12 0.56 1288 30.2 6.3 37.2 A1 (Z1) 0.14 0.53 1286 28.6 6.0 35.6 A2 (Z2) 0.18 0.59 1272 28.0 5.9 35.1 V, DBH 0.06 0.57 1304 33.0 7.0 41.0 1/SI, V 0.36 0.46 1302 25.2 5.4 31.4 V/Age, A0 (Z0) 0.17 0.59 1277 28.9 6.1 35.7 V/Age, A1 (Z1) 0.18 0.56 1279 27.8 5.9 34.6 V/Age, A2 (Z2) 0.23 0.63 1264 27.2 5.7 34.1 1/SI, V/Age 0.38 0.46 1300 24.7 5.3 31.0 1/SI, A0 (Z0) 0.48 0.52 1278 23.2 4.9 28.0 1/SI, A1 (Z1) 0.49 0.52 1275 22.6 4.8 28.0 1/SI, A2 (Z2) 0.56 0.57 1258 20.9 4.4 26.2 1/SI, V, DBH 0.39 0.49 1297 24.6 5.2 30.8 1/SI, V, A0 (Z0) 0.52 0.56 1268 21.8 4.6 27.1 1/SI, V, A1 (Z1) 0.52 0.56 1269 21.8 4.6 27.2 1/SI, V, A2 (Z2) 0.58 0.61 1251 19.9 4.2 25.4 1/SI, V/Age, A0 (Z0) 0.55 0.56 1263 21.0 4.4 0.35 1/SI. V/Age, A1 (Z1) 0.54 0.56 1265 21.1 4.5 26.7 1/SI, V/Age, A2 (Z2) 0.60 0.61 1248 19.4 4.1 25.0 1 Only models for which all explanatory variables, X, are significant are shown. SI = Site Index; V = acoustic velocity; A0 = mean amplitude; A1 = mean adjusted amplitude shifted by the mean rolling minimum amplitude within 100 mm of bark; A2 = mean adjusted amplitude shifted by the mean amplitude within 100 mm of bark; Z0, Z1, Z2 are standardised equivalents of A0, A1, A2 respectively (and can replace the non-standardised values with no change in fit statistics); DBH = diameter at breast height. TABLE 3: Performance metrics of mixed-effects models for estimating wood density, grouped by the number of explanatory variables. The two best metrics, used for candidate model selection, within each group, are shown in bold font. performance metric in the top two of that metric. For these candidate models R2marg ranged from 0.48 to 0.56, MAE from 20.9 to 23.2 kg.m-3, MAPE from 4.4 to 4.9%, and RMSE from 26.2 to 28.0 kg.m-3. One further candidate model with fixed effects X = [V/Age, A2], had the second- best AIC, but other statistics were comparatively low (R2marg = 0.23). Amongst single-variable models, transformed site index had the best R2marg (0.33), MAE (26.0 kg.m -3), MAPE (5.5%) and RMSE (32.2 kg.m-3). However, the best AIC was associated with the model with amplitude (A2, or equivalently Z2) as the sole predictor of wood density. Though also having the second-best R2marg, MAE and MAPE, these statistics, especially R2marg, were somewhat lower (0.18, 28.0 kg.m-3 and 5.9% respectively). Candidate model performance with external data Performance metrics of all candidate models evaluated with the external dataset are shown in Table 4. Metrics for models with non-standardised mean amplitudes (A0, A1, A2) are provided for comparative purposes only. Standardised micro-drill amplitudes recorded substantially better performance statistics than their non-standardised counterparts. This was expected since the standardisation procedure converts values of the two tools to the same scale (whist retaining distribution properties). The comparison clearly demonstrates the need for standardisation. Hereafter, the focus is on the standardised amplitude variants. The optimal model, X = [1/SI, V/Age, Z2], had the best performance metrics amongst the three-variable candidate models, however was second best in terms of the performance metrics MAE, MAPE, and RMSE. The model with the best performance metrics comprised 1/ SI and Z0, the unadjusted amplitude. Based on evaluation with the external dataset, the two candidate models selected were: • The optimal model in the in-sample dataset, with X = [1/SI, V/Age, Z2] and • The simpler 2-variable model, with X = [1/SI, Z0]. Parameters of the two selected candidate models are provided in Table 5. Random effects due to site were negligible (i.e. close to zero) for the optimal model, X = [1/SI, V/Age, Z2], and ranged from -8.9 to 7.7 kg.m-3 for the simpler model, X = [1/SI, Z0]. In contrast, random effects due to plot ranged from -1.2 to 1.4 kg.m-3 for the optimal model and were negligible for the simpler model. Residual plots of the two models shown in Figure 3 appear to be reasonably random. Though there is a small degree of asymmetry in both residual plots, overall there are no clear patterns. The simpler model with just two predictor variables X = [1/SI, Z0], though having an outlier with a residual of -82 kg.m-3, does not appear Todoroki et al. New Zealand Journal of Forestry Science (2021) 51:6 Page 8 Explanatory variables, X MAE (kg.m-3) MAPE (%) RMSE (kg.m-3) 1/SI 53.7 8.8 64.1 A2 182.0 30.6 185.0 Z2 107.8 18.1 111.0 V/Age, A2 179.4 30.2 183.0 V/Age, Z2 55.1 9.1 110.0 1/SI, A0 100.8 16.8 107.0 1/SI, Z0 51.6 8.5 59.3 1/SI, A1 112.5 18.8 117.8 1/SI, Z1 56.4 9.3 63.4 1/SI, A2 130.4 21.8 135.0 1/SI, Z2 53.3 8.8 60.2 1/SI, V, A2 131.5 22.0 136.0 1/SI, V, Z2 55.0 9.1 61.9 1/SI, V/Age, A2 130.3 21.8 135.0 1/SI, V/Age, Z2 52.4 8.7 59.5 TABLE 4: Performance metrics for the external dataset, grouped by the number of explanatory variables in the model. The two best metrics within each group are shown in bold. Grey font indicates the non-standardised amplitudes; shown for comparative purpose only. SI = Site Index; V = acoustic velocity; A0 = mean amplitude; A1 = mean adjusted amplitude shifted by the mean rolling minimum amplitude within 100 mm of bark; A2 = mean adjusted amplitude shifted by the mean amplitude within 100 mm of bark; Z0, Z1, Z2 are standardised equivalents of A0, A1, A2 respectively. Model description Explanatory variable Value SE DF t-value p-value Optimal Intercept 240.8 24.1 86 10.0 0.0000 Z2 18.5 2.2 86 8.3 0.0000 1/SI 6346.0 613.8 3 10.3 0.0019 V/Age 862.3 238.0 86 3.6 0.0005 Best 2-var candidate Intercept 294.3 35.0 87 8.4 0.0000 Z0 14.8 2.5 87 5.8 0.0000 1/SI 6587.9 1216.6 3 5.4 0.0124 TABLE 5: Parameters and statistics of two candidate models for predicting wood density in standing trees. Z0, Z2 = standardised mean resistograph amplitudes (unadjusted and adjusted respectively); SI = site index; V = acoustic velocity. to be markedly worse than the more complex model. Therefore, for operational purposes, the simpler model was selected as the model of choice. Confidence and prediction intervals determined using the model of choice (with X = [1/SI, Z0] and parameters as in Table 5) are demonstrated in Figure 4 for each of the sites. Both Z0 and A0 are indicated on the figures, along with the measured densitometric values. The influence of site index is clear. In general, the model of choice estimated wood density with a good degree of accuracy (Figure 5). Nearly all trees in the in-sample dataset (93%) were within 50 kg.m-3 of their true values, and well over half (62%) were within 25 kg.m-3 of their true values. For external data, extrapolation of the model beyond 570 kg.m-3 resulted in wood density being under-estimated (Figure 5). Estimates beyond 600 kg.m-3 were amongst the worse, particularly for those trees on plots which had been fertilised and/or thinned. All external dataset trees on control plots which had neither been thinned nor fertilised were within 54 kg.m- 3 of their true values. For both in-sample and external datasets, 84% of all predictions were within 50 kg.m-3 of their true values, and 53% within 25 kg.m-3 of their true values. Overall, 87% of all predictions were within 10% of their true values, and more than half (54%) within 5% of their true values. Discussion The ability to non-destructively evaluate and predict wood quality is of great importance and has been reported by many authors for a variety of tools (Cown 1978; Wessels et al. 2011) and a variety of purposes such as assessing young trees for genetic heritability (Gantz 2002; Fundova et al. 2019), determining effects of silvicultural practices on product quality (Wang 1999; Briggs et al. 2008), and evaluating wood composites (Winistorfer et al. 1995). Many studies have done so in stands that are more homogenous in terms of age, geographic location, and genetic composition (Desponts et al. 2017). Because many tools and techniques are labour intensive, not field-based, or are too destructive, the micro-drill for non-destructively determining wood quality attributes has attracted attention. It meets the requirements of being portable, inexpensive to use, and has little impact on tested trees. So, can micro-drill resistance tools alone provide accurate assessments of wood density for a diverse set of trees? In our study, mean resistograph amplitude was highly significant as an explanatory variable for predicting wood density, however, the correlation between wood density and amplitude was weak (0.42, based on R2marg = 0.18). This correlation was a little lower than found in the literature for less diverse cohorts; e.g. R2 = 0.22, based on a correlation of 0.47 for 32-year-old Douglas-fir trees from four comparable sites (El-Kassaby et al. 2011), a range in R2 from 0.21 to 0.44 for 11 year- old loblolly pines (Isik & Li 2003), and R2 = 0.38, based on a correlation of 0.62, for a single stand of 25-year- old Douglas-fir trees (Chantre & Rozenberg 1997). Our performance statistics with in-sample data (Table 3) indicate that estimates of wood density of a tree from within these stands (i.e. within the stands of the in- Todoroki et al. New Zealand Journal of Forestry Science (2021) 51:6 Page 9 FIGURE 3: Residual plots of selected models for predicting wood density, with explanatory variables as indicated. Z0, Z2 = standardised mean resistograph amplitudes (unadjusted and adjusted respectively); SI = site index; V = acoustic velocity. Todoroki et al. New Zealand Journal of Forestry Science (2021) 51:6 Page 10 sample dataset) will, on average, have an error of about 28 kg.m-3, or equivalently 6% (Table 3). However, there will be a lot of scatter around this average, some trees will have greater error, others less. If trees are selected from other stands, i.e. new data, then larger errors could be expected. Statistics for our external dataset (Table 4; models shown in blank font) indicated an average error of about 61 kg.m-3, or equivalently 10%. This is quite large, and possibly too large for assessing wood density in individual trees. Improvements in accuracy were obtained when site index was included in the predictive models. This is consistent with findings of Walker et al. (2019). The best accuracy was achieved with a three-variable model, the optimal model, that in addition to adjusted amplitude, included site index and the quotient of acoustic velocity and age. R2marg was 0.60, equivalent to a correlation of 0.77. With this model we could expect accuracy to be within 20 kg.m-3 (MAE was 19.4 kg.m-3), or equivalently 4.1%. RMSE was 25.0 kg.m-3, indicating good proximity between estimated and true means. In comparison to the model in amplitude alone, statistics for our external dataset were greatly improved and approximately halved (MAE 52.4 kg.m-3, MAPE 8.7%, RMSE 59.5 kg.m-3). As noted in the literature, an optimal model for in- sample data may have a lower predictive ability for new/external data than a sub-optimal model (e.g. Shmueli 2010). This was found to be true in our study, and a simpler model, recorded the best accuracy for the external dataset (MAE 51.6 kg.m-3, MAPE 8.5%, RMSE 59.3 kg.m-3). These statistics, while only marginally better than those for the optimal model, are based on a very simple model, X = [1/SI, Z0], requiring input of just site index and mean unadjusted standardised amplitude. We observe that differences in performance metrics within each A0, A1, A2 model triplet (Table 3) are quite small; similarly, for the two-variable group in Table 4. Thus, with the simpler model we forgo the need to adjust the amplitude signal prior to determining the mean. Parameters of the optimal and simpler models indicate that, as site index increases, wood density decreases (due to the inverse relationship, 1/SI), and as amplitude increases wood density increases. Wood density also increases as the ratio of velocity to age increases. Therefore, within a given stand of trees of the same age, wood density increases as acoustic velocity increases. Parameters of the simpler 2-variable model, X = [1/SI, Z0], indicate that for each unit increase in site index, a decrease in wood density of about 6 kg.m-3 could be expected, and for each 0.1 increase in the z-value associated with the standardised (unadjusted) micro- drill amplitude, an increase in wood density of about 1.5 kg.m-3 could be expected. The simple procedure of standardisation provided large improvements in accuracy of estimates with our external dataset. In the scenario that another Douglas- fir dataset becomes available in the future, then all that is required is the computation of mean amplitude for each FIGURE 4: Confidence and prediction intervals of wood density by site index, augmented with data used in this study. Mean amplitude of the inside-bark resistance profiles are shown with both unadjusted (A0) and standardised (Z0) scales. Todoroki et al. New Zealand Journal of Forestry Science (2021) 51:6 Page 11 tree, and thereupon the mean and standard deviation of amplitudes to calculate the z-score. This data, together with site index, could then be evaluated as a further external dataset for the simpler model with X = [1/SI, Z0]. Therefore, though we anticipate this model to be portable, more data from a greater range of site indices would be required, as are a greater range in tree ages to further test model portability. For all predictive models, evaluating performance with other datasets is a crucial step in gauging practical usefulness, yet despite its importance, the primary focus in predictive modelling studies has centred on how models were developed and on their explanatory accuracy. We suggest a greater focus be given to their predictive accuracy through external verification to assess the usefulness of the model in practice. With all tools, including micro-drills, many factors can cause measurement error. When using micro-drills, care needs to be taken to obtain accurate profiles. The angle of penetration of the bit determines the amplitude profile (Rinn 2012) and due to the bit’s flexibility, can be influenced by operator movement, thereby altering the profile (Ukrainetz & O’Neill 2010). Profiles in turn may be affected by moisture content (Isik & Li 2003; Lin et al. 2003) which increases drilling resistance (Kahl et al. 2009). However, in the case of small conditioned wood specimens, the effect of moisture content above fibre saturation (~32% MC) was not evident (Sharapov et al. 2019b). Air temperature (Ukrainetz & O’Neill 2010) and wood properties such as reaction wood, resin pockets, and branches/knots, can also affect the profiles (Eckard et al. 2010). As forest management objectives shift, there is an increased need to better understand the resource for the products/services it can provide. The focus for industrial landowners has been primarily on volume yield for revenue generation. Wood quality, which has a negative relationship with growth (Jozsa & Middleton 1994; Kennedy 1995), is not often taken into consideration. While it can be challenging to predict future customer needs, land managers interested in marketing trees to mills manufacturing products that require a certain wood density level, the ability to plant, grow, manage, and harvest the trees at the optimal time will be of economic benefit. Understanding the effect of land management regimes and integrating production of selected attributes along the value chain from raw materials to products promotes best use allocation of the forest resource in the future while providing information on how to grow and tend trees for specific end uses. Conclusions We have provided evidence that wood density estimates of individual Douglas-fir trees derived from micro- drilling alone are insufficient for obtaining accurate FIGURE 5: Comparison of measured wood density with predictions. Dark and light grey bands represent differences of up to 25 kg.m-3 and 50 kg.m-3 respectively. The vertical line indicates the maximum density within the in- sample dataset. wood density estimates from a diverse set of trees. Site effects need to be considered too and with the simple inclusion of site index in models, wood density predictions improve considerably. Another simple procedure, that of standardising mean amplitudes to z-scores, extends portability of models to future datasets that may use different micro-drills or may operate micro- drills at different speeds. External data are critical to determining how well a model will perform in practice. Competing interests The authors declare that they have no competing interests. Acknowledgements We thank the field crew, laboratory technicians, and database managers who collected, processed, and organised the data, without which this study would not have been possible. And we extend our appreciation to the Stand Management Cooperative, School of Environmental and Forest Sciences, University of Washington, Seattle, WA, USA. Author contributions CT and EL conceived of the study and set overall objectives. CT participated in the design of the study, performed the statistical analyses, and drafted the manuscript. EL and CF acquired and interpreted data and contributed to the writing of the manuscript. All authors have read and approved the final manuscript. Funding This research was funded by an international joint venture agreement between the New Zealand Forest Research Institute Limited, trading as Scion, and the USDA, Forest Service, Pacific Northwest Research Station. Additional funding was provided by the Forest Innovation Program, Government of Canada. References Akaike, H. (1974). A New Look at the Statistical Model Identification. IEEE Transactions on Automatic Control, 19(6), 716-723. https://doi.org/10.1109/ TAC.1974.1100705 Alexander, D. L. J., Tropsha, A., & Winkler, D. A. (2015). Beware of R2: Simple, unambiguous assessment of the prediction accuracy of QSAR and QSPR models. Journal of Chemical Information and Modeling, 55, 1316-1322. https://doi.org/10.1021/acs. jcim.5b00206 Bouffier, L., Charlot, C., Raffin, A., Rozenberg, P., & Kremer, A. (2008). Can wood density be efficiently selected at early stage in maritime pine (Pinus pinaster Ait.)? Annals of Forest Science, 65(1), 1-8. https:// doi.org/10.1051/forest:2007078 Briggs, D.G., Thienel, G., Turnblom, E.C., Lowell, E., Dykstra, D., Ross, R.J., Wang, X., & Carter, P. (2008). Influence of thinning on acoustic velocity of Douglas-fir trees in western Washington and western Oregon. Paper presented at the 15th International Symposium on Nondestructive Testing of Wood, Sept. 10-12, 2007, Minnesota. Chantre, G., & Rozenberg, P. (1997). Can drill resistance profiles (Resistograph) lead to within-profile and within-ring density parameters in Douglas fir wood. Paper presented at the Timber management toward wood quality and endproduct values. CTIA/ IUFRO International Wood Quality Workshop. Chauhan, S.S., & Walker, J.C.F. (2006). Variations in acoustic velocity and density with age, and their interrelationships in radiata pine. Forest Ecology and Management, 229(1-3), 388-394. https://doi. org/10.1016/j.foreco.2006.04.019 Cown, D. (1978). Comparison of the Pilodyn and torsiometer methods for the rapid assessment of wood density in living trees. New Zealand Journal of Forestry Science, 8(3), 384-391. Dankers, F.J.W.M., Traverso, A., Wee, L., & van Kuijk, S.M.J. (2019). Prediction modeling methodology. In P. Kubben, M. Dumontier, A. Dekker (Eds.), Fundamentals of Clinical Data Science (pp. 101- 120). Springer International Publishing. https:// doi.org/10.1007/978-3-319-99713-1_8 Desponts, M., Perron, M., & DeBlois, J. (2017). Rapid assessment of wood traits for large-scale breeding selection in Picea mariana [Mill.] B.S.P. Annals of Forest Science, 74(3): 53. https://doi.org/10.1007/ s13595-017-0646-x Draper, N., & Smith, H. (1966). Applied linear regression: New York: Wiley. Eberhardt, T.L., & Samuelson, L.J. (2015). Collection of wood quality data by X-ray densitometry: a case study with three southern pines. Wood Science and Technology, 49(4), 739-753. https://doi. org/10.1007/s00226-015-0732-x Eckard, J.T., Isik, F., Bullock, B., Li, B., & Gumpertz, M. (2010). Selection efficiency for solid wood traits in Pinus taeda using time-of-flight acoustic and micro-drill resistance methods. Forest Science, 56(3), 233-241. El-Kassaby, Y.A., Mansfield, S., Isik, F., & Stoehr, M. (2011). In situ wood quality assessment in Douglas-fir. Tree Genetics and Genomes, 7(3), 553-561. https:// doi.org/10.1007/s11295-010-0355-1 Filipescu, C.N., Lowell, E.C., Koppenaal, R., & Mitchell, A.K. (2014). Modeling regional and climatic variation of wood density and ring width in intensively managed Douglas-fir. Canadian Journal of Forest Research, 44(3), 220-229. https://doi. org/10.1139/cjfr-2013-0275 Todoroki et al. New Zealand Journal of Forestry Science (2021) 51:6 Page 12 https://doi.org/10.1109/TAC.1974.1100705 https://doi.org/10.1109/TAC.1974.1100705 https://doi.org/10.1051/forest:2007078 https://doi.org/10.1051/forest:2007078 https://doi.org/10.1016/j.foreco.2006.04.019 https://doi.org/10.1016/j.foreco.2006.04.019 https://doi.org/10.1007/s13595-017-0646-x https://doi.org/10.1007/s13595-017-0646-x https://doi.org/10.1007/s00226-015-0732-x https://doi.org/10.1007/s00226-015-0732-x https://doi.org/10.1007/s11295-010-0355-1 https://doi.org/10.1007/s11295-010-0355-1 https://doi.org/10.1139/cjfr-2013-0275 https://doi.org/10.1139/cjfr-2013-0275 Fundova, I., Funda, T., & Wu, H.X. (2018). Non-destructive wood density assessment of Scots pine (Pinus sylvestris L.) using Resistograph and Pilodyn. PLoS ONE, 13(9). https://doi.org/10.1371/journal. pone.0204518 Fundova, I., Funda, T., & Wu, X.H. (2019). Non-destructive assessment of wood stiffness in Scots pine (Pinus sylvestris L.) and its use in forest tree improvement. Forests, 10(6). https://doi.org/10.3390/ f10060491 Gantz, C.H. (2002). Evaluating the efficiency of the resistograph to estimate genetic parameters for wood density in two softwood and two hardwood species. (MSc), North Carolina State University, Raleigh, NC. Gao, S., Wang, X., Brashaw, B.K., Ross, R.J., & Wang, L. (2012). Rapid assessment of wood density of standing tree with nondestructive methods - A review. Paper presented at the International Conference on Biobase Material Science and Engineering, BMSE 2012. https://doi. org/10.1109/BMSE.2012.6466226 Gao, S., Wang, X., Wiemann, M.C., Brashaw, B.K., Ross, R.J., & Wang, L. (2017). A critical analysis of methods for rapid and nondestructive determination of wood density in standing trees. Annals of Forest Science, 74(2): 27. https://doi.org/10.1007/s13595-017- 0623-4 Gwaze, D., & Stevenson, A. (2008). Genetic variation of wood density and its relationship with drill resistance in shortleaf pine. Southern Journal of Applied Forestry, 32(3), 130-133. https://doi. org/10.1093/sjaf/32.3.130 Howe, G.T., Jayawickrama, K., Cherry, M., Johnson, G., & Wheeler, N.C. (2006). Breeding Douglas-fir Plant Breeding Reviews (Vol. 27, pp. 245-353). https:// doi.org/10.1002/9780470650349.ch6 Isik, F., & Li, B. (2003). Rapid assessment of wood density of live trees using the Resistograph for selection in tree improvement programs. Canadian Journal of Forest Research, 33(12), 2426-2435. https://doi. org/10.1139/x03-176 Johnson, P.C.D. (2014). Extension of Nakagawa & Schielzeth’s R2 GLMM to random slopes models. Methods in Ecology and Evolution, 5(9), 944-946. https://doi.org/10.1111/2041-210X.12225 Jozsa, L., & Middleton, G. (1994). A discussion of wood quality attributes and their practical implications. Vancouver, BC: Forintek Canada Corporation. Accessed 4 October 2019 from: https://pdfs. semanticscholar.org/57b5/7386b0f5029e60812695c7 0ac34eb624069b.pdf Kahl, T., Wirth, C., Mund, M., Böhnisch, G., & Schulze, E.D. (2009). Using drill resistance to quantify the density in coarse woody debris of Norway spruce. European Journal of Forest Research, 128(5), 467- 473. https://doi.org/10.1007/s10342-009-0294-2 Kennedy, R.W. (1995). Coniferous wood quality in the future: concerns and strategies. Wood Science and Technology, 29(5), 321-338. https://doi. org/10.1007/BF00202581 Kimberley, M.O., McKinley, R.B., Cown, D.J., & Moore, J.R. (2017). Modelling the variation in wood density of New Zealand-grown Douglas-fir. New Zealand Journal of Forestry Science, 47: 15. https://doi. org/10.1186/s40490-017-0096-0 Lin, C.J., Wang, S.Y., Lin, F.C., & Chiu, C.M. (2003). Effect of moisture content on the drill resistance value in Taiwania plantation wood. Wood and Fiber Science, 35(2), 234-238. Maclaren, J.P. (2009). Douglas-fir Manual. [FRI Bulletin No. 237], 32 p. Rotorua, New Zealand: Scion. Marshall, D.D., Bell, J.F., & Tappeiner, J.C. (1992). Levels- of-growing-stock cooperative study in Douglas- fir: report no. 10 – The Hoskins Study, 1963-83. Research Paper PNW-RP-448. Portland, OR, USA: USDA Forest Service, Pacific Nortwest Research Station. https://doi.org/10.2737/PNW-RP-448 Mattheck, C., Bethge, K., & Albrecht, W. (1997). How to read the results of resistograph m. Arboricultural Journal, 21(4), 331-346. https://doi.org/10.1080/ 03071375.1997.9747179 Mendenhall, W., & Sincich, T. (2012). Principles of model building. In D. Lynch (Ed.), A Second Course In Statistics: Regression Analysis (pp. 261-325). Upper Saddle River, NJ: Prentice Hall. Nakagawa, S., & Schielzeth, H. (2013). A general and simple method for obtaining R2 from generalized linear mixed-effects models. Methods in Ecology and Evolution, 4(2), 133–142. https://doi. org/10.1111/j.2041-210x.2012.00261.x Newton, P.F. (2017). Acoustic-based non-destructive estimation of wood quality attributes within standing red pine trees. Forests, 8(10): 380. https://doi.org/10.3390/f8100380 Paradis, N., Auty, D., Carter, P., & Achim, A. (2013). Using a standing-tree acoustic tool to identify forest stands for the production of mechanically-graded lumber. Sensors (Switzerland), 13(3), 3394-3408. https:// doi.org/10.3390/s130303394 Pinheiro, J., & Bates, D. (2006). Mixed-effects models in S and S-PLUS: Springer Science & Business Media. Pinheiro, J., Bates, D., DebRoy, S., Sarkar, D., & Team, R.C. (2019). nlme: Linear and nonlinear mixed effects models. R package version 3.1-141. https:// CRAN.R-project.org/package=nlme Todoroki et al. New Zealand Journal of Forestry Science (2021) 51:6 Page 13 https://doi.org/10.1371/journal.pone.0204518 https://doi.org/10.1371/journal.pone.0204518 https://doi.org/10.3390/f10060491 https://doi.org/10.3390/f10060491 https://doi.org/10.1109/BMSE.2012.6466226 https://doi.org/10.1109/BMSE.2012.6466226 https://doi.org/10.1007/s13595-017-0623-4 https://doi.org/10.1007/s13595-017-0623-4 https://doi.org/10.1093/sjaf/32.3.130 https://doi.org/10.1093/sjaf/32.3.130 https://doi.org/10.1002/9780470650349.ch6 https://doi.org/10.1002/9780470650349.ch6 https://doi.org/10.1139/x03-176 https://doi.org/10.1139/x03-176 https://doi.org/10.1111/2041-210X.12225 https://doi.org/10.1007/s10342-009-0294-2 https://doi.org/10.1007/BF00202581 https://doi.org/10.1007/BF00202581 https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf https://pdfs.semanticscholar.org/57b5/7386b0f5029e60812695c70ac34eb624069b.pdf R Core Team. (2018). R: A language and environment for statistical computing. Vienna, Austria. https:// www.R-project.org/ Rinn, F. (1988). A new method for measuring tree- ring density parameters. Physics diploma thesis. Institute for Environmental Physics, Heidelberg University, 85 p. Rinn, F. (2012). Basics of typical resistance-drilling profiles. Western Arborist Winter, 30-36. Rinn, F. (2015). Key to evaluating drilling resistance profiles. Western Arborist Fall, 16-21. Schimleck, L., Dahlen, J., Apiolaza, L. A., Downes, G., Emms, G., Evans, R., Moore, J., Van den Bulcke, J., & Wang, X. (2019). Non-destructive evaluation techniques and what they tell us about wood property variation. Forests, 10(9): 728. https://doi. org/10.3390/f10090728 Searle, S., Casella, G., & McCulloch, C. (1992). Variance components. New York: John Wiley & Sons, Inc. https://doi.org/10.1002/9780470316856 Sharapov, E., Brischke, C., Militz, H., & Toropov, A. (2019a). Impact of drill bit feed rate and rotational frequency on the evaluation of wood properties by drilling resistance measurements. International Wood Products Journal, 10(4), 128-138. https:// doi.org/10.1080/20426445.2019.1688455 Sharapov, E., Brischke, C., Militz, H., & Smirnova, E. (2019b). Prediction of modulus of elasticity in static bending and density of wood at different moisture contents and feed rates by drilling resistance measurements. European Journal of Wood and Wood Products, 77(5), 833-842. Shmueli, G. (2010). To explain or to predict? Statistical Science, 25(3), 289-310. https://doi. org/10.1214/10-STS330 Snee, R. (1977). Validation of Regression Models: Methods and Examples. Technometrics. 19, 415- 428. https://doi.org/10.1080/00401706.1977.10 489581 Spiecker, H., Lindner, M., & Schuler, J. (2019). Douglas-fir – an option for Europe. Joensuu, Finland: European Forest Institute. Accessed 10 October 2019 from: https://www.efi.int/sites/default/files/files/ publication-bank/2019/efi_wsctu9_2019.pdf Tuszynski, J. (2019). caTools: Tools: moving window statistics, GIF, Base64, ROC AUC, etc.. R package version 1.17.1.2. https://CRAN.R-project.org/ package=caTools Ukrainetz, N.K., & O’Neill, G.A. (2010). An analysis of sensitivities contributing measurement error to resistograph values. Canadian Journal of Forest Research, 40(4), 806-811. https://doi. org/10.1139/X10-019 Walker, N.K., & Dodd, R.S. (1988). Calculation of wood density variation from x-ray densitometer data. Wood and Fiber Science, 20(1), 35-43. Walker, T.D., Isik, F., & McKeand, S.E. (2019). Genetic variation in acoustic time of flight and drill resistance of juvenile wood in a large loblolly pine breeding population. Forest Science, 65(4), 469- 482. https://doi.org/10.1093/forsci/fxz002 Wang, X. (1999). Stress wave-based nondestructive evaluation (NDE) methods for wood quality of standing trees. (PhD), Michigan Technological University, Houghton, MI, USA. Wang, X., Ross, R.J., McClellan, M., Barbour, R.J., Erickson, J.R., Forsman, J.W., & McGinnis, G.D. (2001). Nondestructive evaluation of standing trees with a stress wave method. Wood and Fiber Science, 33(4), 522-533. Wessels, C.B., Malan, F.S., & Rypstra, T. (2011). A review of measurement methods used on standing trees for the prediction of some mechanical properties of timber. European Journal of Forest Research, 130(6), 881-893. https://doi.org/10.1007/ s10342-011-0484-6 Wickham, H. (2011). The Split-Apply-Combine strategy for data analysis. Journal of Statistical Software, 40(1), 1-29. https://doi.org/10.18637/jss.v040. i01 Williamson, R.L., & Staebler, G.R. (1971). Levels-of- growing-stock cooperative study on Douglas-fir. Rep. No. 1—description of study and existing study areas. Research Paper PNW-111. Portland, OR, USA: USDA Forest Service, Pacific Nortwest Forest and Range Experiment Station. Winistorfer, P.M., Xu, W., & Wimmer, R. (1995). Application of a drill resistance technique for density profile measurement in wood composite panels. Forest Products Journal, 45(6), 90. Zeidler, A., Borůvka, V., & Schönfelder, O. (2017). Comparison of wood quality of Douglas fir and spruce from afforested agricultural land and permanent forest land in the Czech Republic. Forests, 9(1), 13. https://doi.org/10.3390/ f9010013 Zobel, B.J., & Van Buijtenen, J.P. (2012). Wood Variation: its causes and control: Springer Science & Business Media. Todoroki et al. New Zealand Journal of Forestry Science (2021) 51:6 Page 14