The Probability of Moose Damage at the Stand Level in Southern Finland

The probability of moose damage was studied in sapling stands and young thinning stands in southern Finland. Data from the eighth National Forest Inventory in 1986–92 were used for modelling. The frequency of damage was highest at the height of two to fi ve meters and at the age of ten to twenty years (at the time of measurement). Moose preferred aspen stands the most and least preferred Norway spruce stands. Scots pine and silver birch were also susceptible to damage. Logistic regression models were developed for predicting the probability that moose damage is the most important damaging agent in a forest stand. The best predictive variables were the age and dominant species of the stand. Variables describing the site were signifi cant as cluster averages, possibly characterizing the area as a food source (fertility and organic soil), as well as the lack of shelter (wall stand). When sample plot, cluster and municipality levels were compared, it was found that most of the unexplained variance was at the cluster level. To improve the model, more information should be obtained from that level. The regression coeffi cients for aspen as supplementary species, and for pine as dominant species, had signifi cant variance from cluster to cluster (area to area). It was also shown that the occurrence of aspen is closely connected to the occurrence of moose damage in pine sapling stands.


Introduction
Moose (Alces alces L.) and other cervids cause extensive damage to forestry by browsing tree saplings.On the other hand, they have substantial recreational value via hunting.During recent years, forest owners have demanded reductions in the density of moose because of the severe damage in forest regeneration areas.The damage compensation they receive from the hunting licence fees is considered inadequate.According to the Finnish National Forest Inventory, cervid damage had an occurrence of 1.3 percent on forest land (Statistical yearbook … 1999).The frequency of damage in forest regeneration areas is at least ten-fold, and some damaged trees Silva Fennica 35(2) research articles remain in the forest even in the fi rst commercial thinning (Heikkilä 2000).
In this study, moose and other cervids are called 'moose' since it is diffi cult to make difference between browsing signs of different cervid species.Moose is also the most common damaging agent in this group.In southern Finland, roe deer (Capreolus capreolus L.) is becoming more common but is a more important damaging agent only in Ahvenanmaa county.Moose browsing in Scots pine (Pinus sylvestris L.) sapling stands is common locally in winter habitats where moose density can be even higher than ten animals per 1000 hectares.In the summer, birch (Betula spp.) sapling stands are also heavily browsed, which prevents the regeneration of birch in some areas.
The effect of the areas adjacent to regeneration areas on the occurrence of moose damage was studied by Repo and Löyttyniemi (1985).Habitat requirements and the feeding behaviour of moose have been recently studied by Heikkilä and Mikkonen (1992) as well as Heikkilä andHärkönen (1993, 1996), and browsing in relation to food alternatives by Härkönen (1998a,b).In this report, the hypotheses for examination were mainly based on these Finnish studies.
The forest stands most susceptible to moose damage are considered to be pine sapling stands of up to four meters high in the winter and birch stands up to six meters high in the summer.The risk is assumed to increase with the density of the moose population, a high percentage of the preferred tree species in the stand, the location of the stand with respect to the winter habitat or migration routes, a small regeneration area close to sheltering Norway spruce (Picea abies (L.) Karst.)stands, a high percentage of peatlands in the area, and an advanced growth of deciduous species compared to pine saplings.The risk of browsing may be smaller in summer in pure pine stands and in the winter in mixed stands of high density.In the summer, moist areas, such as peatlands, are favored and in the winter dryer areas.
In forest management planning, quantitative estimates of the risks of forest damage associated with management alternatives would be useful.The aim of this study has been (i) to identify conditions which make forest stands vulnerable for browsing by moose, (ii) to quantify the effect of individual factors, and (iii) to develop modelling methods for studying the variability of moose damage on different scales of observation.The main focus concerns the internal properties of the stand and how they affect the probability that moose damage is the most important damaging agent in a forest stand.

