# Using a segmented logistic model to predict trees to be harvested in forest growth forecasts

Mathieu Fortin

<figure><img src="image_1.png" /></figure>

## To cite this version:

Mathieu Fortin. Using a segmented logistic model to predict trees to be harvested in forest growth forecasts. FOREST SYSTEMS, 2014, 23 (1), pp.139-152. 10.5424/fs/2014231-04824. hal-01195076

HAL Id: hal-01195076

https://hal.science/hal-01195076v1

Submitted on 28 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Using a segmented logistic model to predict trees to be harvested in forest growth forecasts

Mathieu Fortin¹,²*

¹ AgroParisTech. Centre de Nancy. Laboratoire d'Etude des Ressources Forêt-Bois (LERFoB, UMR 1092).
14 rue Girardet. 54042 Nancy, Francy. ² INRA. Centre de Nancy. 1 rue de l'Arboretum.
54280 Champenoux, France

## Abstract

*Aim of the study:* Predicting future harvested trees is a prerequisite to growth forecasts in managed forest stands. While harvest algorithms have been traditionally used, statistical harvest models could be an interesting alternative approach. The objective of this study was to fit statistical harvest models for different partial cutting treatments.

*Area of study:* The study has been carried out in the province of Quebec, Canada.

*Material and methods:* Data from provincial control survey were used to fit harvest models for three different partial cutting treatments. A two-segment logistic modelling approach was used. The harvest models were designed to be compatible with the ARTEMIS growth simulator, which is currently in use in this province.

*Main results:* The results showed that the probability of being harvested is different across the treatments and primarily depends on tree diameter at breast height and species. In selection cutting treatments in particular, trees close to the merchantable limit (e.g., 23 cm in this study) tended to be less frequently harvested than those with smaller or larger diameters, yielding a sinusoidal pattern that was well captured by the segmented approach.

*Research highlights:* Although the models were of average accuracy as indicated by fit statistics, they made it possible to compare different scenarios in terms of productivity and rotation length when coupled with the ARTEMIS growth simulator. Moreover, compatibility requirements between the simulator and the harvest models appeared to be a major limitation in some cases.

**Key words:** harvesting; statistical model; segmented approach; growth forecasts; managed stands.

## Introduction

Growth simulators are widely used in forestry to predict future forest conditions. They make it possible to compare different management scenarios and eventually to select the one that best suits the management goals. In the last two decades, many single-tree growth simulators have emerged and are now commonly used for this purpose. The Forest Vegetation Simulator in the USA (Dixon, 2002), SILVA and BWIN in Germany (Pretzsch et al., 2002; Nagel, 2009), TASS and ARTEMIS in Canada (Goudie, 1998; Fortin and Langelin, 2010, 2012), and MOSES in Austria (Hasenauer, 1994) are some examples.

Single-tree growth simulators easily apply when the only silvicultural option is clearcutting, or as long as there is no partial cutting within the forecast period.

Whenever this happens, forest managers have to face a major issue: which trees are to be harvested and which ones are to be left standing in future partial cuttings? Single-tree growth simulators do not predict these trees by themselves and, consequently, forest managers have to find another way to simulate these partial harvestings.

A common practice has been to design harvest algorithms. First, the harvest rules that underlie a specific treatment are clearly stated and are then turned into an algorithm that aims at reproducing the harvest pattern. For example, a state-of-the-art single-tree selection cutting would consist in harvesting a given number of stems in order to obtain a targeted inverse-J diameter distribution after treatment with a minimum residual basal area and a minimum removal to ensure profitability (Majcen et al., 1990). When dealing with mixed stands, the harvest rules can be further refined to be species-specific. If the growth simulator provides other tree features such as tree vigour and quality,

<footnote>* Corresponding author: mathieu.fortin@agroparistesch.fr
Received: 05-08-13. Accepted: 18-12-13.</footnote>
they can also be part of the algorithm. Finally, a given percentage of stems with smaller diameters is usually discarded in order to reproduce harvesting damage.

Of course, if the harvest algorithm is wrong, the predictions are likely to be biased regardless of the reliability of the growth simulator. Thus, simulating partial harvesting is a major component of growth predictions in managed stands. Considering the importance of this issue, it is surprising that it has been largely overlooked in the scientific literature. Sterba *et al.* (2000), Arii *et al.* (2008) and Thurnher *et al.* (2011) are among the very few studies that have addressed this issue in the last few years. The larger audience and demand for complex growth simulators have probably received more attention than harvest algorithms, which may be seen as a simple forest management issue.

An alternative approach to harvest algorithms is to fit a statistical harvest model to empirical observations. Basically, if permanent-plot data are available, the tree status (cut or alive) can be seen as a binary variable. The individual probability of being harvested can then be modelled through a logistic regression as a function of some tree and plot features such as diameter at breast height (dbh: 1.3 m), tree vigour and quality, plot basal area and density, etc. Whereas the harvest algorithm is most often designed on the basis of theoretical harvest rules, the statistical harvest model is fitted to empirical observations. In the case of departures from the rules, the statistical harvest model can be expected to better represent the way the forest is really harvested. Thurnher *et al.* (2011) recently fitted such a harvest model to be used with the Austrian growth simulator, MOSES.

Regardless of the nature of the harvest procedure, algorithm or statistical model, a major limitation remains the compatibility with growth predictions. In fact, the procedure cannot rely on features that are not provided by the growth simulator. In this context, the objective of this study was to fit statistical harvest models to three major partial cutting treatments in the province of Quebec, Canada. To do this, a segmented logistic model and the control survey that was carried out on public lands in that province were used. These empirical models were designed to be used with ARTEMIS (Fortin and Langevin, 2010, 2012), a single-tree growth simulator that is currently in use in this province. In addition to testing the harvest models, the implications of the coupling with the growth simulator are also discussed.

# Material and methods

## Database

The data that were used to fit the harvest models were taken from the provincial control survey. In the province of Quebec, harvesting on public lands has long been carried out by private companies. To avoid high grading, Quebec's Ministry of Natural Resources (MNR) has set some good practice standards for a large array of silvicultural treatments. Control surveys are carried out by the MNR in cut blocks to check that harvesting was properly done according to these standards.

The survey consists in establishing variable-radius permanent plots using a 2-m² ha⁻¹ basal area factor wedge prism in the cut block a few weeks before harvesting. All trees with a dbh of at least 9.1 cm are tallied from the plot centre in a clockwork direction. For each tree, the species is recorded and the dbh is measured in 2-cm diameter classes. The plot centre is marked and hidden to make sure that the personnel in charge of the harvesting does not behave differently within the plot than in the rest of the cut block. Plots are revisited after harvesting to determine the status (living or cut) of each tree that was previously identified during the pre-harvest measurement. The results of the survey are then compared with the standards associated with the treatment that was planned for this particular stand. If the standards are not met, the private company has to correct for the problem whenever possible or pay a fine otherwise.

The data were recovered from these control surveys for a total of 59,794 plots, which covered the 1988-2005 period. Most plots were established after 1997 and many different treatments were controlled. From this point on, this article will focus on two major treatments: selection cutting and commercial thinning.

