Improved models of harvest-induced bark damage

Abstract

Key message

We provided a precise quantitative analysis of the factors at the origin of bark damage during harvesting operations and developed a model able to predict them accurately. The major factors were the distance of trees to skid trails, the intensity of removals, the harvesting system as well as the interactions between the distance of trees to skid trails with harvesting systems, the average skidding distance, the tree species and tree height.

• Context

During timber harvesting, trees in the remaining stand may suffer bark damage resulting from tree-felling or log manipulation. Although a multitude of case studies and empirical observations provide qualitative and quantitative information with respect to the potential causal factors, the basic quantitative relationship between major factors of influence and the resulting degree of bark damage remains largely unclear.

• Aims

The objective was to provide a precise quantitative analysis of impact factors explaining the occurrence of bark damage during harvesting operations.

• Methods

Three different modelling approaches were tested: boosted regression tree (BRT), a generalised linear mixed effects model (GLMM) and Bayesian Markov chain Monte Carlo generalised linear mixed models (MCMCglmm).

• Results

The major factors with a significant impact on the occurrence of bark damage were the distance of trees to skid trails, the intensity of removals, the harvesting system and the interaction term between the distance of trees to skid trails with harvesting systems, average skidding distance, tree species and tree height.

• Conclusion

The final model includes the relevant major factors impacting on the infliction of bark damage during practical harvesting operations. Furthermore, it discriminates well with respect to the occurrence of bark damage, and it provides managers with a rational and conclusive tool for optimising harvesting operations.

1 Introduction

Since the twentieth century, the increased use of mechanised wood harvesting systems has been associated with the problem of greater levels of harvest-induced bark damage inflicted on the remaining trees of forest stands (Vasiliauskas 2001). Bark injuries can be caused by felling or skidding operations. Injuries caused along the work chain of a specific harvesting system can be summarised as “system damage” (Limbeck-Lilienau and Stampfer 2004). A common indicator used to describe the intensity of such bark injuries is the “degree of damage”, which expresses the rate of damaged trees as the percentage of all trees in the remaining stands.

For instance, in an economic context, a bark injury results in economic damage only if the biological injury is associated with degradation of timber quality or in a reduction in growth or yield. However, bark injuries generally increase the risk of fungus infections and are often associated with considerably deteriorated timber quality followed by losses in yield of marketable timber volume. Bark wounds are reported to deteriorate the timber quality by either the incidence of stain (discoloration) or wood decay (Vasiliauskas 2001). The infection of wounds by fungi may lead to serious wood degradation. Following wood infection, for many tree species, decay may even invade the central portions of the stem and form a typical “heart rot” that expands well above and below the surface of the bark wound (Pawsey and Gladman 1965; Shigo 1966; El Atta and Hayes 1987). Norway spruce (Picea abies) in particular is highly susceptible to wound decays that lead to degradation and serious economic losses (Dimitri and Schumann 1975; Kohnle and Kändler 2007; Metzler et al. 2012).

Since the 1970s, several studies in Germany have examined the impact of mechanised logging operations on harvest-related injuries (Guglhör and Melf 1995; Mahler 1987; Meng 1978; Löffler 1975). Furthermore, the results of the second German National Inventory BWI 2 (reference year 2002) underlined the fact that the degree of harvest-related bark injury in German forests has reached seriously problematic levels with particularly high levels found in Southwest Germany (state of Baden-Wuerttemberg). Here, almost 20 % of the inventoried trees displayed visually discernible bark injuries, most likely inflicted during harvesting operations (Polley and Hennig 2005).

Many case studies describe a variety of factors that influence the injury level in the remaining stands during timber harvesting operations. In general, these studies focus on selected usage scenarios, such as, for instance, certain timber harvesting methods with varying degree of mechanisation and selected pre-skidding machines, certain topographic or stand conditions and harvesting intensities. In a comprehensive literature review and overview, Nill (2011) grouped individual studies by reportedly superior influencing factors. He identified the human factor, for example, as a key factor which comprises training and experience and motivation. Kellog et al. (1986), Knorr and Prien (1988), Vospernik (2004), Limbeck-Lilienau and Stampfer (2004), Stampfer et al. (2002) and Schoettle et al. (1999) show that training is actually an essential basis for careful work and that experience, for instance in operating modern harvesting machines, can have a direct impact on the degree of damage. In some of the abovementioned studies, the additional aspect of work motivation is studied and considered inseparable from other factors. All authors agree on the high impact of workers; however, they state that in case studies, the quantification of the impact of the human factor is rather complicated if not completely impossible.