Material
The data used were from the eighth National Forest Inventory (NFI8), collected during the period 1986-1992 in southern Finland.In the NFI8, southern Finland includes Ahvenanmaa county and the 15 most southern forestry boards in Finland, situated approximately between 60°N and 64°N (Fig. 1).In this area, measurements were made in rectangular clusters, each of which contained 21 sample plots with a distance of 200 m between sample plots.The distance between the clusters was 8 km in north-south and 7 km in east-west direc tion.For descriptions of the NFI-data, see e.g.National Forest Inventory fi eld manual (1991), Nevalainen (1999) or Nuutinen et al. (2000).
For each sample plot, variables were measured at tree and stand levels.In this study, only the stand-level variables were used.The stand description included species proportions, age, basal area, density, the stage of stand development and the status of management actions.The growing conditions on the plots were characterized by variables describing site type (fertility class) and water balance (drainage).
Moose damage occurs mostly as browsing and breakage of branches or terminal shoots.In the NFI8, the economically most important forest damaging agent has been recorded for each forest stand both in the dominant crown layer and in the understorey.Therefore, the recorded damage is not always the most common in the stand.Most moose damage recordings probably concerned the last fi ve years as the signs of moose damage should disappear rapidly as trees grow.
The severity of the damage has been divided into four classes, from mild to total damage.In this study, all sample plots with moose damage in the dominant crown layer were grouped together regardless of the severity of damage, since the aim was to predict the overall occurrence of damage.Most (52 percent) of the detected damage cases were classifi ed as 'mild damage' without reduced stand quality, 38 percent was 'detectable damage' and only 10 percent were classifi ed as 'severe' or 'total damage'.The division into degrees of damage is subjective, depending on the stand density and also other causes of damage besides moose damage, and was thus not utilized in the analysis.Plots with no damage were included as basis for comparison, and plots with unidentifi ed damage were excluded.The modelling data set was defi ned on the basis of the development class of the stand: only the plots in small and advanced sapling stands and young thinning stands were included.Under fi ve years old regeneration areas were also omitted because of the rarity of damage.Thus the basic data set (base population) for damage models consisted of small (dominant height less than 1.3 m) and advanced sapling stands (diameter less than 8 cm) and young thinning stands.
The estimated moose densities for moose economic areas were received from the Finnish Game and Fisheries Research Institute for the years 1985-1990.The densities have been estimated on the basis of observations made by hunters each autumn.The number of moose per 1000 hectares of forest land in each municipality, for the year preceding the forest inventory year, was used as a variable in the model.Thus the same density number was used for all municipalities in a moose economic area, covering a few municipalities.

Statistical Methods
The forest inventory data has a multilevel structure, consisting of sample plot, cluster, municipality and forestry board levels.It can be assumed that the binary responses of sample plots (occurrence of damage on a plot) are correlated within these levels.When this stucture is taken into account, more accurate (usually larger) estimates for the standard errors of the parameters are obtained (Goldstein 1995).The multilevel model provided also information concerning the residual variation, which is divided into variance components representing the different levels.The plot level residuals are assumed to have binomial variance.
The logistic regression models can be expressed as where y ijk is the observed response of plot i in cluster j in municipality k etc., π ijk is corresponding response probability, and v k and u jk are normally distributed errors at the municipality and cluster levels, respectively.X ijk are the observations and b are the regression coeffi cients.If the municipality level is not signifi cant, the model is simplifi ed and the term v k and subscript k are omitted.Probability estimates are obtained as p = exp (logit) / (1 + exp(logit)).
The method for estimating the logistic regression model for forest damage probabilities has been described in Jalkanen andMattila (1998, 2000).The multilevel approach has been used also by Mattila (2001).Basic analyses were carried out with SAS software (SAS Institute Inc.) and the multilevel model was estimated with MLwiN software (Goldstein et al. 1998).
Explanatory variables for the models were fi rst selected by listing the potential candidates from the literature, and then by calculating cross-classifi cations of the occurrence of damage in different classes of the explanatory variables.Ideally, there is a logistic relationship between the explanatory variable and the relative occurrence of damage, since the logit transformation is then linearly related to explanatory variables.To linearize the effect, a square root transformation was made to the average standwise tree height measurements (sapling stands only), and a logarithmic (base 10) transformation to stand age data in all stands.A basic model with all potential effects was estimated stepwise and explanatory variables were selected using the Wald test (p < 0.05) (Hosmer and Lemeshow 1989).
Fixed effects and random parameters on cluster and higher level (u jk and v k ) were estimated using the fi rst order marginal quasi likelihood (MQL) method of the RIGLS procedure in the MLwiN software.The second order PQL method did not converge with this data.Since most of the residual variation was found to be at the cluster level, averages from that level were also employed as explanatory variables.Finally, a mixed model was examined to see whether the regression coeffi cients varied between locations (clusters).
The effects were also expressed as population averaged odds ratios, which were calculated by exp(β) of the estimated parameter values of the fi xed part of the model (Breslow and Day 1980;Hosmer and Lemeshov 1989).The odds ratio approximates how much the ratio of the frequency of damaged plots to undamaged plots changes when the variable changes by one unit (Diggle et al. 1994).For this interpretation to be valid, the occurrence of damage has to be coded as 0 or 1 for each stand, in which case class 1 contains all degrees of damage.
The model fi t was tested graphically by showing the differences in the mean observed frequencies and mean predicted probabilities in a number of classes of the data, when the probabilities were fi rst sorted into ascending order (Hosmer and Lemeshov 1989).

