Modelling the effect of habitat composition and roads on the occurrence and number of moose damage at multiple scales

We modelled the effect of habitat composition and roads on the number and occurrence of moose (Alces alces L.) damage in Ostrobothnia and Lapland using a zero-inflated count model. Models were developed for 1 km2, 25 km2 and 100 km2 landscapes consisting of equilateral rectangular grid cells. Count models predict the number of damage, i.e. the number of plantations and zero models the probability of a landscape being without damage for a given habitat composition. The number of moose damage in neighboring grid cells was a significant predictor in all models. The proportion of mature forest was the most frequent significant variable, and an increasing admixture of mature forests among plantations increased the number and occurrence of damage. The amount of all types of plantations was the second most common significant variable predicting increasing damage along with increasing amount of plantations. An increase in thinning forests as an admixture also increased damage in 1 km2 landscapes in both areas, whereas an increase in pine-dominated thinning forests in Lapland reduced the number of damage in 25 km2 landscapes. An increasing amount of inhabited areas in Ostrobothnia and the length of connecting roads in Lapland reduced the number of damage in 1 and 25 km2 landscapes. Differences in model variables between areas suggest that models of moose damage risk should be adjusted according to characteristics that are specific to the study area.


Introduction
Moose (Alces alces L.) population densities in Fennoscandian countries have been probably the highest in the world since 1970s (Lavsund et al. 2003).High moose population has led to substantial moose damage in forests and the amount of damage has remained high until the 2000s (Lavsund et al. 2003). In Finland in 2004-2008, moose damage was recorded in a total of 741 000 ha corresponding to 19% of all plantations and 24% of pine-dominated plantations (Korhonen et al. 2010).Severe damage was recorded in 3% of pine-dominated plantations.In Sweden, the amount of fresh (1-3 yrs old) moose damage varied between 12-15% of Scots pine and 9-44% of birch plantations during the years 2004-2013(Swedish Statistical Yearbook… 2013).
Moose cause damage to trees by browsing leader and lateral shoots which leads to growth losses and impaired timber quality (Lavsund 1987;Bergqvist et al. 2001;Wallgren et al. 2014).The majority of damage occurs in winter when food sources other than trees are under snow cover.Moose browse on several tree species, such as birch (Betula pendula Roth and B. pubescens Ehrh.), aspen (Populus tremula L.), willows (Salix spp.) and juniper (Juniperus communis L.) (Bergström and Hjeljord 1987;Månsson et al. 2007), but in economic terms, Scots pine (Pinus sylvestris L.) is the most important tree species affected by moose (Lavsund 1987).
Moose have annual home ranges varying from some hundreds of hectares to some thousands of hectares (Sweanor and Sandegren 1989;Cederlund and Sand 1994).At least the northern part of the moose population in Fennoscandia is migratory (Ball et al. 2001), and average distances between summer and winter home ranges have been measured to range from a few kilometres to some tens of kilometres (Sweanor and Sandegren 1988;Ball et al. 2001).Moose have been shown to be selective in their home range selection (Nikula et al. 2004;Cassing et al. 2006), and there is evidence that moose are selective also within their home ranges (Cederlund and Okarma 1988;Ball et al. 2001;Nikula et al. 2004).As a consequence, when assessing factors influencing moose damage risk landscape sizes corresponding to the home range size of moose should be addressed.
In summer moose are capable of using a large variety of habitats due to their ability to browse on numerous plant species (Cederlund et al. 1980).In winter, when fewer plants are available due to snow cover and when deciduous trees are bare, moose are more selective than in summer and use areas characterised by a mosaic of young and older forests (Bjørneraas et al. 2012).Stands with an age class of less than 30 years are usually preferred by moose (Månsson et al. 2011) and especially pine-dominated young stands seem to characterise moose winter pastures (Nikula et al. 2004).From the point of view of moose damage risk, it can thus be hypothesized that if preferred habitats are not evenly distributed at the landscape or regional level, it leads to uneven use of habitats and, consequently, uneven distribution of damage.
Studies of moose damage at national and regional levels have shown that changes in moose populations are also reflected in changes in browsing of preferred tree species and damage to pine stands (Hörnberg 2001a).Moose density index has been found to be positively correlated with browsing pressure on Scots pine also at the local scale but no evidence of density dependent habitat selection (Månsson 2009) nor correlation between population level and damage level at regional scale (Hörnberg 2001b) have been found.Instead, significant correlation between the availability of browse species and browsing intensity of the most preferred food sources have been found at regional (Hörnberg 2001a) and at landscape levels similar to or larger than the home range size of moose (Cassing et al. 2006).In addition, moose habitat selection often includes trade-off situations where moose behaviour is constrained by e.g.human activity (Bjørneraas et al. 2011).Therefore, in addition to the availability of preferred habitats also disturbance factors potentially affect moose damage risk.
Qualitative factors affecting moose resource selection and damage at the local, home range and landscape levels are fairly well known, but only a few studies have assessed damage risk quantitatively.The aim of this study was to develop models for the occurrence and number of moose damage by using located moose damage data as a response variable and habitat composition as well as other landscape features as explanatory variables.In this study, we concentrate on factors ranging from the local to landscape and regional levels influencing the occurrence and number of damage.