Another set of single factors that are discussed is the preparation for harvesting operations and forest spacing. Specifically, the spacing of skid trails determines the pre-skidding distance and, simultaneously, a higher density of skid trails increases the number of trees that are vulnerable because of their positions along these skid trails (Limbeck-Lilienau and Stampfer 2004; Mahler 1987; Han and Kellogg 2000; Sauter and Grammel 1996; Schoettle et al. 1997; Froese and Han 2006). The variable “proximity of a tree to the edge of the skid trail” could also be as a significant factor in models for damage assessment (Bobik 2008; Froese and Han 2006; Fajvan et al. 2002).

Another set of factors that influence the amount of damage is associated with the harvesting operation, in particular the chosen timber harvesting methods. Here, the type of pre-skidding plays the most important role. A direct comparison of completely different types of harvesting operations is usually difficult because different methods often require specific topographic and site-specific conditions and can rarely be viewed as substitutes (Nill 2011). In addition, certain types of harvesting operations were developed for different bucking schemes producing long or short logs with corresponding machinery. This directly influences the degree of damage caused by the operations, whereas in principle, long logs cause greater damage because of their more difficult manipulating processes (Boltz 1987; Nichols et al. 1994; Limbeck-Lilienau 2003; Suwala 1997; Woods et al. 2007). The above-quoted diverse national and international studies on influencing factors on the degree of bark damage to remaining trees in harvesting operations discuss many single factors. However, merging these factors in statistically significant models for the quantification of the influence of the main factors on the degree of damage and the interaction between the factors is still lacking.

Therefore, the objective of this study was to provide a precise quantitative analysis of the impact factors in order to support decision making with respect to economic management of harvest-induced bark damage. For this, we used a bark damage prognosis model previously developed from a comprehensive database (Nill 2011; Nakou et al. 2013) as a “precursor” model and attempted to improve the model by applying and comparing the effect of three different, more complex modelling approaches (boosted regression trees, generalised linear mixed effects model and Bayesian Markov chain Monte Carlo generalised linear mixed models).

2 Material and methods

2.1 Harvesting operations: selection and inventory

The large number of harvest operations included in the database was gathered by the Forest Research Institute of Baden-Wuerttemberg in collaboration with ForstBW (state forest service) and its technology bases, which represent a network of forest enterprise units with specifically trained practitioners for wood harvesting issues and which supervised and conducted the majority of the associated inventories. In previous studies, this network of forest technology bases has been found to be a particularly suitable tool for acquiring empirical data in large-scale studies in practical forestry (e.g. Kohnle et al. 2006).

2.1.1 Selection of harvest operations

The selection of the inventoried harvesting operations took place in a two-stage process (Nill 2011). The accounting and business management system of ForstBW provides data on current or completed harvesting operations in the state and community forests in Baden-Wuerttemberg. Based on previous studies about bark damage and based on the knowledge of experts, the first step involved scanning this database regularly for harvesting operations that were judged to be potentially relevant for the study. Relevant aspects were the tree species harvested, stand age, slope, spatial stand distribution and so on. Further aspects included harvesting-specific characteristics such as harvesting volume, technique, etc. It was thus ensured that all operations that were carried out during one harvesting season (winter 2007/2008 and summer 2008) were included in the base population, and therefore potentially any single harvesting operation could have been selected. To limit the expenditure for the survey and to ensure the validity of the results, it seemed appropriate to restrict the selection of harvesting operations to more or less single-layered stands.

If a potentially interesting harvesting operation was pre-selected in the first step, a questionnaire was sent to the responsible forest district manager asking for more detailed information on characteristics such as harvesting technique, skid trail design, timber assortment, etc. Final selection of a harvesting operation for the study was based on this auxiliary information, and the respective operation was assigned to one of ten predefined harvesting type categories. In the first step, 645 harvesting operations were pre-selected. In the final selection, 183 harvesting operations were used for the study.

Harvesting operation types were essentially classified according to the question: “What was harvested with which technique in which way, and for what length of distance was the timber hauled through the stand?” (Table 1):

Table 1 Type and number of recorded harvesting operations

“what”—classification of the main timber product (short logs <7 m long, long logs >7 m long, tree-length logs, full-trees)

“with which technique”—skidding methods (cable, cable with grip tong,¹ cable yarder)

“in which way”—degree of mechanisation (fully and/or partially² mechanised, motor manual)

“for what length of distance”—mean distance to skid trail³ (20, 40, >40 m)

2.1.2 Inventory method

Seeing as the inventory methods for the harvesting operations were intended to be applied by the trained forestry practitioners from the network of technology bases, the development of the inventory design had to include the following three aspects:

• Applicability to a wide variety of stands/harvesting operations

• Achievement of the measurements by applying a robust and simple design and technology within a narrow time frame