Base Population and the Basic Model
The proportions of the sample plots damaged by moose in different development classes in the original data were as follows: open area 0.30 percent, small sapling stand 3.79, advanced sapling stand 9.89, young thinning stand 2.37, advanced thinning stand 0.09, mature stand 0.05 percent and none in seed-tree and shelterwood stands.In the modelling data, the total amount of sample plots was 19 145, of which 852 were damaged by moose (4.45 percent).
The frequency of moose damage was highest at the height classes of from two to fi ve meters, and aged from ten to twenty years at the time of measurement (Fig. 2a).The proportion of deciduous trees or the density of the stand did not qualify as predictors (Fig. 2b,c).Of the dominant tree species, moose preferred aspen stands most and spruce stands the least.No difference was apparent between pure and mixed sapling stands (Fig. 2d).Pine and silver birch were also very susceptible to browsing.Downy birch and alder were less desirable.
In the basic logistic model (Table 1, Table 2: model I), variables describing the age and dominant species in the stand had the best predictive value.Due to the methodology, the three main tree species pine, birch and spruce, could not be entered into the model at the same time.In this model, the spruce stand is the basic level to which the other tree species are compared.
Of the other recorded variables in the inventory, the effects of other land use classes besides forest land, main site class spruce mire, and cutting method thinning were not signifi cant although were chosen in cross-classifi cations (cf.Table 1).Other studied features were drainage, silver and downy birch as supplementary species, advance growth of birch compared to pine, height above sea level and the distance of stand edge from the sample plot.Stand mean height worked as well      Cajander (1926).c) Forest management actions like thinnings made during the last 10-year-period are recorded.d) A deciduous species is recorded as supplementary species if it is the most important deciduous species in a coniferous stand.
in the model as stand mean age but it was only measured for sapling stands.
An attempt was made to use moose density as an explanatory variable, but in the basic model with plot level variables, it was only nearly signifi cant (Wald test statistic 3.6 and p = 0.06) and clearly inferior to the variables describing the site.Moose density received an odds ratio of 1.002 (model not shown).Since it did not change the regression coeffi cients for other variables, it was not used in other models.