Study areas
The study was conducted in Finland in southern Lapland and western Northern Ostrobothnia (later referred to as Ostrobothnia, Fig. 1).The northernmost border of the study area in Lapland approximately follows the line where mountainous areas and non-private land start to prevail.To the east and west, the Lapland study area is delimited by the Finnish-Russian and Finnish-Swedish  1a,b).Grey areas on the index map show the location of privately owned land in the study area and bolded lines delineate the provinces of Ostrobothnia and Lapland.borders, respectively.In Ostrobothnia, the easternmost border of the study area is approximately delineated by a topographic border formed by the highest shoreline after the last glacial period (200-220 m a.s.l., Ojala et al. 2013), and the Gulf of Bothnia in the west.The area is restricted by the provincial southern border of Ostrobothnia and the Lapland border in the north.The main tree species in both study areas are Scots pine, Norway spruce (Picea abies (L.) Karst.) and birches.The shrub layer consists of different willow species and juniper.
The Ostrobothnia study area lies in the middle boreal region (Ahti et al. 1968).The terrain is rather flat, consisting of the mosaic of mineral soils and peatlands, as well as anthropogenically fragmented areas, mostly agricultural land.About 50% of forestry land comprises peatland, of which about 60% have been drained (Finnish Statistical Yearbook… 2011).About 60% of the forest land is owned by private forest owners and most of the forests are commercially managed.The period of permanent snow cover lasts from mid-November to the end of April, and the maximum depth of snow cover is 40-60 cm (Pirinen et al. 2012).
The southern and south-western parts of the Lapland study area belong to mid-boreal region, whereas the eastern and northern parts belong to the northern boreal region (Ahti et al. 1968).About one third of the forestry land is privately owned, and most of the forests are commercially managed.The majority of privately owned land is located along the Kemi and Tornio rivers and their arteries.The period of permanent snow cover lasts from the beginning of October to mid-May, and the maximum depth of snow cover is 60-90 cm (Pirinen et al. 2012).

Moose damage data
Moose damage data consisted of located data on damaged plantations for which forest owners were compensated during the years 2002-2008.In Finland, the state pays compensation for moose damage to private forest owners for growth and quality losses and for possible regeneration costs (Finlex 2014).The area of damage must be ≥0.1 ha and the calculated damage has to exceed 170 euros.Compensation is paid only if the damage is less than three years old and it cannot be paid for the same plantation until three years has passed since earlier compensation.There were total of 2663 plantations in Ostrobothnia and 1287 plantations in Lapland that had been compensated for in 2002-2008 in our data.The total area of compensated plantations varied annually from 377 ha to 1831 ha in Ostrobothnia and from 367 ha to 1237 ha in Lapland during the study period.The average damage area per compensated plantation varied annually between 1.6-7.9ha in Ostrobothnia and 1.1-5 ha in Lapland.

Landscape data
Land use and forest data for the analysis were provided by the Finnish National Forest Inventory (NFI).In Finland, the NFI utilizes Landsat and Spot satellite images concurrently with field plots, as well as digital maps of roads, agricultural land and other non-forest land, to separate non-forest land from forest land.The Multi-Source National Forest Inventory method (MS-NFI) (Tomppo et al. 2008) uses the k-nn method for producing estimates of e.g.timber volume for each tree species and for each pixel corresponding to a 25 × 25 m land area.The MS-NFI data used in the analyses correspond to the year 2005.
Digital maps of timber volume estimates for pine, spruce and deciduous trees, as well as digital maps of fields, settlements, roads, waters and peatlands, were imported to a geographic information system (GIS) as separate layers.For the analysis, we combined timber volume layers and other land use layers to a single land use and cover layer using total timber volume as a proxy of forest development stage (Tomppo et al. 1998).Drained peatlands were separated from non-drained Silva Fennica vol. 53 no. 1 article id 9918 • Nikula et al. • Modelling the effect of habitat composition and roads… peatlands with the aid of digital maps of peatlands and digitized ditches (source National Survey of Land).In the resulting data each pixel can belong to one of 15 classes (Table 1a).Digital maps of roads (source Finnish Transport Agency, FTA) were classified in four road classes according to the road class definitions of the FTA (Table 1a).