• Mirroring, in principle, the guidelines for the sample technique used in the inventories conducted by ForstBW as the basis for the decade-related forest management plans in state and communal forests in Baden-Wuerttemberg (Kemner and Risse 1994)

The design conceived for the purpose of the study of harvesting-induced bark damage was thus based on a systematic grid sample laid out according to the forest inventories and the concept introduced by Götschmann and Strebel (1996). In principle, the distribution of the sample points within the stands was oriented in relation to the distribution of the skid trails. In contrast to the forest inventories, all trees within the circular plots were inventoried. However, the size of the circular plots within the stands was varied to ensure that at least ten trees per sample plot were measured. Furthermore, the number of sample plots was determined by taking into account expert judgement on the area actually harvested and the stand homogeneity (10–20 sample points per stand).

All trees were measured for diameter at breast height (at 1.3 m of height; diameter at breast height (dbh)). Likewise, all trees that had been removed at the harvesting operation were identified as newly harvested stumps and measured for diameter at the remaining stump. To assess dbh of the removed trees, the remaining trees were measured for diameter at stump height as well, and the relation of dbh to diameter at stump height was used to estimate dbh of the removed trees. Heights were measured for three randomly selected trees and used to construct stand height curves. These stand height curves were used to estimate the height of trees not measured for height (remaining trees as well as removed trees).

For all trees included in the sample plots, the distance to the closest skid trail was recorded and the remaining trees after harvesting were examined for newly inflicted bark injuries. For the removed trees, the felling direction and the skidding processes involved were assessed and noted. The inventory of bark damage was based on the definition according to Meng (1978), and the number and locations of damage (root, stem below breast height, stem above breast height) were recorded for the respective tree and the single largest bark damage location was measured in detail (area, length, width).

2.2 Model improvement

From the database described above, Nill (2011) had developed a model for the probability of the occurrence of bark damage based on the data from the harvesting operations. The database comprised only damage along the stem; damage at the roots was not included. Furthermore, for the model, the variable “occurrence of bark damage” was binary-coded (“yes/no”). Consequently, Nill’s models deliver the probability of occurrence of bark damage. This probability was estimated using mixed models (generalised linear mixed models (GLMMs)). The model was designed as a GLMM with random effects at the level of the stand. Key criteria for deciding which of the initially tested potential variables should remain in the final model were relevant statistical parameters, such as statistical significance, Akaike information criterion (AIC),⁴ pseudo-Bayesian information criterion (BIC)⁴ and area under the curve (AUC)⁴. The selection processes of the variables included in the model and the modelling method are described in detail in Nill (2011). Later on, Nill’s model was evaluated and subsequently slightly modified (Nakou et al. 2013).

In this study, we attempted to improve this model by not only examining the factors that cause damage to bark, as in the original model of Nill, but also by investigating their interactions. In this optimisation process, three different modelling techniques were applied and compared: boosted regression tree (BRT), GLMM and Bayesian Markov chain Monte Carlo generalised linear mixed model (MCMCglmm). All calculations were performed using the statistical software package “R” (R Development Core Team 2005).

2.2.1 BRT

BRT is a combination of two powerful statistical techniques: boosting and regression trees. Boosting is a machine learning technique similar to model averaging where the results of several competing models are merged. However, unlike model averaging, boosting uses a forward, stage-wise procedure where tree models are fitted iteratively to a subset of the training data. Subsets of the training data used at each iteration of the model fit are randomly selected without replacement where the proportion of the training data used is determined by the modeller as the so-called bag fraction parameter. This procedure is known as stochastic gradient boosting and introduces an element of stochasticity into improving model accuracy and reducing overfitting (Elith et al. 2008).

In our study, 50 trees were fitted initially in the normal manner, using recursive binary partitioning of the data. Residuals from the initial fit were then fitted with another set of 50 trees whose residuals were then again fitted with the next set of trees. This procedure was subsequently continued until a specific loss function was minimised, which was verified by n-fold cross-validation. In the case of regression trees, the loss function to minimise is represented by model deviance. The final fit is then based on the entire dataset and computed as the sum of all trees multiplied by the learning rate (Elith et al. 2008).

In fitting a BRT, two parameters must be specified: learning rate and tree complexity. The learning rate represents the contribution of each successive tree to the final model as it proceeds through the iterations. The tree complexity determines whether the model will consist of main effects only (tree complexity = 1) or whether interactions between the factors should be included (tree complexity = 2, 3, …). Ultimately, the combination of learning rate and tree complexity determines the total number of trees in the final model.

In BRTs, the standard output does not include component-wise p values for each of the predictors in the model because term-wise p values are not easily determined. Instead, an index of relative influence may be calculated by summing up the contribution of each variable, which is equivalent to summing up the branch length for each variable in the regression tree.