The Multilevel Model
All variables in the basic model were also tested in the multilevel model context (Table 2, model II).After introducing the random cluster level into the model, the effect of aspen_s (variable names in italics) increased slightly and the regression coeffi cient changed from 1.53 to 1.74, raising the odds ratio from 4.6 to 5.7.Odds ratios are not shown for other than model I since the other changes were even smaller.Standard errors also increased a little.Most of the residual variation was between clusters (u 1jk ), and only 9 percent was between municipalities (v 1k ) in model II.The forestry board level and year were also included but they were not signifi cant.
Since most residual variation was found to be within the cluster level, cluster averages were added as explanatory variables (Table 2, model III).This group of variables (wall, organic, fertility) described the site.Residual variance between municipalities increased to 16 percent in model III.The total residual variation (v 1k + u 1jk ) remained on the same level compared to the previous model with plot level variables.The models behaved correctly in the sense that when the actual occurrence of the damage increased, the predicted probabilities also increased (Fig. 3).
Of the variables tested for randomness at the cluster level in a mixed model, only pine and aspen_s were signifi cant (Table 2, model IV).
The parameters for explanatory variables did not change much with this approach.However, the variances and covariances could be studied more closely.The variance components for the random effects (variance-covariance table) at the cluster level showed that aspen as a supplementary species had signifi cant positive covariance with the cluster level residual (5.70), and negative covariance with pine as dominant species (-9.34).

Discussion
In was an expected result that the best predictive variables were the dominant tree species and the age of individual trees which describe the amount of forage available for moose at a suitable height from the ground.The highest probabilities In the mixed model IV, the coeffi cients for the explanatory variables are also allowed to vary randomly.N.s.not signifi cant, -not estimated.The last 4 rows are the variance and covariance terms: v 1k residual for constant at the municipality level, u 1jk residual for constant at the cluster level, u 3jk residual for pine at the cluster level, u 6jk residual for aspen_s at the cluster level a) .For Wald statistics, see Methods.for moose damage were found in aspen sapling stands since it is a favored food source (Fig. 4).
In addition to aspen, moose is known to prefer rowan, willow and juniper but browsing on these species cannot be predicted due to lack of modelling data.
In NFI, only the economically most important damaging agent is recorded, when in practice two or more damaging agents may occur simultaneously in the same stands.This can cause a systematic underestimation in results compared to frequencies in forests.A lot of moose browsing is directed to less important supplementary species, which may affect the observed differences between susceptibility of different species in the model, as well as the recording method of supplementary species in general (cf.Table 1).Also, the frequency of damage was low in seedlings less than fi ve years old in NFI-data, but moose is known to browse already on small seedlings (Heikkilä and Härkönen 1996).The actual predictive ability of the model with regard to cumulative year-round browsing should be tested with independent data and simulations in planning systems.
First thinning reduced moose damage frequencies signifi cantly but it was omitted since it did not change the parameter estimates of the other explanatory variables.The thinning itself probably did not have any effect on moose damage.Instead, the damaged trees were removed from the plots during thinnings before recording took place.The height of trees should protect them from browsing and breakage at the time of fi rst thinning, but in thinned Scots pine stands, moose may browse on fresh cutting residues (pine crowns).
The situation with aspen as a supplementary species is complicated.In NFI8 in southern Finland, a supplementary species was in coniferdominated forests the most important deciduous species, and in deciduous forests, the most important coniferous species (National Forest Inventory 8, Field manual 1991).This may have an effect on the modelling results.Positive covariance with the constant term in the model was considered to mean if there were a lot of aspen as supplementary species in a cluster, that cluster would be heavily browsed (high cluster averages for both).If aspen was common in a cluster, then pine was less dominant and vice versa (high and low cluster averages).This may be the joint effect of tending, removing aspen, and of pinetwisting rust (Melampsora pinitorqua (Braun) Rostr.), increased by aspen.A large amount of aspen seems to increase both moose damage (this study) and the incidence of pine-twisting rust in a sapling stand (Mattila et al. 2001), and after a few years it may be hard to determine the primary reason for damage.
Moose density turned out not to be a good predictor due to problems in time scale and spatial scale.In the late 1980's, moose densities for game management districts were almost within the recommended level for Southern Finland, usually from three to four animals per thousand hectares, range being two to seven animals.For model improvement, the average density for the last 5 years could be used (with fresh damage recordings) instead of a one-year fi gure since the damage is cumulative.Since only 10 percent of the variance in damage frequency was at the municipality level, more efforts should be put into studying the variation within municipalities, i.e. the location and features of winter habitats and their surroundings where the highest moose densities and the highest damage frequencies occur.The results in this study supported those of Repo and Löyttyniemi (1985) which showed that the surroundings of Scots pine plantations had an effect on moose damage.They found that manmade landscape features reduced the damage, whereas good visibility and the dominance of spruce stands in the area increased it.In the present study, local conditions described by cluster level variables were also signifi cant in explaining moose damage, but the proportion of spruce-dominated plots on the cluster or the effect of other land use classes on the plots were not.However, the study material comprised only of the sample plots on forest land.
In a comparable area in Ontario, Canada, moose tended to choose sheltered closed-canopy coniferous habitats in the winter (Forbes and Theberge 1993).On regional scale, moose were found to favor areas of canopy disturbed by logging or insect damage.Variables from scales of this magnitude, over 100 km 2 , were not studied in the present study, but they might not have improved the model much since most of the variation in damage frequencies was detected at the cluster level.The locations of moose damage could be aggregated at a scale between the cluster and municipality levels (cf.Fig. 1), which may be the same areas in which moose concentrate during the winters.
For predicting probabilities of moose damage, model II or III can be used, depending on the availability of information from the area surrounding the stand, which is necessary for model III.The simpler model II, which has fewer variables, was actually a slightly better predictor than model I and gave slightly higher estimates (cf.Fig. 4).Model III is suited to inventory type data sets with a plot and cluster structure.Since the regression coeffi cients were very close in all the models, the more complicated modelling approaches were not really necessary for predicting moose damage.However, it seemed that the basic model could overestimate probabilities in the highest classes (cf.Fig. 3).
The aims of identifying conditions in which moose damage most often occurs and quantifying their effect were at least partly reached.The presented type of models can also serve as a starting point to combining data from different levels of observation if the studied phenomenon is infl uenced by factors on several levels like forest damage frequencies.The approach helps to direct future studies towards the levels of observation where most of the unexplained variance is.To properly study the appropriate spatial and temporal scale, we need to analyze data where all stands are inventoried in an area, and in the case of moose, where also the local moose populations and their dynamics are known.
With models employing large-scale inventory data, the relationships between the frequencies for different damaging agents are realistic.The predictions may not be meaningful in a single stand but they may be, on average, for a larger area.They can also be helpful in predicting the occurrence of forest damage consequent upon changes in the structure and management of forests.The risk levels associated with forest management plans can therefore be estimated.One problem remains: the consequences of damage at the tree level should also be known and combined with probabilities to provide a risk estimate of expected loss in tree volume or tree quality.