Selection cutting has been, by far, the most popular treatment in uneven-aged broadleaved stands since the early 1990s. Between 1994 and 2003, this treatment was carried out on 40,000 ha of public lands annually (Parent and Fortin, 1999, 2001, 2005). Theoretically, the treatment consists in maintaining an inverse-J diameter distribution in the residual stand, although this is not always the case in practice. In terms of good practice standards, the basal area before harvesting should be equal to or greater than 24 m² ha⁻¹, while the residual basal area should not be lower than 16 m² ha⁻¹ (MRNF, 2004). The percentage of vigorous trees in
<table><caption>**Table 1.** Distribution of the plots as a function of the treatment (stem density and basal area ranges appear in parentheses)</caption><thead><tr><th rowspan="2">Treatments</th><th rowspan="2">Number of plots</th><th colspan="2">Before harvesting</th><th colspan="2">After harvesting</th></tr><tr><th>Stem density (stems ha<sup>-1</sup>)</th><th>Basal area (m<sup>3</sup> ha<sup>-1</sup>)</th><th>Stem density (stems ha<sup>-1</sup>)</th><th>Basal area (m<sup>3</sup> ha<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Selection cutting</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>1997-2004</td><td>37,822</td><td>556.6 (10.4-4227.1)</td><td>25.8 (2.0-70.0)</td><td>426.0 (0.0-3474.2)</td><td>18.0 (0.0-68.0)</td></tr><tr><td>Post-2004</td><td>4,916</td><td>564.5 (38.6-2692.4)</td><td>26.4 (6.0-66.0)</td><td>437.3 (0.0-2274.6)</td><td>18.4 (0.0-50.0)</td></tr><tr><td>Commercial thinning</td><td>2,829</td><td>685.4 (32.5-2516.5)</td><td>28.1 (2.0-74.0)</td><td>470.2 (0.0-2142.7)</td><td>18.4 (0.0-56.0)</td></tr></tbody></table>

terms of basal area should also be higher after the treatment and, consequently, non-vigorous stems should be harvested in priority. The cutting cycle is expected to be between 20 and 30 years.

In spite of all these standards, growth on public lands was found to be lower than that of experimental blocks in the early 2000s (Bédard *et al.*, 2002; Guillemette *et al.*, 2013). It appeared that many vigorous stems were incorrectly classified as non-vigorous and therefore harvested, leaving the residual stands depleted and less vigorous than expected. In 2004, a new vigour classification was adopted (Boulet, 2007) and replaced the former vigour classification of Majcen *et al.* (1990). The new classification was stricter and, as a result, it could be expected that selection cutting after 2004 would be different from what it had been until then. To investigate these possible differences, the selection cutting treatment was divided into two sub-treatments: the 1997-2004 and the post-2004 selection cuttings.

Commercial thinning is mostly carried out in coniferous or mixed stands. According to MNR standards,

this treatment is to be done in stands whose site index is equal to or higher than 15 m at 50 years of age (MRNF, 2004). Stand basal area must be between 23 and 33 m² ha⁻¹ before harvesting, while residual basal area should not be lower than 16 m² ha⁻¹. Removal should be between 28 and 35% of the initial basal area. Thinning is from below since the ratio between the mean quadratic diameter after and before the treatment should be equal to or greater than 1.05.

The dataset that was used to model the three treatments, i.e. the 1997-2004 selection cutting, the post-2004 selection cutting, and the commercial thinning, encompassed 45 567 plots. A summary of these data is provided in Table 1.

There were 593,214 trees in the data for a total of 36 species. Some of these species were rather scarce and, consequently, it was necessary to group them in order to obtain more balanced data. Special care was taken to make sure the species grouping was compatible with the ARTEMIS growth simulator. The number of trees and the dbh range in each species group are reported in Table 2. Yellow birch (*Betula alleghaniensis*)

<table><caption>**Table 2.** Dbh range and tree distribution in the species groups</caption><thead><tr><th>Species group</th><th>Number of trees</th><th>Number of harvested trees</th><th>dbh range (cm)</th></tr></thead><tbody><tr><td>Yellow birch (Betula alleghaniensis Britt.)</td><td>111,228</td><td>29,019 (26.09%)</td><td>10-120</td></tr><tr><td>White birch (Betula papyrifera Marsh.)</td><td>17,517</td><td>5,172 (29.53%)</td><td>10-86</td></tr><tr><td>Spruce spp. (Picea spp.)</td><td>23,196</td><td>7,119 (30.69%)</td><td>10-94</td></tr><tr><td>Red maple (Acer rubrum L.)</td><td>33,946</td><td>10,174 (29.97%)</td><td>10-94</td></tr><tr><td>Sugar maple (Acer saccharum Marsh.)</td><td>220,428</td><td>67,014 (30.40%)</td><td>10-114</td></tr><tr><td>Other broadleaved species</td><td>31,371</td><td>7,107 (22.65%)</td><td>10-112</td></tr><tr><td>American beech (Fagus grandifolia Ehrh.)</td><td>45,864</td><td>14,102 (30.75%)</td><td>10-90</td></tr><tr><td>Poplar spp. (Populus spp.)</td><td>11,547</td><td>5,952 (51.55%)</td><td>10-96</td></tr><tr><td>Pine spp. (Pinus spp.)</td><td>19,712</td><td>5,691 (28.87%)</td><td>10-118</td></tr><tr><td>Other coniferous species</td><td>13,960</td><td>3,376 (24.18%)</td><td>10-112</td></tr><tr><td>Balsam fir (Abies balsamea (L.) Mill.)</td><td>38,871</td><td>19,216 (49.44%)</td><td>10-72</td></tr><tr><td>Eastern cedar (Thuja occidentalis L.)</td><td>25,574</td><td>8,114 (31.73%)</td><td>10-116</td></tr></tbody></table>
niensis Britt.) and sugar maple (Acer saccharum Marsh.) were the most abundant species in the dataset, followed by American beech (Fagus grandifolia Ehrh.). About 30% of the trees were cut in each species group, except for balsam fir (Abies balsamea (L.) Mill.) and poplar spp. (Populus spp.) for which the harvest rate was around 50%.

## ARTEMIS growth simulator

ARTEMIS is a distance-independent growth simulator that is currently used in the province of Quebec for forest management purposes. The simulator was largely described in Fortin and Langevin (2010, 2012). Its architecture is simple, with only four dynamic modules that are all based on 10-year growth intervals: a diameter growth module, a mortality module, a recruitment module and a recruit dbh module. Tree heights and volumes are predicted using a height-diameter relationship model (Fortin et al., 2009a) and a volume model (Fortin et al., 2007).

The four dynamic modules of the simulator were fitted to the data of the provincial permanent-plot network. At the end of a particular growth step, the simulator provides a list of trees with the species, dbh, height, volume and status for each one of them. Tree vigour and quality could not be included in the simulator because these variables had only been measured since the early 2000s and, consequently, there were not enough remeasurements to predict these features. In addition to tree characteristics, some plot features also enter into the different dynamic modules that compose the simulator, including plot basal area, stem density and the 1971-2000 mean annual temperature and precipitation.

The simulator is designed to work either in a deterministic or a stochastic fashion. When simulating a stand or a stratum, it computes growth predictions for each sample plot individually. The forecast for the whole stand or stratum is calculated as the mean of the individual plot projections. For more details about the simulator and its components, readers can refer to Fortin and Langevin (2010, 2012).

## Modelling tree harvest

Tree harvest can be considered as a binary random variable. In the aforementioned context, this variable

could be defined as $y_{ij}$ where i and j are the plot and tree indices, respectively. The response variable $y_{ij}$ takes the value of 1 if the tree was harvested, and 0 if it was still alive after harvesting. Logistic models are adapted to the modelling of binary variables. This kind of model can be defined as follows:

$$ \Pr(y_{ij} = 1 | \mathbf{x}_{ij}, \boldsymbol{\beta}) = \frac{e^{\mathbf{x}_{ij} \boldsymbol{\beta}}}{1 + e^{\mathbf{x}_{ij} \boldsymbol{\beta}}} \quad [1] $$