The fitted BRT models compared in this study were obtained using the BRT script provided by Elith et al. (2008) implemented in R’s “gbm” library (Ridgeway 2007). The parameters applied for fitting the BRT were set as follows: learning rate 0.01, tree complexity 10, cross-validation 10-fold and bag fraction 0.8.

2.2.2 GLMM

GLMMs are regression models that allow for selecting, as with generalised linear models (GLM), from various distributions and link functions in order to model a wide range of types of dependent variables through linear combinations of one or multiple predictor variables (“fixed effects”). Additionally, GLMMs include so-called random effects. These random effects quantify the variation of regression intercept or slopes amongst different grouping levels of the dependent variable (“grouping variable”) by a probability distribution instead of estimating a fixed regression coefficient for each level.

A grouping variable should be used as a random effect if its levels may be conceived as representing a random sample from a larger group (O’Hara 2009). Often, the levels of grouping variables are meaningless but need to be taken into account in order to obtain valid p values and estimates. Therefore, such “nuisance variables” should be used as random effects whenever their variation is of interest. In all other cases, variables should generally only be used as fixed effects (Robinson 1991). Thus, it would be proper to use grouping or nuisance variables as fixed effects. However, random effects have the practical advantage of using fewer degrees of freedom, particularly if random effects comprise numerous levels. In this study, the random effect is at the stand level (183 different harvesting operations).

The first step of the modelling process of GLMM is to identify an appropriate distribution and link function for the data, if necessary. The purpose of the link function is to transform the values of the dependent variable so that they match the scale of the linear predictor (i.e. [−∞, ∞]) and to linearise the relationship with the predictor variables. For each distribution, a canonical (“natural”) link function exists. However, there are also less commonly used alternatives that may suit the data better in some cases. As the response variable in our study was binary, we choose a logit link function.

Before actually calculating the model, we needed to consider which estimation method to use. This choice generally depends on the dependent variable and on the random effects that are to be included in the model (cf. Bolker et al. 2009). For our study, the Laplace approximation method appeared to be best suited. The computation of p values for the significance of fixed effects was done with Wald chi-square tests and for the random effects with likelihood ratio (LR) tests. The computation of p values for the significance of the single parameters was done with Wald Z tests. The model selection criteria were AIC, pseudo-BIC and AUC value as discussed above.

The subsequent calculation of the relative influence of each predictor (fixed effect) was executed by fitting a model for each predictor without the respective predictor. The degree of deterioration of the statistical fit triggered by the omission of the respective predictor served as an indicator of the predictor’s impact. Dividing this value by the sum of all deteriorations aggregated for all predictors rendered a relative proportion of the respective predictor’s impact.

Overdispersion is an important concept in the analysis of discrete data. Overdispersion is the situation that occurs most frequently in Poisson and binomial regression when variance is much higher than the mean (whereas it should ideally be the same). It is evident with a high (>2) residual mean deviance (which should normally be around one) and the presence of too many outliers. The reasons for overdispersion may be outliers, misspecification of the model, variation between the response probabilities or correlation between the binary responses. It results in incorrect estimations of standard error and confidence interval. In this study, the degree of overdispersion from the residual deviance (the residual scaled deviance should be roughly equal to the residual degrees of freedom) was examined.

The GLMMs were fitted using the “glmer” function of the “lme4 package” in R.

2.2.3 MCMCglmm

As outlined above, GLMMs are widely used for regression analysis in different research areas, as they are flexible and efficient for the analysis of grouped data. Several estimation methods have been applied to these models. However, for the non-Gaussian distributed response variables, such as in this study, there is the restriction that the likelihood cannot be obtained in closed form. For instance, McCulloch and Searle (2001) demonstrated that generalising linear mixed models to non-Gaussian distributed data is questionable because integrating over the random effects proves intractable.

This challenge may be solved with MCMCglmm methods by sampling from a series of simpler conditional distributions that can actually be evaluated. This provides an alternative strategy for marginalising the random effects and may be more robust (Zhao et al. 2006; Brown and Draper 2006).

In principle, Markov chain Monte Carlo (MCMC) methods represent a bridge between traditional frequency-based methods and Bayesian methods: from an iterative estimation technique, MCMC builds an empirical distribution of statistical parameter estimates. MCMC provides initial parameter estimates from a priori information that is used to evaluate subsequent iteratively remodelled parameter estimates (e.g. using Bayesian computational methods). This a priori information is called “the prior” and is generally extracted from the data under investigation in the form of a specified distribution or several distributions if multiple parameters are involved. The prior represents a reference point against which iteratively produced parameter estimates are evaluated, thus ensuring convergence on the best set of estimates possible. The empirical distribution of estimates built through the MCMC method is called the posterior distribution because it is created after the data has been fitted.