Grid based modelling
We superimposed three grids sized 1 × 1 km, 5 × 5 km and 10 × 10 km on the study area and then calculated the number of compensated moose damage plantations in each grid cell and used this as the response variable.Hereafter, the cell sizes are referred to as 1 km 2 , 25 km 2 and 100 km 2 cells and landscapes, respectively.Due to there being several owner groups other than private land owners, we used the following criteria for selecting grid cells for inclusion in modelling: 1) privately owned forest land had to be >50% of the area, 2) there had to be >0.5 ha of privately  2012).For modelling we also calculated the sum of different combinations of original variables (later referred to as sum variables) in each grid cell (Table1b).Both the original as well as the sum variables were used as explanatory variables in modelling.

Zero-inflated negative binomial modelling
The number of non-damage grid cells (zero cells) in the 1 km 2 and 25 km 2 cell sizes were excessive in both study areas (52-96% of the cells) indicating zero-inflated data (McCullagh and Nelder 1989).Also, the mean number of damage per grid cell, i.e. the number of compensated damage was smaller than the variance for all cell sizes, which makes our data possibly overdispersed (McCullagh and Nelder 1989).Due to these two factors, we assumed the data to fulfil a zero-inflated nega- tive binomial distribution (ZINB) when modelling the occurrence (binary response, compensated damage occurs or not in a grid cell) and the number of moose damage (Zuur et al. 2009).Out of all 10 km grid cells the number of moose damage cells were higher than non-damage cells in both study areas, but to preserve the comparability of the models among cell sizes we used a ZINB-based modelling procedure also for 10 km cells.However, we also constructed models with exactly the same variable combinations for the number of damage by using a negative binomial (NB) distribution, and compared the performance of the ZINB and NB models using a Vuong test (Vuong 1989).Zero-inflated count models are two-component mixture models in which the distribution of counts is represented by two separate components: the zero component models the probability of false zeros and a count component models true zeros and non-zero values (Zuur et al. 2009).Count models that use zero-inflated negative binomial (ZINB) distribution assume that zeros can come from both the binomial process and the count processes.Basically, our data can contain zero (non-damage) grid cells for at least two reasons.First, so-called false zeros (Martin et al. 2005) might occur because a land owner had not applied for compensation, even though moose damage had occurred in the grid cell.Second, zero cells can also be so-called true zeros when there is no moose damage in a grid cell, even though it contains plantations susceptible to damage.Therefore, the probability (Pr) that there is no moose damage Y in a cell i, can be written as Pr(Y i = 0) = Pr(False zeros) + (1 -Pr(False zeros)) × Pr(Count process gives a zero).
In ZINB regression, the response variable Y i (i = 1,…n) has a probability mass function, i.e. the observation is zero or non-zero, given by: where π i is the logistic link function with regressors z 1 …z n and 0 ≤ π i ≤ 1, µ i is the mean of negative binomial distribution with regressors x 1 …x n and µ ≥ 0, θ is the dispersion parameter with θ > 0 and Γ(.) is the gamma function.The logistic regression with logit link function is given by: The quantities in negative binomial distribution modelled with regression can be written as: where α is an intercept and β 1 , β 2 , … β k are unknown parameters for regressors (landscape variables) x that are estimated from data.
Zero-inflated models use binomial general linear models (GLM) to model the probabilities of measuring zeros, and the count process is modelled by a Poisson, or as in our case, negative binomial distribution (Zuur et al. 2009).The idea of ZINB modelling is to first estimate the zero part of the model with a binary model, and after that, to calculate estimates for the count part of the model.The probability for a zero to belong to either of the groups is taken into account when modelling count part of the model.As a result, two models are produced.The first that gives the best combination of variables for predicting zero-cases, i.e. non-damage grid cells in our data.The second model gives variables that best predict the number of damage as a function of the variables.
Grid cells that are adjacent or close to each other are likely to have more similar characteristics in terms of moose damage and landscape characteristics than more distant ones and are thus spatially correlated.To account for spatial autocorrelation we used a simple autocovariate model (Dormann et al. 2007) by calculating the mean number of moose damage plantations per year in adjacent grid cells, NNDAMAGE, for each grid cell and forced NNDAMAGE to model as an explanatory variable at each step.
The modelling was made using R and package pscl for zero-inflated negative binomial and binomial modelling and a Vuong test (Zeileis et al. 2008).We used R-package MASS (Venables and Ripley 2002) for calculating a 95% confidence interval for predictions.In calculating confidence intervals, we assumed a normal distribution and used 2500 replicates/random samples for each estimate per unit.Finally, because we only modelled main effects, we also calculated marginal effects of variables and produced graphs of the predicted number of moose damage for the variables included in the final models by using 25%, 50% and 75% quartile levels for the main co-variables in the models.Due to highly skewed distributions, 90 or 95% quantile levels were used for some variables.In this study, we only present model outputs for plantations, because moose damage mostly occurs in these.For zero models we only present models calculated with the mean values of each significant variable estimates.

