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

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.


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 AR-TEMIS in Canada (Goudie, 1998;Langevin, 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, 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 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.

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 2 ha -1 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. 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 2 ha -1 , while the residual basal area should not be lower than 16 m 2 ha -1 (MRNF, 2004). The percentage of vigorous trees in 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 2 ha -1 before harvesting, while residual basal area should not be lower than 16 m 2 ha -1 . 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 allegha-Predicting tree harvest using a segmented logistic model 141

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 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 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: where x ij is a row vector of explanatory variables and β β 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 β β, 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 threshold of 23 cm in dbh appeared to be a natural choice since it matched the merchantable limit for broadleaved 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: 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: [3a] -For the post-2004 selection cutting: [3b] -For commercial thinning: where N i is the stem density at the plot level (stems ha -1 ) 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 . This statistic follows a Chi-square distribution with 8 degrees of freedom. Large values for X 2 HL 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 Predicting tree harvest using a segmented logistic model 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 ŷ 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 f it, 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 CAP-SIS 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 2 . 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 2 , 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 simu-lations 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 > χ 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)  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 3 ha -1 . 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 3 ha -1 , the earliest thinning scenario led to a production of 176 m 3 ha -1 , the second thinning scenario led to a production of 172 m 3 ha -1 , and the last thinning scenario led to a production of 159 m 3 ha -1 . 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 2 ha -1 (30.0% of the preharvest basal area). A 40-year rotation was necessary to reach the 24-m 2 ha -1 threshold that made the stratum eligible for another selection cutting. This second cut resulted in a removal of 7.3 m 2 ha -1 (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 2 ha -1 , Predicting tree harvest using a segmented logistic model 145 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 f it 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 146 M. Fortin / Forest Systems (2014) 23 (1): 139-152 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 . 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 . 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 managment 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 . The increase in the probabilities of being harvested for this species is clearly a response to the new standards that were Predicting tree harvest using a segmented logistic model 147 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 speciesdependent 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-  del improved the model likelihood by almost 10% . 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.
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 Predicting tree harvest using a segmented logistic model 149  shown in Fig. 5. Because they represent the mathematical 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 simulated. For a single stand, there will always be a difference between the predicted growth and the observed one. A confidence interval around the predicted growth curves could eventually be computed using Monte Carlo techniques. However, such Monte Carlo-simulated uncertainty 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 standlevel 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 2 ha -1 . However, the plot-level residual basal areas ranged from 0 to more than 40 m 2 ha -1 (Table 1). Applying a stand-level criterion such as a target residual basal area of 16 m 2 ha -1 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 standlevel 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 regulation (Fig. 5). For the moment, there is no strict regulation regarding stand composition and diameter distribution following selection cutting or commercial thinning. Should this change in the future, the consistency of the model predictions would need to be checked again. Although the absence of explicit stand-level criteria in the models might seem to be a limitation, I advocate that maintaining the plot-to-plot spatial variability is more realistic than applying a unique standlevel criterion across the plots.

Conclusions
Modelling the individual probabilities of being harvested 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-mulate new alternative silvicultural treatments that have not been implemented at this time since they require data to be fitted on. However, this should not be a limitation. The data required to fit these models are not longterm 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 spite 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. Whether these variables had been available or not in the dataset 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 simulator. 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 unevenaged management might underestimate stand growth.