where $\mathbf{x}_{ij}$ is a row vector of explanatory variables and $\boldsymbol{\beta}$ is a column vector of unknown population parameters. It is worth mentioning that model [1] is based on a logit link function (McCullagh and Nelder, 1989), which makes it possible to restrain model predictions within the range [0,1]. The model parameters, i.e., the elements of $\boldsymbol{\beta}$, can be estimated using a maximum likelihood estimator.

A visual check of the data revealed that smaller trees exhibited a different trend than those in larger diameter classes (Fig. 1). There was a u-shaped pattern in smaller diameter classes, followed by a quadratic trend, yielding a sinusoidal-like curve. This complex trend justified the use of a segmented approach. If each segment can fit a quadratic trend, the model would require at least two segments to fit the trend shown in Fig. 1. However, the exact location of the joint between the two segments remained unclear. Selecting the diameter class with the lowest proportion of harvested trees, i.e., 18 cm, was not necessarily a good choice because it would create two segments of unequal complexity. The left segment would be almost linear, whereas the right segment would have two inflexion points. Selecting a diameter class larger than (but close enough to) 18 cm would create two segments with a single inflexion point. In these conditions, the

<figure><img src="image_3.png" /><figcaption>Figure 1. Proportions of harvested trees as a function of 2-cm dbh classes in the dataset.</figcaption></figure>
threshold of 23 cm in dbh appeared to be a natural choice since it matched the merchantable limit for broad-leaved species in the province of Quebec.

To account for the different trends above and below the 23-cm threshold, the continuous variable $d_{ij} = dbh_{ij} - 23$ was created, where $dbh_{ij}$ is the dbh of tree j in plot i (cm), and a binary variable $m_{ij}$, which takes the value of 1 if tree dbh is larger than 23, and 0 otherwise. These two variables were used in the following segmented model:

$$ \Pr(y_{ij} = 1 | \mathbf{x}_{ij}, \mathbf{\beta}) = \frac{e^{\beta_0 + (\beta_1 + \beta_2 m_{ij})d_{ij} + (\beta_3 + \beta_4 m_{ij})d_{ij}^2} + \beta_6 \ln(N_{ij}+1)}{1 + e^{\beta_0 + (\beta_1 + \beta_2 m_{ij})d_{ij} + (\beta_3 + \beta_4 m_{ij})d_{ij}^2}} \quad [2] $$

If $dbh_{ij} = 23$, then $d_{ij} = 0$ and the two segments of model [2] converge to the same value, which is actually the joint between the two segments. Variable $m_{ij}$ makes it possible to consider different trends for the two segments. For the sake of the example, it was assumed that tree dbh was error-free.

Model [2] was independently fit to each treatment, and some additional plot-level variables, namely the basal area and stem density, were also tested. For a particular treatment, the models were compared using the Akaike Information Criterion (AIC, Pinheiro and Bates, 2000, p. 84). The different explanatory variables were kept in the model only if they were significant and if they yielded lower AIC values.

After some preliminary trials based on model [2], the following specific versions were obtained

— For the 1997-2004 selection cutting:

$$ \Pr(y_{ij} = 1 | \mathbf{x}_{ij}, \mathbf{\beta}) = \frac{e^{\beta_0 + \beta_{1,0} + (\beta_{2,0} + \beta_{3,0} m_{ij})d_{ij} + (\beta_{4,0} + \beta_{5,0} m_{ij})d_{ij}^2 + \beta_6 \ln(N_{ij}+1)}}{1 + e^{\beta_0 + \beta_{1,0} + (\beta_{2,0} + \beta_{3,0} m_{ij})d_{ij} + (\beta_{4,0} + \beta_{5,0} m_{ij})d_{ij}^2 + \beta_6 \ln(N_{ij}+1)}} \quad [3a] $$

— For the post-2004 selection cutting:

$$ \Pr(y_{ij} = 1 | \mathbf{x}_{ij}, \mathbf{\beta}) = \frac{e^{\beta_0 + \beta_{1,0} + (\beta_{2,0} + \beta_{3,0} m_{ij})d_{ij} + (\beta_{4,0} + \beta_{5,0} m_{ij})d_{ij}^2 + \beta_6 \ln(N_{ij}+1)}}{1 + e^{\beta_0 + \beta_{1,0} + (\beta_{2,0} + \beta_{3,0} m_{ij})d_{ij} + (\beta_{4,0} + \beta_{5,0} m_{ij})d_{ij}^2 + \beta_6 \ln(N_{ij}+1)}} \quad [3b] $$

— For commercial thinning:

$$ \Pr(y_{ij} = 1 | \mathbf{x}_{ij}, \mathbf{\beta}) = \frac{e^{\beta_0 + \beta_{1,0} + (\beta_{2,0} + \beta_{3,0} m_{ij})d_{ij}}}{1 + e^{\beta_0 + \beta_{1,0} + (\beta_{2,0} + \beta_{3,0} m_{ij})d_{ij}}} \quad [3c] $$

where $N_i$ is the stem density at the plot level (stems ha⁻¹) and is a group species index.

## Model fit and evaluation

When a particular covariate was found to be significant and to yield a lower AIC value, its fit was further assessed using standardised binomial residuals.

The method consists in averaging the predicted probabilities over some classes of the explanatory variables. For example, the predicted values can be averaged by 2-cm dbh classes. Within each class of the explanatory variable, the sum of observed harvested trees can be calculated as well. Let $p$ and $o$ be these averaged predictions and the sum of observed values, respectively. A standardised binomial residual can be derived as $(n^{-1}o - p) / [p(1 - p)n^{-1}]^{1/2}$, where $n$ is the number of observations that served to calculate $p$ or $o$. As $n$ increases, these standardised residuals are expected to approximately follow a normal distribution of mean 0 and variance 1. Roughly, any standardised residual larger than 2 or smaller than -2 might indicate a lack of fit. These standardised residuals are particularly useful to identify unaccounted for trends in the models.

Once the fit of the models was found to be satisfactory, a 10-fold cross-validation was also performed. The plots were split into 10 groups. The models were then fit 10 times, each time omitting a group. At each run, the model predictions were produced for the observations of the omitted group. After the 10 iterations, each observation of the dataset has a prediction as if it had been generated from an independent dataset. It is worth mentioning that this process also made it possible to detect over-parameterisation, which actually resulted in a lack of convergence when one of the 10 groups was omitted. Once the cross-validation proved to be satisfactory, the final models were fit with the whole dataset.

The Hosmer-Lemeshow test (Hosmer and Lemeshow, 2000) is a popular method for assessing the goodness of fit for logistic models. The test consists in ordering the observations by their predicted values. The observations are then divided into 10 even groups. Like the binomial residuals above, the predicted values are averaged within each group. If $o_k$ and $p_k$ are the sum of observed values and the average predicted probability in group $k$, respectively, the Hosmer-Lemeshow statistic

is calculated as $\chi_{HL}^2 = \sum_{k=1}^{10} \frac{(o_k - n_k p_k)^2}{n_k p_k (1 - p_k)}$. This statistic follows a Chi-square distribution with 8 degrees of freedom. Large values for $X_{HL}^2$ indicate a lack of fit. This test was performed on predictions that resulted from the cross-validation.