Fig. 1 .
Fig. 1.The location of moose damage in the NFI8 data for the years 1986-94.The size of the circle presents the number of damaged plots of the total of 21 plots in a cluster.The modelling data consists of the basic measurements from the years 1986-92, not including the updated recordings from part of southern Finland in 1994.© Map data: National Land Survey, permission no.MYY/54/99-p.
Fig. 2a.Moose damage by stand height and stand age.Height (rounded to meters): only sapling stands in base population (N = 7574).Age (rounded to 10 year classes): all stands in base population (N = 19 145).

Fig. 2b .
Fig. 2b.Moose damage by the proportion of deciduous species in the total density/volume in percent in base population.

Fig. 3 .
Fig. 3. Mean observed frequencies and predicted probabilities from models I-III in ten classes of predicted probabilities ('risk deciles').

Fig. 4 .
Fig. 4. Examples of predicted probabilities from modelII with aspen, pine or silver birch as dominant species, and from models I-II with pine as dominant species and aspen as supplementary species (P+A).

Table 1 .
The explanatory variables in the models for moose damage, signifi cant in italics.
a) A sample plot is in a 'wall' stand if it is close to the edge of the forest, edge formed by a change in development class or land use class (clearcutting, agricultural fi elds, lakes etc.).b) Medium high fertility class; Myrtillus and Vaccinium type in

Table 2 .
Parameter estimates for the models: β are the regression coeffi cients and S.E.their standard errors.The basic model I is fi xed.Random variance components for residuals are added in model II with plot level variables only and model III with added cluster level variables.
Note that the terms v 1k -u 6jk in models II-IV consist of variance and covariance terms and the column headers β and S.E. are not applicable.