Developing MCMC methods for GLMM has currently been under intensive investigation (e.g. Zeger and Karim 1991; Damien et al. 1999; Sorensen and Gianola 2002; Zhao et al. 2006), and several software packages are now available which implement such techniques. For the purposes of our study, we used R’s MCMCglmm package (Hadfield 2010). This package implements Markov chain Monte Carlo routines for fitting multi-response GLMM. To test the statistical significance of variance components in MCMCglmm, we used the deviance information criterion (DIC). The implementation of DIC in MCMCglmm is further described in R’s reference manual. By default, DIC values are calculated by MCMCglmm. DIC balances the model and the parameter number simultaneously, and small values of DIC are preferred.

In our study, because the response variable was binary, we used a binomial model with logit link function. Based on the nature of the database at our disposal (large number of data, rather complex structure), we choose an uninformative type of prior. To examine the significance of fixed effects as main effects, we used MCMC tests. We calculated the relative influence for each predictor of our model by fitting a model for each predictor without the respective predictor. The deterioration of the fit statistic (DIC value) due to omitting the respective predictor was used as an indicator of the impact. Dividing this value by the sum of all deteriorations, we obtained the relative proportion.

3 Results

3.1 BRT

The optimal number of regression trees was reached at 1,150 regression trees. The final BRT model accounted for 20 % of the mean total deviance; the value of AUC value was 0.8 indicating that the model discriminates the occurrence of bark damage very well.

For the BRT fitted here, the six most influential variables were (1) intensity of removals (28 %), (2) distance of trees to skid trails (25 %), (3) harvesting system (16 %), (4) mean skidding distance (13 %), (5) tree height (12 %) and (6) tree species (6 %).

The information about the distribution of the fitted values in relation to each predictor, in the case of BRT, was obtained from term-wise plots of fitted functions versus observed values (Fig. 1).

Fig. 1

BRT fitted functions for each term in the BRT main effects model ranked by relative influence value

The first plot shows that the probability of bark damage rises with the increasing intensity of removals (Fig. 1a). The same is true for mean skidding distance as soon as mean skidding distance exceeds 20 m (Fig. 1d). The second most influential variable was the distance of damaged trees from skid trails. The respective plot (Fig. 1b) clearly illustrates that trees close to skid trails (less than 10 m) are more likely to be damaged than trees situated further away from a skid trail. Although there is a trend that trees taller than 25 m are damaged more often than smaller trees, tree height in general contributes little to explaining the likelihood of bark damage occurring (Fig. 1e). The least influential, significant variable was tree species. Norway spruce and beech display the same sensitivity to bark damage, whereas fir, oak, Scots pine and larch are significantly less affected (Fig. 1f).

The BRTs identified the harvesting system (Fig. 1c) as one of the most influential variables for the occurrence of bark damage. Harvesting systems that integrate the use of a cable yarder or a cable skidder (cable winch only) for skidding logs are shown to cause more bark damage than skidders using a grip tong. Likewise, systems harvesting tree-length logs or full-trees as main products appeared in BRTs as more damaging to the remaining stand than systems harvesting mainly short-length logs. Furthermore, fully and/or partially mechanised harvesting operations showed generally relatively low damage levels. From the performance of an interaction analysis, it became clear that the most important interaction was the interaction between “harvesting systems” and “distance of trees to skid trails”.

In most respects, the results of the BRTs are similar to the results of Nill’s model (2011). The only major difference is that in Nill’s model, the harvesting systems were not considered as having a significant influence on the probability of bark damage occurrence. That clearly happened because in BRT analysis, the interactions between the predictors are automatically included, whereas in GLMM, used in Nill’s models, they are not. To quantify these interactions, a function of BRTs was used that creates for each possible pair of predictors a temporary grid of variables representing combinations of values at fixed intervals along each of their ranges. This analysis showed that the most important interaction was the one between “harvesting systems” and “distance of trees to skid trails”.

3.2 GLMM

Fitted fixed effects comprised overall effects of the harvesting system, mean skidding distance, distance of tree to skid trail, intensity of removals, tree species and stand height. Random effects were quantified for the variation of the fixed-effect parameters across the stands. The final “full” model is presented in Table 2.

Table 2 Parameters of the GLMM model for bark damage on the basis of harvesting operations