Another well-known measure of fit for logistic models is the area under the receiver operating characteristic (ROC) curve. The ROC curve consists in
plotting the sensitivity of the model against 1 minus its specificity for the full range of cutpoints between 0 and 1. Let us divide the observations into two groups, a first group $g_1$ with all the observations with a positive outcome, and a second group $g_2$ with the negative outcomes. For any cutpoint $z$, we can consider the event to occur when the predicted probability $\hat{y}_{ij}$ is larger than or equal to $z$. The model sensitivity is calculated as the ratio between the number of observations in $g_1$ for which this condition is met and the total number of observations in $g_1$. Likewise, 1 minus the specificity is calculated as the ratio between the number of observations in $g_2$ for which this condition is met and the total number of observations in $g_2$. Once the ROC curve is plotted for the full range of $z$, the area under the curve (AUC) provides an assessment of the fit. An AUC value of 1.0 represents a perfect fit, while random chance gives an AUC of 0.5 (Lasko *et al.*, 2005). Models with an AUC larger than or equal to 0.8 are considered to be accurate. Readers can refer to Lasko *et al.* (2005) for further details about this fit statistic. The AUC of each model was computed using the predictions from the cross-validation.

## Model implementation

The harvest models were implemented in the CAPSIS platform (Dufour-Kowalski *et al.*, 2012; see http://capsis.cirad.fr/), which already hosts the ARTEMIS growth simulator. The simple interface makes it easy to run simulations either in deterministic or full stochastic mode.

To test different management scenarios, two strata of Forest Management Unit 6152 (46° 51' N, 74° 26' W), which is located northwest of Montréal, Quebec, Canada, were selected. A stratum is defined as a group of similar stands that are not necessarily continuous in space. The two strata were chosen to represent typical coniferous and broadleaved stands of this region. The coniferous stratum consisted of 15 plots, each with an area of 400 m². It was dominated by spruce species (68% of basal area), followed by balsam fir (21% of basal area). Some deterministic simulations were run to test different commercial thinning schedules.

The broadleaved stratum consisted of 32 plots, each with an area of 400 m², and was dominated by yellow birch (35% of basal area) and sugar maple (33% of basal area). These mixed stands are usually managed through selection cutting. We ran deterministic simulations

to determine the appropriate cutting cycle for this stratum using the post-2004 selection cutting harvest model.

## Results

The Hosmer-Lemeshow statistics and the area under the curve computed from the cross-validation are shown in Table 3. All three models exhibited Hosmer-Lemeshow statistics with probabilities $Pr > \chi^2_{HL}$ smaller than 0.01. However, when plotting the average predicted probabilities against the average proportion of events for each group, no major departures could be found. The AUC value ranged from 0.684 to 0.734, with the maximum and minimum values observed for the 1997-2004 selection cutting. These AUC values indicated an average accuracy. In other words, the models were better than randomness but could not be considered as highly accurate. The parameter estimates resulting from the final fit using the whole dataset are reported in Tables 4, 5 and 6.

To illustrate the different trends, predicted probabilities were generated for the three models and the most abundant species in each treatment. For all three species in the selection cutting treatments, the pattern with respect to tree dbh was roughly the same (Fig. 2). The probabilities of being harvested decreased from 10 cm to 20 cm in dbh and increased afterwards. For American beech, there was a major difference between the two selection cuttings (Fig. 2c). As of 24 cm in dbh, the probability of being harvested tended to be higher in the post-2004 selection cutting. The difference increased with dbh. With regard to the effect of the stem density, it induced a slight increase in the probabilities of all three species (Fig. 3). Although this effect was significant, its magnitude remained small.

The predicted probabilities for the commercial thinning treatment exhibited a different pattern than those of the selection cuttings (Fig. 4). Balsam fir

<table><caption>**Table 3.** Hosmer-Lemeshow statistics ($\chi^2_{HL}$) and areas under the curve (AUC) for the different models (the probability $Pr > \chi^2_{HL}$ appears in parentheses)</caption><thead><tr><th>Model</th><th>χ<sup>2</sup><sub>HL</sub></th><th>AUC</th></tr></thead><tbody><tr><td>Selection cutting<br>1997-2004</td><td>67.14 (&lt; 0.0001)</td><td>0.684</td></tr><tr><td>Post-2004</td><td>22.94 (0.0034)</td><td>0.734</td></tr><tr><td>Commercial thinning</td><td>25.59 (0.0012)</td><td>0.704</td></tr></tbody></table>
<table><caption>**Table 4.** Parameter estimates of the 1997-2004 selection cutting model (standard errors appear in parentheses)</caption><thead><tr><th rowspan="2">Species group</th><th colspan="6">Parameter estimates</th></tr><tr><th>β̂<sub>0</sub></th><th>β̂<sub>1,s</sub></th><th>β̂<sub>2,s</sub></th><th>β̂<sub>3,s</sub></th><th>β̂<sub>4,s</sub></th><th>β̂<sub>5,s</sub></th><th>β̂<sub>6</sub></th></tr></thead><tbody><tr><td>Yellow birch</td><td>-0.9805<br>(0.0491)</td><td>0.1328<br>(0.0164)</td><td>-0.0329<br>(0.0172)</td><td>0.0149<br>(0.0014)</td><td>-0.0161<br>(0.0014)</td><td></td><td></td></tr><tr><td>White birch</td><td>-0.6752<br>(0.0674)</td><td>0.0578<br>(0.0272)</td><td>0.1146<br>(0.0317)</td><td>0.0061<br>(0.0024)</td><td>-0.0090<br>(0.0023)</td><td></td><td></td></tr><tr><td>Spruce spp.</td><td>-0.2337<br>(0.0624)</td><td>0.0834<br>(0.0220)</td><td>0.0225<br>(0.0258)</td><td>0.0069<br>(0.0018)</td><td>-0.0081<br>(0.0018)</td><td></td><td></td></tr><tr><td>Red maple</td><td>-0.0492<br>(0.0538)</td><td>0.0664<br>(0.0152)</td><td>0.0306<br>(0.0183)</td><td>0.0062<br>(0.0013)</td><td>-0.0079<br>(0.0012)</td><td></td><td></td></tr><tr><td>Sugar maple</td><td>-0.4718<br>(0.0457)</td><td>0.0821<br>(0.0084)</td><td>0.0172<br>(0.0092)</td><td>0.0089<br>(0.0007)</td><td>-0.0102<br>(0.0007)</td><td></td><td></td></tr><tr><td>Other broadleaved<br>species</td><td>-2.4113<br>(0.0633)</td><td>-0.0150<br>(0.0209)</td><td>0.0855<br>(0.0231)</td><td>0.0037<br>(0.0017)</td><td>-0.0044<br>(0.0017)</td><td></td><td>0.1853<br>(0.0070)</td></tr><tr><td>American beech</td><td>-0.0596<br>(0.0534)</td><td>0.1432<br>(0.0159)</td><td>-0.0609<br>(0.0186)</td><td>0.0119<br>(0.0013)</td><td>-0.0133<br>(0.0013)</td><td></td><td></td></tr><tr><td>Poplar spp.</td><td>0.2362<br>(0.0812)</td><td>0.2207<br>(0.0516)</td><td>-0.1292<br>(0.0546)</td><td>0.0158<br>(0.0047)</td><td>-0.0172<br>(0.0046)</td><td></td><td></td></tr><tr><td>Pine spp.</td><td>-0.2776<br>(0.0975)</td><td>0.1712<br>(0.0690)</td><td>-0.1372<br>(0.0712)</td><td>0.0161<br>(0.0062)</td><td>-0.0164<br>(0.0062)</td><td></td><td></td></tr><tr><td>Other coniferous<br>species</td><td>-0.7741<br>(0.0909)</td><td>0.1094<br>(0.0490)</td><td>-0.0346<br>(0.0519)</td><td>0.0122<br>(0.0041)</td><td>-0.0134<br>(0.0041)</td><td></td><td></td></tr><tr><td>Balsam fir</td><td>1.9240<br>(0.0544)</td><td>0.3418<br>(0.0121)</td><td>-0.2092<br>(0.0209)</td><td>0.0149<br>(0.0009)</td><td>-0.0183<br>(0.0010)</td><td></td><td></td></tr><tr><td>Eastern cedar</td><td>0</td><td>0.0354<br>(0.0253)</td><td>0.0170<br>(0.0278)</td><td>0.0023<br>(0.0022)</td><td>0.0034<br>(0.0022)</td><td></td><td></td></tr></tbody></table>

trees were harvested in priority. The predicted probability sharply increased between 10 and 23 cm in dbh. After the 23-cm threshold, the probabilities flattened out at around 0.90. Compared with balsam fir, spruce species were less frequently harvested on average. The probabilities were rather flat between 10 and 23 cm in dbh. Beyond the threshold, they tended to increase up to 0.5 when the dbh reached 40 cm.

Some management scenarios are compared in Fig. 5. For the coniferous stratum, the absence of thinning led to an asymptotic trend (Fig. 5a). The different thinning schedules resulted in removals between 35 and 48 m³ ha⁻¹. Moreover, the responses to thinning were different: an earlier thinning yielded a positive growth after harvesting, whereas the volume decreased immediately after the last thinning. All volumes considered, i.e., standing and removed volumes, the absence of thinning led to a production of 158 m³ ha⁻¹, the earliest thinning scenario led to a production of 176 m³ ha⁻¹,

the second thinning scenario led to a production of 172 m³ ha⁻¹, and the last thinning scenario led to a production of 159 m³ ha⁻¹. In this particular case, the earliest thinning would be the scenario that maximises the production.

The second example was this broadleaved stratum with high initial basal area (Fig. 5b). With no management at all, the basal area of this stratum was expected to slowly but steadily decrease over time. A first selection cutting at the very beginning of the simulation led to a removal of 8.3 m² ha⁻¹ (30.0% of the pre-harvest basal area). A 40-year rotation was necessary to reach the 24-m² ha⁻¹ threshold that made the stratum eligible for another selection cutting. This second cut resulted in a removal of 7.3 m² ha⁻¹ (30.9% of the pre-harvest basal area). However, the response to this second removal was rather slow. At the end of the projection, i.e., 40 years after the second selection cutting, the stratum had a basal area of 20.7 m² ha⁻¹,
<table><caption>**Table 5.** Parameter estimates of the post-2004 selection cutting model (standard errors appear in parentheses)</caption><thead><tr><th rowspan="2">Species group</th><th colspan="6">Parameter estimates</th></tr><tr><th>β̂<sub>0</sub></th><th>β̂<sub>1,s</sub></th><th>β̂<sub>2,s</sub></th><th>β̂<sub>3,s</sub></th><th>β̂<sub>4,s</sub></th><th>β̂<sub>5,s</sub></th><th>β̂<sub>6</sub></th></tr></thead><tbody><tr><td>Yellow birch</td><td>-1.5033<br>(0.1230)</td><td>-0.0181<br>(0.0479)</td><td>0.1177<br>(0.0504)</td><td>0.0090<br>(0.0043)</td><td>-0.0038<br>(0.0043)</td><td></td><td></td></tr><tr><td>White birch</td><td>-1.0608<br>(0.2034)</td><td>0.0822<br>(0.0891)</td><td>0.1221<br>(0.1076)</td><td>0.0130<br>(0.0076)</td><td>-0.0157<br>(0.0073)</td><td></td><td></td></tr><tr><td>Spruce spp.</td><td>-1.2236<br>(0.2091)</td><td>0.0360<br>(0.0871)</td><td>-0.0190<br>(0.1030)</td><td>0.0030<br>(0.0073)</td><td>-0.0022<br>(0.0071)</td><td></td><td></td></tr><tr><td>Red maple</td><td>-0.5770<br>(0.1357)</td><td>0.1373<br>(0.0485)</td><td>0.0022<br>(0.0566)</td><td>0.0131<br>(0.0041)</td><td>-0.0154<br>(0.0040)</td><td></td><td></td></tr><tr><td>Sugar maple</td><td>-0.6732<br>(0.1055)</td><td>0.1955<br>(0.0226)</td><td>-0.1294<br>(0.0247)</td><td>0.0181<br>(0.0019)</td><td>-0.0183<br>(0.0019)</td><td></td><td></td></tr><tr><td>Other broadleaved<br>species</td><td>-1.5882<br>(0.1793)</td><td>-0.0956<br>(0.0703)</td><td>0.1693<br>(0.0775)</td><td>-0.0004<br>(0.0057)</td><td>0.0000<br>(0.0057)</td><td></td><td>0.0929</td></tr><tr><td>American beech</td><td>-1.6013<br>(0.1709)</td><td>-0.3691<br>(0.1288)</td><td>0.2349<br>(0.0430)</td><td>-0.1171<br>(0.0501)</td><td>-0.0218<br>(0.0034)</td><td></td><td>(0.0207)</td></tr><tr><td>Poplar spp.</td><td>-0.0393<br>(0.2371)</td><td>0.2151<br>(0.1464)</td><td>0.0140<br>(0.1582)</td><td>0.0202<br>(0.0126)</td><td>-0.0253<br>(0.0125)</td><td></td><td></td></tr><tr><td>Pine spp.</td><td>-0.9233<br>(0.4537)</td><td>0.1154<br>(0.2992)</td><td>-0.0995<br>(0.3114)</td><td>0.0213<br>(0.0257)</td><td>-0.0216<br>(0.0256)</td><td></td><td></td></tr><tr><td>Other coniferous<br>species</td><td>-1.8361<br>(0.2435)</td><td>0.0290<br>(0.1225)</td><td>0.0780<br>(0.1304)</td><td>0.0105<br>(0.0100)</td><td>-0.0116<br>(0.0100)</td><td></td><td></td></tr><tr><td>Balsam fir</td><td>3.4432<br>(0.1883)</td><td>0.4345<br>(0.0498)</td><td>-0.4542<br>(0.0891)</td><td>0.0111<br>(0.0034)</td><td>-0.0105<br>(0.0035)</td><td></td><td></td></tr><tr><td>Eastern cedar</td><td>0</td><td>0.2063<br>(0.0640)</td><td>-0.1542<br>(0.0689)</td><td>0.0151<br>(0.0057)</td><td>-0.0156<br>(0.0056)</td><td></td><td></td></tr></tbody></table>

which was not enough to pursue a third selection cutting. In this particular case, even a 40-year cutting cycle did not appear to be sustainable.

## Discussion

Predicting future forest conditions in managed stands is challenging because it involves the coupling of single-tree growth simulations and harvest predictions. Whenever the harvesting is partial, i.e., not clearcut, the selection of the trees to be harvested in growth forecasts has to be predicted in some way. In this study, three statistical harvest models, which stood for two selection cuttings and a commercial thinning treatment, respectively, were fitted to be compatible with the ARTEMIS growth simulator. Two reasons motivated the use of statistical models over harvest algorithms. First, a large dataset was available

to fit these harvest models. Second, the statistical approach provides more information about the model uncertainty. Because the model relies on a parametric distribution, the different sources of variability are known and inferences can be made, for example, on the errors in the parameter estimates. This facilitates the implementation of stochastic simulations, which is a great advantage considering that ARTEMIS also implements full stochastic simulations (see Fortin and Langevin, 2012).

It appears that tree dbh and species, as well as stem density in two out of three models, significantly affect the probability of being harvested. Most broadleaved species exhibited this sinusoidal pattern that can be seen in Fig. 1: the probabilities of being harvested were initially higher for smaller diameters, decreased around 18 cm, and then increased again with tree dbh. This pattern appears to be typical of selecting cuttings (Fig. 2 and 3) as opposed to commercial thinning for
<table><caption>**Table 6.** Parameter estimates of the commercial thinning model (standard errors appear in parentheses)</caption><thead><tr><th rowspan="2">Species group</th><th colspan="4">Parameter estimates</th></tr><tr><th>β̂<sub>0</sub></th><th>β̂<sub>1,5</sub></th><th>β̂<sub>2,5</sub></th><th>β̂<sub>3</sub></th></tr></thead><tbody><tr><td>Yellow birch</td><td></td><td>-1.3086 (0.1011)</td><td>-0.0247 (0.0174)</td><td>0.0816 (0.0187)</td></tr><tr><td>White birch</td><td></td><td>-0.8393 (0.1271)</td><td>-0.0170 (0.0174)</td><td>0.1818 (0.0255)</td></tr><tr><td>Spruce spp.</td><td></td><td>-0.2645 (0.1074)</td><td>-0.0026 (0.0122)</td><td>0.0777 (0.0168)</td></tr><tr><td>Red maple</td><td></td><td>-0.3257 (0.1166)</td><td>-0.0447 (0.0111)</td><td>0.1104 (0.0180)</td></tr><tr><td>Sugar maple</td><td></td><td>-0.4900 (0.1083)</td><td>-0.0079 (0.0143)</td><td>0.0561 (0.0172)</td></tr><tr><td>Other broadleaved species</td><td>-0.9458 (0.0741)</td><td>-0.5316 (0.1070)</td><td>-0.0093 (0.0157)</td><td>0.0683 (0.0185)</td></tr><tr><td>American beech</td><td></td><td>-0.5061 (0.1935)</td><td>-0.0158 (0.0288)</td><td>0.0657 (0.0393)</td></tr><tr><td>Poplar spp.</td><td></td><td>0.4036 (0.1084)</td><td>0.0541 (0.0198)</td><td>-0.0062 (0.0221)</td></tr><tr><td>Pine spp.</td><td></td><td>-0.2306 (0.0852)</td><td>0.0300 (0.0134)</td><td>-0.0138 (0.0141)</td></tr><tr><td>Other coniferous species</td><td></td><td>0.4489 (0.2186)</td><td>0.1121 (0.0755)</td><td>-0.1024 (0.0789)</td></tr><tr><td>Balsam fir</td><td></td><td>3.3417 (0.1167)</td><td>0.2600 (0.0104)</td><td>-0.2117 (0.0283)</td></tr><tr><td>Eastern cedar</td><td></td><td>0</td><td>0.0306 (0.0197)</td><td>-0.0045 (0.0214)</td></tr></tbody></table>

which the probabilities exhibited quadratic trends (Fig. 4).

In selection cuttings, the higher probabilities for larger trees are not surprising. First, large trees tend to be less vigorous on average than smaller trees (Pothier et al., 2013). Since the current regulation aims at maintaining and even increasing the stand vigour, non-vigorous trees are usually harvested in priority. Secondly, large trees have more quality volume on average than smaller trees (Fortin et al., 2009b) and consequently they tend to be more harvested to ensure the profitability of this treatment.

On the other hand, the increase in the probabilities for the smaller trees cannot be explained by profitability or non-vigorous trees. In fact, this trend seems to be the result of the spatial distribution of the trees. Heavy machinery such as harvesters, forwarders and skidders cannot reach larger trees without cutting or damaging smaller ones along their path (Delisle-Boulianne, 2013). This hypothesis has been confirmed in other selection cuttings where most damaged trees had a dbh smaller than 23 cm (Lamson et al., 1985; Tavankar et al., 2013). These losses seem to be reflected in this increase in the small diameters.

Stem density was found to be a significant explanatory variable in the two selection cutting models. Although the effect was rather small, it indicated that a higher stem density increased the probability of being harvested. This effect can be explained by “collateral damage”. In fact, as the stem density increases, there is a higher probability for a tree close to or on the path of a marked tree to be harvested. Thereby, harvesting

the marked tree implies that some other unmarked trees need to be harvested as well.

The sinusoidal pattern was not observed in the commercial thinning treatment (Fig. 4). Because the thinning is theoretically prescribed from below, a higher probability in smaller diameters was expected. A bias towards profitability might be reflected in the increasing probabilities associated with large dbh. In terms of species, there is a clear trend to favour spruce over fir which is in accordance with the management policy on public lands. As a matter of fact, balsam fir is a short-lived species that is vulnerable to the well-known spruce budworm (Choristoneura fumiferana) (Frank, 1990) whereas spruce species are more resistant and have a longer life span (Blum, 1990; Nienstaedt and Zasada, 1990; Viereck and Johnston, 1990). In order to reduce the vulnerability of coniferous stands to spruce budworm outbreaks, it is good practice to harvest balsam fir trees in priority during commercial thinning.

In selection cuttings, the species-specific trends were more obvious in the post-2004 treatment. Whereas the three species showed a similar trend before, the adoption of a new vigour classification in 2004 (see Boulet, 2007) resulted in an increase in the probability for American beech (Fig. 2c) and a decrease in the probability for the largest stems of the other two species (Fig. 2a,b) to be harvested. As a result of its thin bark and its sensitivity to bark disease, American beech is a species that tends to be less vigorous on average than yellow birch and sugar maple, all things being considered (Pothier et al., 2013). The increase in the probabilities of being harvested for this species is clearly a response to the new standards that were
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Predicted probabilities of being harvested for yellow birch (a), sugar maple (b), and American beech (c) (stem density is set to 500 stems ha⁻¹; grey line: 1997-2004 selection cutting; black line: post-2004 selection cutting).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Predicted probabilities of being harvested for selection cutting after 2004 as a function of three stem densities for yellow birch (a), sugar maple (b), and American beech (c) (grey line: 1,000 stems ha⁻¹; dashed line: 500 stems ha⁻¹; and black line: 250 stems ha⁻¹).</figcaption></figure>

adopted in 2004: non-vigorous stems have been marked and harvested in priority since that time.

Such statistical harvest models have also been fitted in Europe. As mentioned in the introduction of this paper, Thurnher *et al.* (2011) fitted a harvest model for Austrian uneven-aged stands using a logistic model with a logit link function. The results of this study are similar to theirs. The same increase in the probabilities could be observed for smaller diameters, and the magnitude of the probabilities was species-dependent as well. This study differs from theirs in that a segmented approach was used and different models were fitted for different partial cuttings. The segmented approach was clearly necessary since the trend

shown in Fig. 1 could not be fitted with a classical logistic regression.

In spite of all the predictors that were included in the models, the AUC values indicated average accuracy (Table 3). Two major reasons explain this result. First, tree selection in partial harvesting is not only based on dbh and species, but on many other characteristics as well such as tree vigour and quality. In other words, two trees with the same dbh and species in the same sample plot may be different with respect to these features and it may be clear that one is to be harvested in priority. In a context of selection cutting, the inclusion of tree vigour and quality in a harvest mo-
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Predicted probabilities of being harvested for balsam fir (black) and spruce species (grey) in the commercial thinning treatment.</figcaption></figure>

del improved the model likelihood by almost 10% (Delisle-Boulianne, 2013). The other reason is that harvesting is a spatial phenomenon by definition. Including tree coordinates as was done in some other studies (e.g., Arii *et al.*, 2008; Fortin *et al.*, 2013) would have improved the accuracy of the predictions. Regardless of whether it was the unconsidered tree features or tree coordinates, there was clearly an issue related to the compatibility with the growth simulator.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Two examples of harvesting simulations: a coniferous stand with different commercial thinning schedules (a) and a broadleaved stand under selection cutting management (b).</figcaption></figure>

In this case study, these variables could not be included in the harvest models because the simulator did not provide them.

It appears that this lack of information did not greatly affect treatments such as commercial thinning for which other tree features may not be so important. The growth forecasts in the coniferous stratum (Fig. 5a) seemed to be consistent with previous studies: as the stand ages, the response to thinning decreases (Pretzsch, 2009, p. 410). However, the growth simulations for the broadleaved stratum were of greater concern. In that case (Fig. 5b), even a 40-year cutting cycle was not sustainable. In the province of Quebec, tree marking in selection cuttings heavily relies on tree vigour (Majcen *et al.*, 1990; Boulet, 2007). Tree vigour is also known to affect tree mortality as well as diameter growth of broadleaved species (Fortin *et al.*, 2008a,b). According to the standards for this treatment, the residual stands should have a higher proportion of vigorous stems than before selection cutting. Consequently, after harvesting, it could be reasonably expected that the probabilities of mortality would decrease, whereas the diameter increments would increase. In its current form, the growth simulator cannot account for any change in the proportion of vigorous trees. Because the mortality module in ARTEMIS does not include tree vigour, the trees keep dying at the same rate after harvesting, all other things being considered. Knowing that the probabilities of mortality for non-vigorous trees are about six times those of vigorous trees for sugar maple (Fortin *et al.*, 2008a), changing the proportion of vigorous stems in a sugar maple-dominated stand or stratum cannot be taken lightly in this case. This remains a major limitation of this coupling between harvest models and ARTEMIS.

Of course, the reader can rightfully wonder why tree vigour has not been taken into account in the simulator. In fact, modelling the evolution of tree vigour requires permanent-plot data in which this variable was measured. Tree vigour has only been assessed since the early 2000s in the provincial permanent plot network that served to fit the component of ARTEMIS. We can expect ARTEMIS to include a vigour and a quality module in the next few years since the plots will be revisited, but for the time being, these features cannot be accounted for.

Such harvest models actually aim at representing the mathematical expectation of stochastic phenomenon, namely tree harvesting. As such, their purpose is not to guide the harvest in the field. They rather serve to predict future harvests in growth simulations as
shown in Fig. 5. Because they represent the mathema-
tical expectation of the harvest, they should be used at
a strategic level for large scale predictions. Under the
assumption that both the harvest models and the growth
simulator are unbiased, the curves shown in Fig. 5 are
actually those we would observe on average if we had
a large number of stands like those that were simula-
ted. For a single stand, there will always be a differen-
ce between the predicted growth and the observed one.
A confidence interval around the predicted growth cur-
ves could eventually be computed using Monte Carlo
techniques. However, such Monte Carlo-simulated un-
certainty needs to be tested (cf. Fortin et al., 2009c) and
this issue alone could be the topic a new study.

The harvest models that were fitted in this study apply at the plot level. They do not include any features that ensure a particular composition or diameter distribution at the stand level. The implementation of such stand-level criteria is complex and may eventually lead to levelling off the spatial variability across the plots. According to the current regulation, the residual stand basal area following selection cutting should not be lower than 16 m² ha⁻¹. However, the plot-level residual basal areas ranged from 0 to more than 40 m² ha⁻¹ (Table 1). Applying a stand-level criterion such as a target residual basal area of 16 m² ha⁻¹ at the plot level would imply a huge loss of spatial variability. The same rational applies for stand composition and diameter distribution.

The consistency of the models with respect to stand-
level criteria can only be checked a posteriori. In this
case study, the residual basal areas as predicted by the
harvest models were consistent with the current regu-
lation (Fig. 5). For the moment, there is no strict re-
gulation regarding stand composition and diameter dis-
tribution following selection cutting or commercial
thinning. Should this change in the future, the consis-
tency of the model predictions would need to be chec-
ked again. Although the absence of explicit stand-le-
vel criteria in the models might seem to be a limitation,
I advocate that maintaining the plot-to-plot spatial va-
riability is more realistic than applying a unique stand-
level criterion across the plots.

Conclusions

Modelling the individual probabilities of being har-
vested through a segmented logistic regression has the
advantage of better representing the harvesting as it is
done in the field. Such models are not designed to si-

ulate new alternative silvicultural treatments that ha-
ve not been implemented at this time since they requi-
re data to be fitted on. However, this should not be a li-
mitation. The data required to fit these models are not long-
term monitoring data. After the implementation of a new
silvicultural treatment, a control survey can provide enough
data on the short term to fit a new harvest model.

The trend in the harvest probabilities may be hard to
fit with a single quadratic trend, as shown in Fig. 1.
Using a two-segment logistic model, as was done in this
study, might be a better alternative in that case. In spi-
te of all the explanatory variables that were included in
the models, their accuracy was not as high as expected.
In fact, the tree marking is often based on their vigour
and quality. Unfortunately, it was impossible to take
these features into account in our harvest models. Whe-
ther these variables had been available or not in the da-
taset was of lesser importance, knowing that the AR-
TEMIS growth simulator did not consider them and
that the harvest models had to be used with this simu-
lator. The monitoring of quality and vigour, as currently
assessed in the province of Quebec (see Boulet, 2007),
will eventually provide the data for fitting vigour and
quality modules in ARTEMIS and for including these
variables in our harvest models. In the meantime,
growth forecasts in broadleaved stands under uneven-
aged management might underestimate stand growth.

Acknowledgements

Daniel Pin, Gordon Weber and Gyna Gagnon from
*Le Bureau du Forestier en Chef du Québec* provided
the database and support for this study. Special thanks
are due to Catherine Collet (LERFoB-INRA Nancy)
and Jean-Pierre Saucier (DRF-MRNQ) for their com-
ments on a preliminary version of this paper. The
author is grateful to the *Direction de la recherche fo-
restière du Ministère des Ressources naturelles du Qué-
bec* who provided the funding for the implementation
of the harvesting models in the CAPSIS platform. The
UMR 1092 LERFoB is supported by the French Na-
tional Research Agency through the Laboratory of Ex-
cellence ARBRE (ANR-12-LABXARBRE-01).

References

Arii K, Caspersen JP, Jones TA, Thomas SC, 2008. A selec-
tion harvesting algorithm for use in spatially explicit in-
dividual-based forest simulation models. Ecological Modelling 211, 251-266.

Bédard S, Brassard F, Boulet B, Gagnon R, Majcen Z, Pinard S, Vézina S, Blais L, Boulay P, Frenette J-M, 2002. Les effets reels des coupes de jardinage dans les forêts publiques du Québec en 1995 et 1996. Gouvernement of Québec. Ministère des Ressouces naturelles. Direction de la planification et des communications. Rapport. 15 pp.

Blum BM, 1990. Red spruce. In: Silvics of North America, Vol 1. USDA Forest Service Agriculture Handbook no 654. 877 pp.

Boulet B, 2007. Défauts et indices de la carie des arbres: guide d'interprétation, 2ᵉ éd. Les Publications du Québec, Québec, Canada.

Delisle-Boulianne S, 2013. Mise au point d'un modèle de prélèvement par tiges individuelles pour simuler les coupes de jardinage en Outaouais. Université Laval, Québec, Canada. Mémoire de maîtrise. 27 pp.

Dixon GE, 2002. Essential FVS: a user's guide to the Forest Vegetation Simulator. Internal Report Fort Collins, CO: USDA Forest Service. Forest Management Service Center. 226 pp. [revised January 2013].

Dufour-Kowalski S, Courbaud B, Dreyfus P, Meredieu C, De Coligny F, 2012. Capsis: an open software framework and community for forest growth modelling. Annals of Forest Science 69: 221-233.

Fortin M, Langevin L, 2010. ARTEMIS-2009: un modèle de croissance base sur une approche par tiges individuelles pour les forêts du Québec. Gouvernement of Québec (Canada), Ministère des Ressources naturelles et de la Faune, Direction de la recherche forestière. Mémoire de recherche forestière no 156. 48 pp.

Fortin M, Langevin L, 2012. Stochastic or deterministic single-tree models: is there any difference in growth predictions? Annals of Forest Science 69: 271-282.

Fortin M, DeBlois J, Bernier S, Blais G, 2007. Mise au point d'un tarif de cubage général pour les forêts québécoises: une approche pour mieux évaluer l'incertitude associée aux prévisions. The Forestry Chronicle 83: 754-765.

Fortin M, Bédard S, DeBlois J, Meunier S, 2008a. Predicting individual tree mortality in northern hardwood stands under uneven-aged management in southern Québec, Canada. Annals of Forest Science 65: 205.

Fortin M, Bédard S, DeBlois J, Meunier S, 2008b. Accounting for error correlations in diameter increment modeling: a case study applied to northern hardwood stands in Quebec, Canada. Canadian Journal of Forest Research 38: 2274-2286.

Fortin M, Bernier S, Saucier J-P, Labbé F, 2009a. Une relation hauteur-diamètre tenant compte de l'influence de la station et du climat pour 20 espèces commerciales du Québec. Ministère des Ressources naturelles et de la Faune, Direction de la recherche forestière. Mémoire de recherche forestière no 153, 22 pp.

Fortin M, Guillemette F, Bédard S, 2009b. Predicting volumes by log grades in standing sugar maple and yellow birch trees in southern Quebec, Canada. Canadian Journal of Forest Research 39: 1928-1938.

Fortin M, Bédard S, DeBlois J, Meunier S, 2009c. Assessing and testing prediction uncertainty for single tree-based models: a case study applied to northern hardwood stands in southern Québec, Canada. Ecological Modelling 220: 2770-2781.

Fortin M, Delisle-Boulianne S, Pothier D, 2013. Considering spatial correlations between binary response variables in forestry: an example applied to tree harvest modelling. Forest Science 59: 253-260.

Frank RM. 1990. Balsam fir. In: Silvics of North America, Vol 1. USDA Forest Service Agriculture Handbook no 654. 877 pp.

Goudie JW, 1998. Model validation: a search for the magic grove or the magic model. In: Stand density management conference: planning and implementation, November 6-7, 1997, Edmonton, AB, Bamsey, C. pp: 45-58.

Guillemette F, Gauthier M-M, Lambert M-C, Bédard S, 2013. Effets réels décennaux des coupes de jardinage pratiquées de 1995 à 1999 dans un contexte opérationnel. Gouvernement of Québec. Ministère des Ressources naturelles, Direction de la recherche forestière. Mémoire de recherche forestière no 168. 34 pp.

Hasenauer H, 1994. Ein Einzelbaumwachstumssimulator für ungleichaltrige Fichten-Kiefern- und Buchen-Fichtenmischbestände. Österreichische Gesellschaft für Waldkosystemforschung und Experimentelle Baumforschung, Vienna, Austria. 152 pp.

Hosmer DW Jr, Lemeshow S, 2000. Applied logistic regression, 2nd ed. John Wiley & Sons, Inc, New York, USA.

Lamson NI, Smith HC, Miller GW, 1985. Logging damage using an individual-tree selection practice in Appalachian hardwood stands. Northern Journal of Applied Forestry 2: 117-120.

Lasko TA, Bhagwat JG, Zou KH, Ohno-Machado L, 2005. The use of receiver operating characteristic curves in biomedical informatics. Journal of Biomedical Informatics 38: 404-415.

Majcen Z, Richard Y, Ménard M, Grenier Y, 1990. Choix des tiges à marquer pour le jardinage d'érablières inéquivalentes, Guide technique. Ministère de l'Énergie et des Ressources du Québec, Canada, Direction de la recherche forestière. Mémoire de recherche forestière n° 96. 96 pp.

McCullagh P, Nelder JA, 1989. Generalized linear models, 2nd ed. Monographs of Statistics and Applied Probability 37. Chapman&Hall. New York. 511 pp.

MRNF, 2004. Manuel d'aménagement forestier, 1ère révision. Gouvernement du Québec (Canada), Ministère des Ressources naturelles et de la Faune, Direction de la planification et des communications.

Nagel J, 2009. Waldwachstumssimulation mit dem Java Software Paket TreeGross. Nordwestdeutsche Forstliche Versuchsanstalt. 62 pp.

Nienstaedt H, Zasada JC, 1990. White spruce. In: Silvics of North America, Vol 1. USDA Forest Service Agriculture Handbook no 654. 877 pp.

Parent B, Fortin C, 1999. Ressources et industries forestières. Portrait statistique, édition 1999. Gouvernement du Québec, Ministère des Ressources naturelles et de la Faune.
ne du Québec, Direction du développement de l'industrie des produits forestiers.

Parent B, Fortin C, 2001. Ressources et industries forestières. Portrait statistique, édition 2001. Gouvernement du Québec, Ministère des Ressources naturelles du Québec, Direction du développement de l'industrie des produits forestiers.

Parent B, Fortin C, 2005. Ressources et industries forestières. Portrait statistique, édition 2005. Gouvernement du Québec, Ministère des Ressources naturelles et de la Faune du Québec (Canada), Direction du développement de l'industrie des produits forestiers.

Pinheiro JC, Bates DM, 2000. Mixed-effects models in S and S-PLUS. Springer, New York, USA. 528 pp.

Pothier D, Fortin M, Auty D, Delisle-Boulianne S, Gagné L-V, Achim A, 2013. Improving tree selection for partial cutting through joint probability modelling of tree vigor and quality. Canadian Journal of Forest Research 43: 288-298.

Pretzsch H, 2009. Forest dynamics, growth and yield. Springer-Verlag Berlin Heidelberg, Germany. 664 pp.

Pretzsch H, Biber P, Ďurský J, 2002. The single tree-based stand simulator SILVA: construction, application and evaluation. Forest Ecology and Management 162: 3-21.

Sterba H, Golser M, Moser M, Schadauer K, 2000. A timber harvesting model for Austria. Computer and Electronics in Agriculture 28: 133-149.

Tavankar F, Majnounian B, Bonyad AE, 2013. Felling and skidding damage to residual trees following selection cutting in Caspian forests of Iran. Journal of Forest Science 5: 196-203.

Thurnher C, Klopf M, Hasenauer H, 2011. Forests in transition: a harvesting model for uneven-aged mixed species forests in Austria. Forestry 84: 517-526.

Viereck LA, Johnston WF, 1990. Black spruce. In: Silvics of North America, Vol 1. USDA Forest Service Agriculture Handbook no 654. 877 pp.