Modelling strategy
The total of 32 landscape and road variables (Tables 1a,b) that were calculated for each grid cell yield numerous plausible variable combinations and, consequently, too many candidate models to be all tested.In addition, not all the variable combinations are interesting or relevant with respect to the study question.Therefore, we formed ten models (later referred as Model #1, Model#2, ...Model#10) consisting of different combinations of variables (Table 2).Variables that were used in modelling consist of either one or more original landscape variables or the sums of these.
Models were formed in an essentially hierarchical manner, where Model#1 consists of two variables, forest (sum of the proportion of all forests) and non-forest (sum of the proportion of all non-forest areas), and Model#2 of forests and human-modified land (sum of all human modified area).For Model#3, only forest variables were used and forests were divided according to forest age into five classes.For Model#4, young forests were divided according to their age into plantations and thinning forests but also according to the dominance of pine, yielding a total of four variables.For Model#5 pine-dominated young forests were divided according to soil type (4 variables).For model#6, all pine-dominated young forests were divided into peatland and mineral soil forests (2 variables).For Model#7, young peatland forests were divided into drained and nondrained peatland forests (4 variables).In Model#8, all young pine-dominated forests were divided into either drained or non-drained forests (2 variables).In Model#9, the proportions of different types of human-modified non-forest land were used in the modelling (4 variables).In Model#10, we included the length of the different road classes in model (4 variables).
At each step we ran the models as many times as necessary to find only significant variables in both the count and zero models.Finally, we included all the significant variables from Models #1-10 in a final model and calculated estimates for the model as many times as needed to obtain only significant variables in both the zero and count models.We used Akaike´s information criteria (AIC) for model comparisons in each model step and to select the best combination of variables in the models.Models with ∆AIC > 2 were considered to be significantly different (Burnham and Anderson 2002).

Differences among models
Both in Ostrobothnia and Lapland and for all the landscape sizes, several model combinations with significant variables in both count and zero models could be constructed.As judged by ∆AIC, there were, however, large differences among models because the ∆AIC values were more than 10 units, in most cases tens of units and the largest ∆AIC was 276.16 (Suppl.files S3 and S4).Comparison of ZINB and NB models with the Vuong test indicated that ZINB models were better than NB models for all the models.The number of moose damage per year in neighbouring grid cells, NNDAMAGE, was also a significant predictor in all the models constructed and in both the count and zero parts of the models (showing opposite signs of the estimate, though).Therefore, when discussing significant variables in the models below, we focus only on variables other than the number of moose damage in neighboring grid cells.
Table 2. Model steps and the effect of different landscape features on moose damage that were tested in each step.At each step, models were constructed as many times as was necessary to find a combination of only significant variables in both the count and zero models.A final model was constructed by entering all significant variables from each step to the same model and by finding variable combinations with the lowest AIC and where all variables in the both count and zero models were significant.At each step, also the average number of damage in neighboring cells (NNDAMAGE) was forced to both the count and zero models.

Ostrobothnia 1 km 2
In the final count model for Ostrobothnia 1 km 2 landscapes, pine-dominated plantations on mineral soil, non-pine-dominated plantations, thinning forests, mature forests and inhabited areas were the best predictors (Table 3).In the zero model increasing proportion of mature forests reduced the probability for non-damage cells (Table 3, Suppl.file S5).
The predicted number of damage exceeded the average number of damage in the actual data already at very low proportions of both mature forests and thinning forests (data not shown).A similar effect was found for pine-dominated plantations on mineral soil when the proportion of nonpine-dominated plantations was used as a covariate (data not shown).When different proportions of non-pine-dominated plantations were used as covariates their effects on the predicted number of damage as a function of pine-dominated plantations were rather mild.Instead, the negative effect of inhabited areas on the number of damage predicted with non-pine-dominated plantations is clear and the predicted number of damage decreases by factor 0.125 when the proportion of inhabited areas increases from 1.6% to 7.1%, i.e., from one damage per 178 ha to one damage per 1428 ha (Fig. 2).For pine-dominated plantations on mineral soil, the decrease in the number of damage was by a factor of 0.166, i.e., from one damage per 178 ha and to one damage per 1071 ha.However, the model overestimates the number of damage for Q25 and Q50, therefore, the data is not shown.

Ostrobothnia 25 km 2
The proportion of all plantations, mature forests and inhabited areas were the best predictors for the count model and for the Ostrobothnia 25 km 2 landscapes (Table 3).The probability for zero damage cells decreased with increasing proportions of thinning and mature forests (Table 3, Suppl.file S5).
The predicted number of damage as a function of the sum of all plantations increased by factor 1.6 when the proportion of mature forests increased from 30.5% to 47.9% (Fig. 3).In practice, this means that where there is an average of one damage per 195 ha of plantations with 30.5% mature forests the model predicts one damage per 120 ha plantations when the proportion of mature forests is 47.9%.Instead, an increase in inhabited areas decreased the number of damage by factor 0.8 when the proportion of inhabited areas increased from 0.2% to 5%, i.e. the model predicted one damage per 153 ha and 186 ha plantations, respectively (Fig. 3).

Ostrobothnia 100 km 2
The proportion of all plantations and mature forests were the best predictors for the count model for 100 km 2 landscapes in Ostrobothnia (Table 3).For the zero model, mature forests and open land in mineral soils were the best predictors and an increase in both variables decreased the probability of finding a non-damage grid cell (Suppl.file S5).However, due to very large confidence intervals the proportion of open land on mineral soils did not yield feasible predictions (data not shown).The number of predicted moose damage increased by factor 1.4 when the proportion of mature forests increased from 31.8% to 45.3% (Fig. 3), i.e. one damage per 211 ha and 151 ha of all plantations were predicted, respectively.

Lapland 1 km 2
Pine plantations on mineral soil, non-pine-dominated plantations and the sum of thinning forests were the best predictors of the number of damage for 1 km 2 landscapes in Lapland (Table 4).The increasing length of connecting roads also showed a negative effect on the number of moose damage.For the zero model, open areas on minerals soils, non-pine-dominated plantations, the sum of thinning forests and mature forests predicted that the probability of finding a non-damage grid cell decreases as the proportions of these variables increase (Table 4, Suppl.file S6).Marginal effect plots of the count models showed that for both the proportion of pine plantations on mineral soils and for the proportion of non-pine-dominated plantations on all soils the models predict two to three times more damage than shown on average in the actual data (data not shown).The effect of connecting roads had a negative effect on the predicted number of damage (Fig. 2) but, as with other covariates, the model predictions were about two-fold compared with actual data.
Zero models either predicted a non-damage grid cell to be found with very high proportions of variables or they had confidence levels up to 60-70% indicating that the model performance was not very good for any of the variables (Suppl.file S6).

Lapland 25 km 2
The sum of all plantations together, thinning forests on minerals soil, and mature forests best explained the number of moose damage in 25 km 2 landscapes in Lapland (Table 4).The occurrence of non-damage grid cells was best explained by the proportion of pine plantations on mineral soils and non-pine-dominated thinning forests (Suppl.file S6).
According to the mean model, the predicted number of damage as a function of the proportion of all plantations showed a rather slow increase until the proportion reached about 40-50 % of the 25 km 2 landscape (Fig. 4).However, the 95% confidence intervals start to grow substantially at the same point.The predicted number of damage increased along with increasing proportion of mature forests (Fig. 4).When the proportion of mature forests increased from 12.8% to 25.4% the number of damage was predicted to increase by factor 1.5, i.e. from one damage per 630 ha to one damage per 407 ha, respectively.An increase in pine-dominated thinning forests, instead, reduced the predicted number of damage by factor 0.86, i.e. from one damage per 463 ha plantations to one damage per 537 ha plantations.Similarly, an increase in connecting roads from 1.9 km to 8.1 km reduced the predicted number of damage by factor 0.83, i.e. from one damage per 463 ha to one damage per 537 ha (Fig. 4).

Lapland 100 km 2
The sum of all plantations and mature forests were the only significant predictors for 100 km 2 landscape count model in Lapland (Table 4).The number of moose damage increased along with increasing percentage of both variables.The sum of all plantations and the length of forest roads predicted that the probability of non-damage grid cells decreases when the proportions of these variables increase.However, for both variables 95% confidence levels were large, up to ±73% of the mean, indicating poor model performance (Suppl.file S6).
The mean model showed that the predicted number of damage increases slowly as a function of all plantations until the proportion of plantations is about 40-50 % of the landscape (Fig. 4).At the same time, however, the 95% confidence levels of the predicted damage also increase strongly.An increase in mature forests from 14.3% to 25.5% increased the predicted number of damage per plantation areal unit by a factor of 1.4, i.e. from one damage per 191 ha to one damage per 134 ha.

Mature forests
According to our models "Mature forest" was a significant variable in five out of six final count models and in four out of six zero models.An increase in the proportion of mature forests increased the number of damage in all models and decreased the probability of no moose damage in a grid cell.Our results are in accordance with Olsson et al. (2011) and Bjørneraas et al. (2011), who found that moose home ranges included more mature forests than generally found in landscapes.In winter, mature forests do not provide much forage for moose, but they might provide better cover and, perhaps, lower energy costs due to lower snow depth or the structure of the snow (Bjørneraas et al. 2011).
The effect of mature forests on moose browsing risk might also be related to landscape structure in terms of grain size (Cassing et al. 2006).The average size of regeneration area in privately owned forests in Finland is 1-2 hectares (source Finnish Forestry Centre).Over time, regeneration areas become embedded in the mosaic of forests of different successional stages.Due to small overall compartment size, the distance between cover (mature forest) and food resource (plantations) is rather short, especially when related to the home range size of moose.In Fennoscandia there is no evidence of the effect of patch size on the intensity of browsing in different parts of plantations (Andren and Angelstam 1993), which might imply that due to increased cover-food adjacencies moose might perceive forestry modified landscapes as fine grained.In fine grainedlandscapes moose browsing is more evenly distributed than in landscape with large variation in patch size (Edenius et al. 2002).
In our study areas mature forest is the most common habitat type out of all original habitat types studied, at just under 40% in Ostrobothnia and about 20% in Lapland, and thus the probability of finding mature forest in a grid cell is higher than all other habitat types.It should be remembered, however, that according to our habitat classification mature forests included all forests that had >81 m 3 ha -1 of total volume of trees (Table 1) and thus our forest class "Mature forests" category includes not only old-growth forests, sensu stricto, but also advanced thinning forests.According to definition (Metsäntutkimuslaitos 2002) advanced thinning forests comprise some timber-sized trees and, in our study areas, these forests were 40-60 years old.

Thinning forests
The effect of thinning forests on the number of damage varied according to scale and type of thinning forests.In the 1 km 2 models, an increase in all thinning forests predicted the number of damage to increase in both Ostrobothnia and Lapland.The mechanisms behind this are probably the same as for mature forests; thinning forests do not provide much food in winter but they might provide better cover and lower energy costs due to lower snow depth or the structure of the snow (Bjørneraas et al. 2011).Furthermore, the food-cover adjacencies might have the same effect in thinning forest plantation mosaics as mature forests.
In the 25 km 2 model in Lapland, the number of damage decreased when the proportion of pine-dominated thinning forests increased, although a post hoc analysis showed a significant positive correlation between the proportion of pine-dominated thinning forests and pine-dominated plantations (Spearman's rho = 0.685, p < 0.001).One plausible explanation would be the legacy of post WWII forestry in Lapland, when large forest areas were clear cut and regenerated mainly with Scots pine.Deciduous trees are a common admixture in plantations regenerated for pine and cleaning and soil scarification are used to reduce competing vegetation in plantations.In the least productive sites plantations also develop more or less as pine-dominated monocultures throughout the rotation period.Dense pine-dominated thinning forests thus provide only a limited amount of browse for moose in winter.Therefore, moose probably do not favour areas dominated by pinedominated thinning forests but establish their winter ranges in areas with a more diverse mosaic of forests comprising other species in addition to pine.

Inhabited areas and roads
According to our 1 km 2 and 25 km 2 Ostrobothnia models, the number of moose damage decreased as the amount of inhabited areas increased.The effect was clear especially when the number of damage was predicted as a function of pine-dominated plantations on mineral soils and inhabited areas were used as a covariate.In general, the average proportion of inhabited areas in our data was rather low, from 1.6% to 2.1%, depending on the grid cell size (Suppl.file S1).However, the proportion of inhabited areas was on average two to six times higher in no-damage grid cells than in damage grid cells, the highest difference being in 1 km 2 grid cells.The difference between damage grid cells and the overall landscape was marginally smaller, but still about the same magnitude (Suppl.file S1).There were similar differences in Lapland, although inhabited area was not a significant variable in any of the models.
Although moose can exploit variations in food and cover caused by human land use (Bjørneraas et al. 2011), the vicinity of inhabited areas and human activities causes a trade-off in habitat use (Lykkja et al. 2009).In addition to avoidance of areas close to human habitation (Tinoco Torres et al. 2011), the time during which moose are most active, e.g.browsing, might also be shorter in areas close to human habitation as compared with more remote areas (Lykkja et al. 2009).In this sense, the response of moose to human disturbance has been suggested to be analogous to the response to predator risk (Frid and Dill 2002), according to which human disturbance can induce predator avoidance behaviour and cause animals to shift their habitats farther from the disturbance source.According to this hypothesis, one pre-requisition for habitat shift is that there are available habitats for foraging also farther from inhabited areas.Since WWII, forestry has strongly changed the structure of the forest landscape in Finland in terms of the amount and pattern of different age classes, and has especially increased the amount of young successional stages.Therefore, moving farther from inhabited areas is unlikely to lower the probability of moose finding good quality habitats also in there.
Our finding regarding the effect of roads on damage risk is in line with some earlier studies that suggest that the vicinity of roads decreases damage risk (Repo and Löyttyniemi 1985).However, Ball and Dahlgren (2002) found that major highways might increase moose damage close to roads by inhibiting moose from migrating across roads and by increasing moose density near to roads.According to our findings, the effect of roads was the opposite; the higher the length of connecting roads, the less damage incurred in 1 km 2 and 25 km 2 landscapes in Lapland.One explanation for the discrepancy between these findings might be that roads in our study area do not generally form barriers for moose to cross.Ball and Dahlgren (2002) did not attribute the barrier effect to the presence of roadside fences, but rather to the open area formed by major highways.In Lapland, connecting roads, the most common road type after forest roads, are 4-5 m wide and generally asphalted.Moose vehicle collisions occur more frequently on connecting roads than any other class of road in Lapland (Finnish Transport Agency 2014), which indicates that connecting roads are not barriers for moose.However, according to our results connecting roads do seem to have an inhibitive effect on moose damage at a distance of up to 2.5-3.5 km from the road.This, again, could be explained by the predator risk hypothesis (Frid and Dill 2002).

Plantations
The number of damage as a function of the proportion of plantations increased quite slowly until the proportion of plantations was about 18% in Ostrobothnia and about 32% in Lapland, which correspond to average landscapes in both areas, respectively.However, at the same time, the uncertainty of the predictions starts to increase strongly, probably due to fewer observations with a high proportion of plantations.The amount of any type of plantations or the total amount of plantations were not good predictors alone for the number of damage.Instead, the number of damage increases when the landscape forms a mosaic of plantations and mature forests or thinning forests.This is in line with (Bjørneraas et al. 2012;Herfindal et al. 2015) who found moose to select areas with habitat types that provide both good forage and cover.

Transferring models between regions
Our results partly support the idea that the transfer of habitat selection models from one area to another should be made with care (Mabille et al. 2012).Our study areas, Lapland and Ostrobothnia differed from each other in several aspects.In Ostrobothnia, there were more drained peatlands, mature forests, inhabited areas and agricultural land than in Lapland.Instead, there were more regeneration areas, plantations and young thinning forests in Lapland.The road network is also denser and there are more inhabited areas in Ostrobothnia than in Lapland.However, the variables in the models for respective landscape sizes were quite similar for the two areas, indicating that moose damage risk is related to similar habitat composition in both areas.Some differences, like the effect of inhabited areas, roads and pine-dominated thinning forests, still imply that moose damage risk models should include landscape features typical to each region.

Count models and zero models
There is empirical evidence that the abundance of moose can to some degree be predicted on the basis of predicted occurrence and by partly using the same environmental covariates used in both occurrence and abundance models (Nielsen et al. 2005).If this holds true also for moose damage, data of moose damage and the location of damage in relation to available habitats might be used in predicting the number of damage on the basis of damage probability.However, Michaud et al. (2014) found that the same indicators were ineffective in concurrently predicting both the occurrence and abundance of moose.According to our models, only a few variables explained both the number of damage and the occurrence of damage.It should be remembered that zero models describe conditions where moose damage was not recorded, not conditions where damage was found.Our data results from several mechanisms that can lead to non-damage grid cells and some non-damage grid cells can be so-called false zeros (Martin et al. 2005), i.e. moose damage are present but not recorded.The simplest explanations for this are that existing damage have not been observed, or that damage has been observed but no damage compensation has been sought by the land-owners.Due to the large number of compensated damage spanning several years, we believe, however, that our data are sufficient to describe at least the variation in relative occurrence of damage in our study area.

General model performance and evaluation
As our modelling shows, it is possible to construct numerous models with reasonable variables and with several variable combinations (Suppl.files S3 and S4).However, when arranged by their AIC, some models are tens or even hundreds of times more informative than others.We do not claim that those models with variables that were not included in the final models are necessarily not useful or uninformative, but instead denote the view of Burnham and Anderson (2002) who regarded model construction as a process of testing several options rather than building the most parsimonious model.Our results do, however, urge careful consideration of the results in relation to the type of landscape/habitat data used in modelling.
Habitat models might perform poorly if habitat selection for a species is strongly density dependent and there is no information available about the population size in relation to long-term cycles or carrying capacity (Boyce et al. 2016).For moose in Finland, information on estimated population size after hunting is available only at the game management area level at scales of 10 2 -10 4 km 2 (Nygrén and Pesonen 1993), while estimates for smaller areas are missing.In winter, moose tend to aggregate to certain winter pastures which means that population densities can be very high and damage can be the result of several individuals.Although studies in Fennoscandia have no found correlation between population density and damage level (Hörnberg 2001b) or density-dependent habitat selection (Månsson 2009), high densities of moose can locally lead to the utilization of sub-optimal habitats.This, again, might level out the effect of varying habitat structure on the amount of browsing.
Finally, our study was based on only one measure of landscape structure: composition of habitats.There is evidence that other indices of landscape structure, like such as patch size, are also linked to browsing intensity (Cassing et al. 2006).Therefore, using other measures of landscape structure, such as mean patch size, distance to other similar patches, patch density and length of edges, as explanatory variables in modelling might provide further insight into moose damage risk in a mosaic landscape (Cassing et al. 2006).Additional variables not measured in this study might also improve the performance of the models at the smallest scale.Furthermore, resources at landscape and region levels might vary due to biogeographical variation, topography or land use (Michaud et al. 2014;Nikula 2017), which should be assessed when evaluating moose damage risk.

Fig. 1 .
Fig. 1.An illustration of the study setup.Grid cells of 1 km, 5 km and 10 km were overlaid on located moose damage and land use and forest classification (for a complete list of variables see Tables1a,b).Grey areas on the index map show the location of privately owned land in the study area and bolded lines delineate the provinces of Ostrobothnia and Lapland.

Fig. 2 .
Fig. 2. Predicted number of moose damage per 1 km 2 landscape size in relation to the proportion of non-pine-dominated plantations and different proportions of inhabited areas in Ostrobothnia (upper row) and the length of connecting roads in Lapland (lower row).Predictions were calculated for different proportions of non-pine-dominated plantations by using Q1 (left column), Q2 (middle column) and 90% quartile (right column) for the proportion of inhabited area and length of roads while keeping all other possible covariates at their means.Arrow lines show the average number of damaged plantations during years 2002-2008 (Y-axis) and change in moose damage risk per plantation area (X-axis) as a function of different levels of co-variates.Grey area shows 95% confidence limits.

Fig. 3 .
Fig. 3. Predicted number of moose damage per 100 km 2 and 25 km 2 landscape sizes in Ostrobothnia.Predictions were calculated for different proportions of plantations using quartiles of Q1 (left column), Q2 (middle column) and Q3 (right column) for the proportion of mature forests and inhabited areas while keeping all other possible covariates at their means.However, due to the highly skewed distribution 95% was used for inhabited areas as the highest value.Arrow lines show the average number of damaged plantations in Ostrobothnia during years 2002-2008 (Y-axis) and the change in moose damage risk per plantation area (X-axis) as a function of different levels of co-variates.Grey area shows 95% confidence limits.

Silva
Fennica vol.53  no. 1 article id 9918 • Nikula et al. • Modelling the effect of habitat composition and roads…

Fig. 4 .
Fig. 4. Predicted number of moose damage per 100 km 2 and 25 km 2 landscape sizes in Lapland.Predictions were calculated for different proportions of plantations using quartiles of Q1 (left column), Q2 (middle column) and Q3 (right column) for the proportions of mature forests and thinning forests and for the length of connecting roads while keeping all other possible covariates at their means.Arrow lines show the average number of damaged plantations in Lapland during years 2002-2008 (Y-axis) and the change in moose damage risk per plantation area (Xaxis) as a function of different levels of co-variates.Grey area shows 95% confidence limits.

Table 1a .
(McGarigal and Marks 1995) classification was formed of 19 classes whose criteria are shown in right column.Middle column gives a description of habitat type, land use and cover type and road class.The total number of damaged plantations used in modelling varies among cell sizes due to the criteria used in selecting cells for modelling.For each grid cell we calculated the proportion of each land use and forest class (Table1a) using Fragstats(McGarigal and Marks 1995)and the length of each road class with ArcMap 10.1 (Environmental Systems Research…

Table 1b .
In the different steps of modelling further variables #20-32 were formed by summing the original classes (right column).Middle column gives a description of habitat type, land use and cover type and road class.The average number of damage per year in neighboring cells (NNDAMAGE) was forced to both the count and zero models at each step and to the final models.

Table 3 .
Coefficients, standard error of mean, Z-test value, and significance of the best combination of model variables explaining the amount of moose damage (count model) and the occurrence of a non-damage landscape (zero model) in Ostrobothnia for 1 km 2 , 25 km 2 and 100 km 2 landscape sizes.See Tables 1a,b for variable descriptions.

Table 4 .
Coefficients, standard error of mean, Z-test value, and significance of the best combination of model variables explaining the amount of moose damage (count model) and the occurrence of non-damage landscape (zero model) in Lapland for 1 km 2 , 25 km 2 and 100 km 2 landscape sizes.See Tables 1a,b for variable descriptions.