Judged by the proportion of their relative impact, the seven most influential variables in GLMM were (1) distance of trees to skid trails (33 %), (2) intensity of removals (18 %), (3) harvesting system (14 %), (4) interaction distance of trees to skid trails with harvesting system (10 %), (5) tree species (10 %), (6) mean skidding distance (9 %) and (7) tree height (6 %). Thus, the results of the GLMMs are quite similar to the BRTs. In both modelling approaches, the variables describing the intensity of removals, distance of trees to skid trails, harvesting systems and the interaction term “distance of trees to skid trails with harvesting systems” together have the relatively strongest impact (ca. 75 %) on the occurrence of bark damage. The only difference in these two modelling approaches is the relative influence of the variable “tree species”.

The model fit had an adjusted R ² of 0.24 and an AUC value of 0.745, indicating that the GLMM discriminates well with respect to the occurrence of bark damage (Table 2). Examining overdispersion indicated that there was no overdispersion in the model. The validation of the model disclosed a success index⁵ of 0.75, a specificity of the model of 85 % (percentage of correctly predicted undamaged trees), a sensitivity of 65 % (proportion of correctly predicted damaged trees) and that a total 85 % of the trees were correctly classified (Table 3).

Table 3 Validation of the GLMM model for bark damage on the basis of harvesting operations

Figure 2 depicts the influence of the different predictors on the occurrence of bark damage. This information was obtained by examining the term-wise plots of fitted values versus the predictors. The fit of the model showed that not only the variable “distance of tree to skid trail” was significant as a main effect (being the most influential single variable), but the interaction of this term with two different harvesting systems was highly significant too.

Fig. 2

GLMM fitted functions for each effect in the model

Although trees that are situated close to skid trails usually carry a higher probability for occurrence of bark damage (Fig. 2c), the first plot (Fig. 2a) shows that in particular for harvesting system n o. 6, the opposite effect exists, which explains in this case the significance of the interaction between this harvesting system and the “distance of tree to skid trail”. The other significant interaction is the one found for harvesting system no. 9 and the “distance of tree to skid trail”. Harvesting system no. 9 generally causes less bark damage than other systems, but especially for this system, the probability of bark damage occurring was very high for trees with positions closer than 10 m to skid trails. This effect is clearly illustrated by the interaction plot (Fig. 2b). Figure 2d shows that that the probability of bark damage rises with increasing intensity of removals. The same can be stated for increasing mean skidding distance (Fig. 2e). Although statistically significant, tree height has no important impact on the magnitude of the occurrence of bark damage (Fig. 2f). Likewise, the variable tree species delivers no clear differentiation related to bark damage: spruce and beech appeared about equally vulnerable, whereas fir, oak, pine and larch showed slightly less bark damage (Fig. 2g).

The harvesting system (Fig. 2h) as a main effect (without interaction with other variables) was one of the most influential variables for the occurrence of bark damage. Harvesting systems that include a cable yarder or a cable skidder (cable winch only) for skidding cause significantly more bark damage than systems employing skidders using a grip tong. Likewise, harvesting systems based on tree-length logs or full-trees as the main assortment appeared in GLMMs as clearly more damaging to the residual stand than systems harvesting mainly short-length logs. Furthermore, fully and/or partially mechanised harvesting operations were generally associated with lower damage levels.

3.3 MCMCglmm

The fixed effects in the MCMCglmm approach quantify the overall effects of harvesting systems, mean skidding distance, distance of tree to skid trail, intensity of removals, tree species and stand height. The random effects quantify the variation across the stands of the fixed-effect parameters. After computation of variance components and confidence intervals for the random effects and the use of DIC and AUC value as model selection criteria, we decided on the final “full” model presented in Table 4.

Table 4 Parameters of the MCMCglmm model for bark damage on the basis of harvesting operations

The results of MCMC tests show that all fixed effects incorporated in the full model were statistically significant. The seven most influential variables, as determined by the calculation of the relative influence, were (1) interaction distance of trees to skid trails with harvesting systems (20 %), (2) distance of trees to skid trails (18 %), (3) intensity of removals (15 %), (4) harvesting system (15 %), (5) tree species (14 %), (6) mean skidding distance (12 %) and (7) tree height (6 %).

The model fit has an AUC value of 0.7 (Table 4). Furthermore, the summary plot of the variance component associated with the random effect in stand level (ID) clearly showed that the variance values of the random effects do not meet the zero line, meaning that the random effect observed for the stand level is significant. The validation of the model results in a percentage of correctly predicted undamaged trees of about 85 % and a proportion of correctly predictedAlthough trees that are damaged trees of about 40 %. Altogether, a total of about 85 % of the trees can be correctly classified (Table 5).

Table 5 Validation of the MCMCglmm model for bark damage on the basis of harvesting operations

Figure 3 depicts how the occurrence of bark damage is influenced in the MCMCglmm by the different predictors. This information was gained by examining the term-wise plots of fitted values versus the predictors. The fit of the model shows, like the GLMM, that not only the “distance of tree to skid trail” as a main effect was significant (the most influential variable) but also the interaction of this term with two different harvesting systems, no. 6 and no. 9, was highly significant as well. Further results are very similar to those stated already for the GLMM model. They are visualised in detail in Fig. 3.

Fig. 3

MCMCglmm fitted functions for each effect in the model

The MCMCglmm approach also disclosed the predictor “harvesting system” as one of the most influential variables for the occurrence of bark damage either in interaction with the predictor “distance of tree to skid trail” or as a predictor on its own. Several effects associated with “harvesting systems” become evident:

• Systems including either a cable yarder or cable tractors equipped only with cable winches for skidding cause more bark damage than skidding operations using skidders additionally equipped with a grip tong.

• Systems based on either tree-length logs or full-trees as the main assortments cause significantly more stand damage than those focused on short-length logs.

• The lowest damage levels appear to be associated with fully mechanised harvesting systems and systems harvesting mainly short-length logs employing skidders using a grip tong (Fig. 3).

4 Discussion

The objective of the presented study was to provide a precise quantitative analysis of impact factors relevant to the occurrence of bark damage during harvesting operations. The final goal is to provide forest managers with decision support for an effective and economically rational management with regard to the occurrence of newly inflicted bark damage during harvesting. For this purpose, three methodologically different modelling approaches were applied to develop and test an improved model from an existing database and precursor model. The main focus was to investigate if these different modelling approaches would yield consistent results and if the newly developed model would actually improve model performance compared to the existing precursor model developed by Nill (2011).

In the present study, models were evaluated comparatively, with the newer BRT and MCMCglmm approaches compared to traditional GLMM.

The boosted regression tree approach was applied because of its advantages (Elith et al. 2008):

• Robust parameter estimation, using the stochastic gradient boosting algorithm to minimise variance and bias

• Reduced risk of misspecification, as the model “learns” from the data

• All the benefits of tree-based models in handling model predictors (e.g. automatic fitting of complex interactions and unaffected by multicollinearity, missing predictor values or outliers)

The GLMMs are widely used for regression analysis. They are flexible and efficient for the analysis of grouped data and they are widely applied. However, McCulloch and Searle (2001) showed that generalising linear mixed models to non-Gaussian data may prove difficult because integrating over the random effects is intractable. Therefore, techniques that approximate these integrals have become popular, and the convenience of using the Laplace approximation method for estimating GLMM (Breslow and Clayton 1993) is routinely practiced in most software, such as the lme4-package in R. However, such approximation can lead to estimators that are asymptotically biased towards zero (Breslow and Lin 1995). Markov chain Monte Carlo methods (MCMCglmm) help solve this problem by sampling from a series of simpler conditional distributions that can be evaluated. Thus, the MCMCglmm procedure provides an alternative strategy for marginalising the random effects that may be more robust (Zhao et al. 2006; Brown and Draper 2006) and was, for this reason, applied in this study along with the traditional GLMM approach.

All three modelling approaches showed quite similar results and thus fulfilled the required consistency. The fit of the models under all three approaches showed satisfactory results with AUC values ranging from 0.7 (MCMCglmm) to 0.8 (BRT), which means that all the models can successfully discriminate the occurrence of bark damage. Likewise, the validation of the models showed that GLMM as well as MCMCglmm could correctly classify 85 % of the trees. If one had to choose from the three tested modelling approaches, we would recommend the GLMM approach based on its superior fit (compared to the MCMCglmm approach) and the more robust model predictions (compared to the BRT approach) (Tables 3, 4 and 5). As a consequence, the interpretation of the GLMM model’s outputs is less complicated.

Summarising the assessment of the models allows for an evaluation of the significance of those factors that are most influential for the infliction of bark damage during harvesting operations. In the GLMM and MCMCglmm approaches, the variable “distance of trees to skid trails” was not only significant as a main effect but also as an interaction term with the variable “harvesting systems”. In particular, GLMM and MCMCglmm showed that the interaction “distance of trees to skid trails” with harvesting system no. 6 (partially mechanised, full tree-length, cable and grip tong, >40 m) and no. 9 (partially mechanised, long logs, grip tong, 40 m) was highly significant and proved to be one of the most influential variables. In the results of the BRT analysis, interactions terms could not appear as significant for the principal reasons that in BRT analysis consideration of interactions between predictors is automatically included and cannot be specified as a separate source of influence. Nevertheless, after performing the interactions analysis function, it became clear that in BRT, like in GLMM or MCMCglmm, the most important interaction was the interaction between “harvesting systems” and “distance of trees to skid trails”.

The high significance of the interaction terms between these two harvesting systems (no. 6 and no. 9) and the variable “distance of tree to skid trail” can be explained by the technical characteristics of these harvesting systems:

Harvesting system no. 6 is a cable logging system in which ground hauling is used for full-trees. The felled tree is pulled by a cable winch through the stand and at the stand edge is manipulated with a grip tong onto the skid trail. Consequently, the trees that were located far away from skid trails suffered more in this system from bark damage than the ones closer to skid trails.

Harvesting system no. 9 is used for long logs. For this system, the probability for the occurrence of bark damage was very high for trees with locations closer than 10 m to skid trails. That happened because when the long logs had to be extracted with the harvester, some of them had to be moved into the stand located on the other side of the skid trail before depositing them parallel to the skid trail.

All three modelling approaches show that increasing intensities of removals lead to an increased probability for bark damage. This result was not unexpected, as it concurs with empirical knowledge of forest practitioners and research studies (Nill et al. 2011; Han and Kellogg 2000; Meng 1978). However, the consistency of this result throughout all three modelling approaches should be emphasised as well as the fact that thinnings which remove >30 % are associated with high probabilities of bark damage.

The effect of the influence of the mean pre-skidding distance (up to the skid trail) on bark damage shows a similar trend. Increasing pre-skidding distances led to higher probability of bark damage. In particular, when the mean pre-skidding distance exceeded 20 m, there was a clear trend of rising probability of newly inflicted bark damage, but it should be acknowledged that in this study, there was only one combination of harvesting system with a mean pre-skidding distance of 20 m, and for this reason maybe, there are limitations to our results. But the relationship described is a clear indication not only from our modelling approaches but was also observed in several studies like the one by Froese and Han (2006) of CTL thinning in mixed conifer stands in Idaho, USA.

The variable “harvesting systems” was one of the most influential variables for the occurrence of bark damage throughout all three models. Harvesting systems that include the use of a cable yarder or a cable skidder (cable winch only) for hauling logs were clearly identified to cause more bark damage than skidders using a crane with a grip tong for smoothly manipulating logs between the trees remaining in the stand. This principal finding is confirmed by various case studies including forest machines combining forwarder and classical hauling solutions with strong cranes and grip tongs. Also Stuhlmann and Findeisen (2009) as well as Sauter et al. (2004) emphasise that these advantages become even more obvious under difficult conditions like with large trees, steep terrain, mixed stands or stands with already established advanced regeneration. Likewise, systems which harvest tree-length logs or full-trees as the main assortment appeared more detrimental to the remaining stand than systems harvesting mainly short-length logs. Furthermore, fully and/or partially mechanised harvesting operations were generally associated with low damage levels. Herewith, the results of the models underline results from several harvesting studies from different countries like Austria (Limbeck-Lilienau and Stampfer 2004; Limbeck-Lilienau 2003; Stampfer et al. 2002), Oregon, USA (Han and Kellogg 2000), Maine, USA, (Ostrofsky et al. 1986) as well as from Germany (Sauter and Grammel 1996). The higher occurence of bark damage was obviously caused by the length of the logs or whole trees independent of the harvesting machine used.

Tree height does not show a clear impact on the occurrence of bark damage. The variable “tree species” in GLMM and MCMCglmm show a significant global effect on the occurrence probability of bark damage. Comparison between tree species indicates that spruce and beech appeared equally vulnerable to suffering bark damage, whereas fir, oak, pine and larch were significantly less vulnerable. These differences between the species could be because of their different bark characteristics (e.g. thin bark in spruce and beach).

Finally, the results of the three introduced modelling approaches and their results provided important information for optimising Nill’s prototype GLMM-based bark damage prognosis model (Nill 2011; Nakou et al. 2013). Specifically, the predictor variable “sampling point is located on skid trail”, representing in Nill’s model proximity to a skid trail, displayed significant collinearity in the new model with the variables “mean skidding distance” and “distance of tree to skid trail” and was therefore subsequently removed. The improved modelling approach made the relevance of the chosen harvesting systems to the occurrence of bark damage evident and so far provides better fits with higher AUC value and adjusted R ² than that of Nill (2011). According to the new model, up to 85 % of the bark damage can be correctly classified instead of only 76 %.

Thus, we believe the final GLMM model presented in this study includes the relevant major factors impacting on the infliction of bark damage during practical harvesting operations. Furthermore, as the model allows for assessing the magnitude of different factors’ impacts in a quantitative manner, it provides managers with a rational and conclusive tool for optimising harvesting operations with respect to bark damage. Thus, the current situation, which appears rather unsatisfactory, might be advanced by selecting, designing and/or applying improved harvesting operations.