Skip to main content
  • Research article
  • Open access
  • Published:

A spatially-explicit count data regression for modeling the density of forest cockchafer (Melolontha hippocastani) larvae in the Hessian Ried (Germany)

Abstract

Background

In this paper, a regression model for predicting the spatial distribution of forest cockchafer larvae in the Hessian Ried region (Germany) is presented. The forest cockchafer, a native biotic pest, is a major cause of damage in forests in this region particularly during the regeneration phase. The model developed in this study is based on a systematic sample inventory of forest cockchafer larvae by excavation across the Hessian Ried. These forest cockchafer larvae data were characterized by excess zeros and overdispersion.

Methods

Using specific generalized additive regression models, different discrete distributions, including the Poisson, negative binomial and zero-inflated Poisson distributions, were compared. The methodology employed allowed the simultaneous estimation of non-linear model effects of causal covariates and, to account for spatial autocorrelation, of a 2-dimensional spatial trend function. In the validation of the models, both the Akaike information criterion (AIC) and more detailed graphical procedures based on randomized quantile residuals were used.

Results

The negative binomial distribution was superior to the Poisson and the zero-inflated Poisson distributions, providing a near perfect fit to the data, which was proven in an extensive validation process. The causal predictors found to affect the density of larvae significantly were distance to water table and percentage of pure clay layer in the soil to a depth of 1 m. Model predictions showed that larva density increased with an increase in distance to the water table up to almost 4 m, after which it remained constant, and with a reduction in the percentage of pure clay layer. However this latter correlation was weak and requires further investigation. The 2-dimensional trend function indicated a strong spatial effect, and thus explained by far the highest proportion of variation in larva density.

Conclusions

As such the model can be used to support forest practitioners in their decision making for regeneration and forest protection planning in the Hessian Ried. However, the application of the model for predicting future spatial patterns of the larva density is still somewhat limited because the causal effects are comparatively weak.

1Background

The forests in southern Hessen in the vicinity of the Rhine-Main urban agglomeration are among the most problematic forest areas for forest management in Central Europe. Here extraordinary demands are made on forests and forest enterprises in view of the high population density, highly concentrated industrialization and dense road and traffic infrastructure. Urban agglomeration has led to the acquisition of considerable land area for development, which, in turn, has led to an unusually high fragmentation of the forest area. Furthermore environmental impacts, including an increase in emissions of air pollutants and the serious lowering of the groundwater table arising from high water usage, are evident. Forestry management options are severely restricted in this region due to the exceptional importance of forests, particularly old oak forests, for recreation and nature conservation (NW-FVA, Nordwestdeutsche Forstliche Versuchsanstalt Hrsg. [2013]): p. 32 ff.). The abiotic pressures and especially the wide-spread drop in groundwater levels have already degraded many forests. Additionally, the massive outbreaks of biotic pests like the forest cockchafer (Melolontha hippocastani), the European oak leafroller (Tortrix viridana) and the gypsy moth (Lymantria dispar) have in part led to a total destruction of forests in some areas. It is anticipated that site conditions will degrade further in the future if the groundwater depletion intensifies and if climate change leads to higher temperatures and lower precipitation as projected (NW-FVA, Nordwestdeutsche Forstliche Versuchsanstalt Hrsg. [2013]): p. 40 ff.). Already the Rhine-Main region is one of driest and warmest regions in Hessen.

The Hessian Ried is located in the southern part of the Rhine-Main region, and encompasses an area of about 100 000 hectares, which includes about 30 000 hectares of forest area (Figure 1A). Although the Hessian Ried is a less densely populated part of the Rhine-Main region, impacts from the surrounding economic and industrial centers are evident. Within the framework of a larger research project, changes in the forests in the Hessian Ried have been projected assuming different groundwater management (NW-FVA, Nordwestdeutsche Forstliche Versuchsanstalt Hrsg. [2013]): p. 37 ff. and climate change scenarios (NW-FVA, Nordwestdeutsche Forstliche Versuchsanstalt Hrsg. [2013]): p. 40 ff.). In this context standard methods for growth and yield simulation were used. Additionally, models for quantifying the effects of biotic and abiotic risks were developed to generate more realistic prognoses.

Figure 1
figure 1

Geographical location (Gauß-Krüger projection) of the Hessian Ried and layer of the forest cover (A), inventory grid of sample excavation plots for recording larva population data in the Hessian Ried (B).

One of the most serious regional biotic risks, especially in the forest regeneration phase, is the forest cockchafer (Melolontha hippocastani). The northern Upper Rhine region is populated by different tribes of the forest cockchafer each of which has a specific quadrennial life cycle and swarming year. The forest cockchafer population in the Hessian Ried belongs to the South Hessian tribe. Since the early 1980s a massive outbreak of this tribe has been observed in some parts of the Hessian Ried. The area of forest populated by the forest cockchafer has increased constantly and this trend is expected to continue. An area covering 13 000 hectares of the total forest area of 30 000 hectares is assumed to be populated by this species. Of this area, 4 000 hectares are assumed to be infested by extremely high larva densities.

Generally the damage caused by the adult beetle browsing the leaves, which occurs at the end of each quadrennial live cycle, is not critical to forest health (Schwerdtfeger [1981]). In contrast, damage caused by the larvae eating the roots is most serious. Natural regeneration and the planting of seedlings fail repeatedly: middle-aged and old stands are increasingly at risk. In some parts of the Hessian Ried, the current high larva densities prohibit the implementation of any regeneration measures in the absence of prior implementation of forest protection measures such as the application of chemical or biological pesticides (Ott et al. [2006]). Therefore, given the outstanding recreation and natural conservation functions of the Hessian Ried, in particular, a detailed risk analysis is essential. Consequently, the density and dispersal behavior of the larva population over time need to be assessed. For this purpose a systematic sample plot inventory of larva density was initiated in 2009, in which the population of larvae in excavated pits was recorded for each plot. Based on this excavation inventory a specific regression model for count data was developed that permitted an area-wide prediction of the larva densities taking into account significant causal covariate effects and the spatial autocorrelation of the data. Thus two major aims of the model development were to:

  1. 1)

    predict larva density by area across the Hessian Ried for decision support;

  2. 2)

    identify significant causal variables and quantify their effects to enhance the generality of the model and to gain greater insight into the suitability of a site to serve as a habitat for larvae.

Moreover, the development of a generalized model approach that could be extended for the investigation of spatial population densities of the forest cockchafer larvae in the future as time series data become available was envisaged. The approach would ensure the future relevance of the model to forest managers in the planning and implementation of optimal spatial forest regeneration and pest control measures and would allow for the investigation of the change in the spatial pattern of larvae density over time.

2Methods

2.1 Data

In the summer of 2009 counts of larvae and additional potential covariate data were sampled over a regular quadratic grid with the sample plots arranged 500 m apart. In a small area of 994 hectares, a denser sample grid of 250 m × 250 m was employed to obtain additional data needed for further investigations in forest protection measures (Figure 1B). Each sample plot on the larger and smaller scale grid system consisted of 4 subplots each 50 cm × 50 cm (0.25 m2) in size, which were excavated to count the number of larvae present. These subplots were located 10 m away from the plot center in the 4 cardinal directions.

The excavation of a subplot was only conducted if it was located in a forested area. In few cases one or more of the subplots in a sample plot were ignored because they were located in a bog, in settlement areas etc. In total, data from 1 276 sample plots could be used for constructing the model. The excavated soil was carefully searched for forest cockchafer individuals at all relevant life stages that were expected to swarm in 2010. Hence, not only the number of 3-year-old larvae, but also the number of pupae and imagines were recorded.

For modeling the sum of larvae in the 4 × 0.25 m2 subplots at each sample point, that is the larva density per square meter, was used. The subplots were treated as one sample plot for larva excavation 1 m2 in area since the distances between the 4 subplots were small compared to the distances between the sample points on the sample grid. The maximum depth of excavations was 1 m. However, at several subplots, the mechanical resistance of the soil prevented excavations reaching this depth. In this case, it was assumed that, if the mechanical resistance was too high for excavating, it was also too high to serve as a habitat for larvae.

Geographic location (Gauß-Krüger projection), clay thickness to 1 m soil depth (CTH), distance to water table (DWT) (Table 1), available from other data sources, were used as predictor variables. Data for forest stand structure were not recorded or available. CTH was defined as the percentage of pure clay layer up to 1 m soil depth. CTH was calculated using a regional area-wide soil substrate mapping (NW-FVA, Nordwestdeutsche Forstliche Versuchsanstalt Hrsg. [2013]): p. 25 ff., and DWT was derived from a regional simulation of the groundwater level in the Hessian Ried (NW-FVA, Nordwestdeutsche Forstliche Versuchsanstalt Hrsg. [2013]): p. 37 ff.). However, based on the soil substrate map, a clay layer was present in 54 sample plots only (Table 1). DWT varied considerably (Table 1). DWT was predicted to be less than 4 m at 460 (37%) sample plots, less than 3 m at 267 (21%) sample plots and less than 2 m at 140 (11%) sample plots. The distribution in larva densities was extremely skewed to the right, with 52% of sample plot excavations comprising no larvae. Densities higher than 5 larvae · m−2, and higher than 10 larvae · m−2 were recorded at 20.7%, and 11.0% of the sample plots. A maximum larva density of 56 larvae · m−2 was recorded.

Table 1 Distribution of larva density at 1 276 sample plots recorded in the larva density inventory and assigned covariates distance to water table ( DWT ) and clay thickness ( CTH)

2.2 Methodology

2.2.1 Regression models for count data

For the estimation of larva density per square meter, different generalized additive regression models (gam) were parameterized and compared during model selection. The larval densities per square meter are count data that can have only non-negative integer values. Hence, discrete distribution functions that take into account the specific characteristics of count data were selected for modeling. In exceptional cases categorical regression models (Yee and Wild [1996], Yee [2010]) offer an alternative if only few different low counts are observed (Fahrmeir et al. [2007]). Hence, in this investigation different methods had to be applied because of the rather wide distribution of larval counts (Table 1). The Poisson distribution is often used for modeling count data (McCullagh and Nelder [1989]): 194 ff. The natural logarithm is used as a link function, giving a general model defined by:

g λ i = η i =x ' i β= β 0 + β 1 x i1 ++ β k x ik
(1)

or

λ i =exp η i =exp x i β =exp β 0 exp β 1 x i 1 exp β k x ik

with g(.): link-function: natural logarithm

and y i | x i ~ Poisson λ i with E y i | x i = λ i and Var y i | x i = λ i ; y i = 0 , 1 , 2 , .

the probability density function of the Poisson distribution is given by:

f y ; λ = λ y e λ y ! y = 0 , 1 , 2 , .

In the Poisson distribution, the conditional mean equals the conditional variance and is determined by the multiplication of the exponentials of the predictor variables. By applying the inverse link-function, positive predictions result from the model. However in many applications, the empirical variance is found to be higher than the value assumed by the Poisson distribution. In this case, an overdispersion parameter is introduced, and the conditional variance is defined as follows: Var (y i |x i ) = ϕλ i with the unknown constant overdispersion parameter ϕ.

Alternatively overdispersion in count data can be dealt with by using different probability functions such as the negative binomial distribution (McCullagh and Nelder [1989]: 199; Venables and Ripley [2002]: 206 ff; Rigby and Stasinopoulos [2009]: 269). In this case the overdispersion is modeled explicitly. However, the procedures for the estimation are more extensive, particularly if the dispersion parameter ϕ is unknown and must be estimated. The negative binomial distribution can be derived from a two-stage model for the distribution of a discrete variable Y (Venables and Ripley [2002]). The variable Y follows a Poisson distribution, where the mean is regarded as a random variable itself that has a gamma distribution, with shape parameter ϕ , and scale parameter 1/ϕ, or alternatively with mean = 1 and variance 1/ϕ. ϕ equals the dispersion parameter of the negative binomial distribution:

g 1 μ i = η 1 i =x ' i β 1 = β 01 + β 11 x i 1 ++ β k 1 x ik
(2a)

with g 1(.): link-function: natural logarithm

and y i | x i ~ Poisson μ i with μ i ~ Gamma shape = ϕ , scale = 1 / ϕ = Gamma mean = 1 , variance = 1 / ϕ and y i | x i ~ NegBin μ i , ϕ with E y i | x i = μ i and Var y i | x i = μ i + μ i 2 / ϕ ; y i = 0 , 1 , 2 , .

the probability density function of the negative binomial distribution is given by:

f y , θ , μ = Γ ϕ + y Γ ϕ y ! μ y ϕ ϕ μ + ϕ ϕ + y

, with Γ denoting the gamma function. y = 0, 1, 2, ….

An additional flexibility can be achieved if not only the mean μ i but also the dispersion parameter ϕ i is estimated as a function of the predictor variables by employing a second linear predictor (Eq. 2b). Instead of estimating the dispersion parameter ϕ i directly, usually the scale parameter or the variance 1/ϕ i of the gamma distribution is estimated (Rigby and Stasinopoulos [2009]: 269):

g 2 1 / ϕ i = η 2 i =x ' i β 2 = β 02 + β 12 x i 1 ++ β k 2 x ik
(2b)

with g 2(.): link-function: natural logarithm

Since the larva density database revealed a high percentage (52%) of sample plots with no larvae present (Table 1), as a third alternative, a regression model assuming a zero-inflated Poisson distribution was fitted and validated (Rigby and Stasinopoulos [2009]: 270). If Y = 0 with probability ω, and Y ~ Poisson(λ) with probability (1 − ω), then Y has a zero-inflated Poisson distribution. Both parameters λ i and ω i are estimated via specific linear predictors (Eqs. 3a and 3b):

g 1 λ i = η 1 i =x ' i β 1 = β 01 + β 11 x i 1 ++ β k 1 x ik
(3a)

with g 1(.): link-function: natural logarithm

g 2 ω i = η 2 i =x ' i β 2 = β 02 + β 12 x i 1 ++ β k 2 x ik
(3b)

and g 2(.): link-function: logistic link function

and

y i | x i ~ ZIP λ i , ω i with E y i | x i = 1 ω i λ i and Var y i | x i = λ i 1 ω i ( 1 + λ i ω i ) y i = 0 , 1 , 2 , .

The probability density function of the zero-inflated Poisson distribution is given by:

f y ; ω , λ = ω + 1 ω e l , y = 0 1 ω λ y y ! e y , y = 1 , 2 , .

2.2.2 Generalized additive models

Generalized additive regression models gam (Hastie and Tibshirani [1990]; Wood and Augustin [2002]) were employed to test the effects of predictor variables for potential non-linearity. By way of example, the general definition of a gam that contains a spatial trend function is given below for the Poisson regression model (Eq. 4):

g λ i = η i = β 0 + f 1 x i 1 ++ f k x ik + f n eas t i , nort h i
(4)

with g(.): link-function: natural logarithm

and y i | x i ~ Poisson λ i and E y i | x i = λ i and Var y i | x i = λ i ; y i = 0 , 1 , 2 , .

f 1, f 2, f 3,…, f k : 1-dimensional smooth functions

f n : 2-dimensional smooth function for modeling a spatial trend.

Assuming negative binomial and zero-inflated Poisson distributions, both linear predictors (Eqs. 2a/2b and 3a/3b) were checked for non-linear model effects. Subsequently it was proofed if significant non-linear effects could be adequately approximated by segmented linear effects. All additive regression models were fitted with a 2-dimensional smoothing function of the geographic location to test the data for spatial autocorrelation, which could not be described by observed causal predictors. To parameterize the model, the statistical language and environment R (R Development Core Team [2010]) was adopted using the two libraries gamlss (Rigby and Stasinopoulos [2005]) and mgcv (Wood [2006]). If possible the models were fitted using the library mgcv only. If the distribution assumptions or specific linear predictors could not be specified with the functions of the library mgcv, functions of both libraries were combined: functions of gamlss were used to specify the distribution functions for the response variable, and functions of mgcv were used to apply specific smoothing techniques, including the 2-dimensional smooth function. The library gamlss was adopted because it enables a variety of continuous and some discrete distributions to be modeled and extends the classical generalized linear (McCullagh and Nelder [1989]) and additive (Hastie and Tibshirani [1990]) models, which are otherwise limited to the distributions of the exponential family (Fahrmeir et al. [2007]): 218.

2.2.3 Validation

In the validation procedure randomized quantile residuals were calculated (Dunn and Smyth [1996]) to determine which distribution assumption and model specification described the pattern of larval counts best. Randomized quantile residuals have been found to be superior in the validation of regression models for continuous and discrete response variables to Pearson and deviance residuals (Dunn and Smyth [1996]).

Let y 1,…, y n be the responses that are independently distributed following a probability density distribution f(y i , μ i , ϕ) where μ i = E(y i ) and ϕ is a parameter vector common to all y i . F(y i , μ i , ϕ) denotes the corresponding probability distribution function, Φ denotes the distribution function of the standard normal distribution, and Φ−1 its inverse denotes the corresponding quantile function. For continuous response variables, the quantile residual is then easily defined by:

r i = Φ 1 F y i ; μ ^ i , ϕ ^
(5)

Given the hypothesis that f(y i , μ i , ϕ) is the correct model for the observations y 1,…,y n , all r i follow approximately a standard normal distribution. Deviations only result from sampling variability in μ ^ i and φ ^ . Hence standard methods such as the ‘normal quantile-quantile plot’ qqplot or the Kolmogorov-Smirnov test can be used to validate the regression models. An application of this approach in forestry is given by Zucchini et al. ([2001]).

If F is not continuous, however, a more general definition of quantile residuals is required leading to the definition of randomized quantile residuals:

r i = Φ 1 u i
(6)

where u i is a uniform random variable in the interval (a i , b i ] with a i =li m y yi F y i , μ ^ i , ϕ ^ and b i =F y i , μ ^ i , ϕ ^ .

Randomized quantile residuals are also approximately standard normal distributed if the model assumed is correct. Further insight into the adequacy of the models can be gained by wormplots that are detrended ‘normal quantile-quantile plots’ (qqplots) (Buuren and Fredriks [2001]). Here ‘detrended’ means that the empirical quantiles are subtracted from their corresponding standard normal quantiles. In a wormplot, these detrended quantiles are plotted against their corresponding (theoretical) standard normal quantiles. Wormplots highlight deviations from an assumed theoretical distribution more clearly than qqplots. Additionally, conditional wormplots were calculated for different intervals of DWT. Conditional wormplots are employed to detect potential intervals of predictor variables where the models do not fit adequately.

Furthermore the assumption of spatially independent model residuals and the methodology of 2-dimensional surface fitting for modeling spatially correlated data were validated. Therefore simple empirical semi-variograms of the randomized quantile residuals were calculated for various models with spatial trend functions of different complexity. To construct 95% confidence intervals for the null hypothesis of spatial independency, 1 000 random permutations of the coordinates were conducted. Subsequently the semi-variograms of the randomized quantile residuals were estimated and the pointwise 2.5% and 97.5% quantiles calculated.

Finally model computations were conducted for specific settings of predictor variables to illustrate the sensitivity of model predictions to varying conditions.

3Results

3.1 Model selection

The specific Poisson regression model selected for predicting larva density per m2 using all available predictor variables was defined as follows:

g λ i = β 0 + f 1 DW T i + f 2 CT H i + f 3 eas t i , nort h i
(7)

with g(.) : link-function: natural logarithm, and GD i ~ Poisson (λ i ) with E(GD i ) = λ i and Var (GD i ) = λ i ; GD i  = 0, 1, 2,….

GD i : larva density at plot i (n · m−2)

DWT i : simulated distance to water table in October 2007 at plot i (m)

CTH i : modeled clay thickness at plot i (%)

east i , north i : Gauß-Krüger east and north coordinates of plot i defined in relation to the 3rd meridian

f 1, f 2 : 1-dimensional smooth functions (penalized thin plate regression splines)

f 3 : 2-dimensional smooth function (penalized thin plate regression spline)

Based on a model that employed only the average larva density (Eq. 7.1), model selection was conducted by the stepwise inclusion of predictors to determine the extent to which they improved the model (Table 2). In the first simple validation, the Akaike information criterion (Burnham and Anderson [2004]) was used. In the process of model selection, it became evident that the geographic location of observations affected larva density to a much higher extent than the causal predictors CTH and DWT (Table 2). By increasing the dimension of the 2-dimensional spatial trend function, the model was improved significantly. However, the empirical variance was significantly larger than the mean. The dispersion parameter decreased with increasing model complexity, but, for all models, the overdispersion was highly significant (Eqs. 7.1–7.7). Therefore, the assumption that larva density distribution could be described by a Poisson distribution was not valid.

Table 2 Various Poisson regression models validated during model selection

In modeling the overdispersion, analogous regression models were fitted assuming that the larva density followed a negative binomial distribution (Eqs. 8a/8b), and a zero-inflated Poisson distribution (Eqs. 9a/9b). Again model improvement was validated by a stepwise integration of predictors, and by increasing the basis dimension of the 2-dimensional spatial trend function (Tables 3 and 4). The specific negative binomial regression model for predicting larva density per m2 using all available predictor variables was defined as follows (The abbreviations correspond to those given for model 7):

g 1 μ i = β 01 + f 11 DW T i + f 21 CT H i + f 31 eas t i , nort h i
(8a)
g 2 1 / ϕ i = η 2 i =x ' i β 2 = β 02 + f 12 DW T i
(8b)
Table 3 Various negative binomial regression models validated during model selection
Table 4 Various zero-inflated Poisson regression models validated during model selection

with g 1(.) and g 2(.): link-functions: natural logarithm

and G D i ~ NegBin μ i , ϕ i with E G D i = μ i and Var G D i = μ i + μ i 2 / ϕ i ; G D i = 0 , 1 , 2 ,

The AIC values of the negative binomial models were much lower than those of Poisson models with approximately the same complexity of the linear predictor. An optimum value was reached at approximately 120 degrees of freedom (edf) for the spatial trend function (Eq. 8.5). Hence the optimum model was much more parameter parsimonious than the best Poisson model (Eq. 7.7). The model was improved slightly when the dispersion parameter was modeled as a function of the DWT as well (Eq. 8.6). In the best Poisson model, the dimension of the spatial trend was found to be ~ 320 edf, but AIC was much higher (Table 2) than the one of the best negative binomial model. Even with this high level of complexity, an optimum degree of smoothing was not reached (Eq. 7.7). Model selection using the Poisson regression approach was discontinued at this stage because more complex models no longer converged. Using all available predictor variables, a specific negative zero-inflated Poisson regression model for predicting larva density per m2 was defined as follows:

g 1 λ i = β 01 + f 11 DW T i + f 21 CT H i + f 31 eas t i , nort h i
(9a)

with g 1(.): link-function: natural logarithm

g 2 ω i = β 02 + f 12 DW T i
(9b)

with g 2(.): link-function: logistic link function

andG D i ~ZIP λ i , ω i
with E G D i = 1 ω i λ i and Var G D i = λ i 1 ω i 1 + λ i ω i G D i = 0 , 1 , 2 , .

(For a description of abbreviations see legend for model 7)

The AIC values in the zero-inflated Poisson models were also considerably lower than in the Poisson regression models, but higher than in the associated negative binomial models (Table 4).

The model effects and their linear approximations were only illustrated for the best- suited model (Eq. 8.5). Model 8.5 was chosen as the best model because further improvement by employing model 8.6 was negligible. The integration of the covariates DWT and CTH in the model both resulted in significant model improvements (Table 3). The effect of DWT showed a significantly non-linear pattern (Figure 2A). The effect of CTH indicated a non-linear pattern also (Figure 2B). However, in the latter case the deviation from a simple linear model was not significant. Additional stratification when using different substrate groups like sand, loam and loess as predictors did not improve the model. Between ~ 0 m and 4 m DWT, the model effect increased (Figure 2A). Hence the predicted larva density increased as the groundwater level dropped. From a DWT of ~ 4 m, an approximately asymptotic pattern resulted, indicating that DWT had no further effect on larva density once it exceeded 4 m. The mean larva density decreased with increasing clay thickness. However, the confidence intervals for the effect of CTH were large, especially for sites with high CTH (Figure 2B).

Figure 2
figure 2

Model effects (Eq. 8.5) of distance to water table ( DWT ) (A) and clay thickness ( CTH ) (B) on the mean of the conditional negative binomial distribution for describing larva density per m 2 . Dashed lines mark the 95% confidence intervals of the model effects.

For simplification the model (Eq. 8.5) was reparameterized by using segmented linear functions to approximate the effects of DWT and CTH on mean larva density (Eq. 8.51, Table 5). In doing so the original effect of DWT was approximated by a segmented function that resulted in an increasing effect of DWT down to a depth of 4 m, and a constant effect at depths below 4 m (Figure 3A). The effect of CTH was approximately constant for values below 30%, and decreased for values higher than 30% (Figure 3B). The values of exponents of the non-constant segments were selected iteratively, and the models compared using the AIC (Eq. 8.51). The assumption of partially constant effects is statistically justified since the confidence intervals would include straight lines (Figure 2). The simplification resulted in an AIC of 4150.15, which means that the model was improved further by the linear approximations. Hence the gam (Eq. 8.5) was simplified to a generalized linear model (glm) (Eq. 8.51):

g 1 μ i = β 01 + f 11 DW T i + f 21 CT H i + f 31 eas t i , nort h i , g 2 1 / ϕ = β 02
(8.51)
Table 5 Statistical characteristics of the generalized linear model (Eq.8.51) for estimating larva density/m 2
Figure 3
figure 3

Linear approximations of model effects (Eq8.51) of distance to water table ( DWT ) (A) and clay thickness ( CTH ) (B) on the mean of the conditional negative binomial distribution for describing larva density per m 2 . Dashed lines mark the 95% confidence intervals of the model effects.

with g 1(.) and g 2(.) : link-functions: natural logarithm

with f 11 DW T i : β 11 DW T i 4 1.5 , for DW T i < 4 m , 0 otherwise and f 21 CT H i : b 2 1 CT H i 30 1.8 , for CT H i > 30 % , 0 otherwise and G D i ~ NegBin μ i , ϕ with E G D i = μ i and Var G D i = μ i + μ i 2 / ϕ ; G D i = 0 , 1 , 2 , .

(For a description of abbreviations see legend for model 7)

3.2 Validation of the distribution assumption

The models with a spatial trend function of approximately 125 edf for each of the different distribution assumptions (Eqs. 7.5, 8.5, 9.6) were validated in more detail by qqplots (Figure 4) of the randomized quantile residuals. Additionally the negative binomial model which included linearly approximated model effects (Eq. 8.51) was validated (Figure 4D).

Figure 4
figure 4

qqplots of randomized quantile residuals for different models for quantifying larva density per m 2 (Eqs. 7.5, A; 8.5, B; 9.6, C;8.51, D).

The qqplots indicated that the negative binomial regression (Eqs. 8.5/8.51) models were best for quantifying larva density, and they confirmed the ranking resulting from the AIC-values (Tables 2, 3, 4 and 5). For both negative binomial regression models, the randomized quantile residuals lay on the bisecting line, indicating that they follow approximately a standard normal distribution, and hence they represent an almost optimal fit to the data. The distribution of quantile residuals of model 8.5 was caused only marginally by the linear approximations adopted (Eq. 8.51). The qqplots for the Poisson and zero-inflated regression models displayed major deviations from the bisecting line. Hence, neither the assumption that larva density follows a Poisson distribution nor the assumption that larva density follows a zero-inflated distribution could be validated in this case.

The wormplots for models 8.5 and 8.51 showed that the distributions of randomized quantile residuals deviated marginally from a standard normal distribution (Figure 5). An optimal model fit would be characterized by a horizontal line through the origin. However in the case of model 8.5 none of the deviations were significant to the 95% level since all points lay within the dashed boundary lines of the confidence interval of the standard normal quantiles (Figure 5A). For model 8.51 the overall pattern of deviations was even better, since the line flattens out, but some deviations were significant to the 95% level (Figure 5B).

Figure 5
figure 5

Overall wormplots of negative binomial regression models (Eq. 8.5, A; Eq.8.51, B) for quantifying larva density. Curved dashed lines indicate the boundaries of the 95% confidence interval of the standard normal quantiles. The curved solid gray line describes a polynomial trend curve that highlights a potential pattern in the randomized quantile residuals.

Conditional wormplots for model 8.51 showed that the model fit was good, even with respect to different intervals of the predictor variable DWT (Figure 6). Marginal deviations from standard normal quantiles occurred mainly for DWT values less than 2 m (lower left panel) and between 2 m and 3 m (lower right panel). Conditional wormplots for groups of CTH were omitted since the number of plots with a clay occurrence was low.

Figure 6
figure 6

Conditional wormplots of the negative binomial regression model (Eq.8.51) for different groups of the predictor variable DWT . Deviations from standard normal quantiles are displayed for DWT values less than 2 m, 2–3 m, 3–4 m and more than 4 m in depth (from lower left panel to upper right panel). Curved dashed lines indicate the boundaries of the 95% confidence intervals of the standard normal quantiles. The curved solid gray lines are polynomial trend curves that highlight potential pattern in deviation from standard normal quantiles.

Even when the negative binomial regressions (Eqs. 8.5/8.51) were found to fit the data almost perfectly, in the extensive validation procedure some artefacts in the spatial trend function became evident in graphical representations of the model predictions (Figure 7A). These artefacts with very high values occurred in areas comprising only very few or no sample plots, and therefore can be characterized as edge effects. To resolve this problem a different type of thin plate spline was employed that uses 1st order derivatives in the penalty (Duchon [1977]) for the 2-dimensional spatial trend function (Wood [2006]). Additionally the dimension of the spatial trend function was decreased slightly (edf = 93.38). The resulting trend surface was smoother and the artefacts of extreme high-trend values disappeared (Figure 7).

Figure 7
figure 7

Contour plots displaying the spatial pattern of predicted larva density assuming a CTH of 0% and a DWT of 4 m employing negative binomial regression models (Eq.8.51) with conventional thin plate splines (A) and thin plate splines with 1 st order derivative penalty and slightly reduced dimension of the spatial trend (B). The transect highlighted by a dashed line connects the three geographic locations used in the model computations. The light gray triangle represents a location with a low, the dark gray circle with a medium, and the black square with a high spatial effect. The outermost isolines of the three distinct areas connect locations with predictions of 2 larvae · m−2.

However, these constraints led to a slight increase in the AIC (4174.69 versus 4150.15), but the qqplot of the randomized quantile residuals still characterized this more parsimonious model variant (Eq. 8.51 Duchon) as an almost perfect fit to the data (Figure 8).

Figure 8
figure 8

The qqplot of randomized quantile residuals for a negative binomial model for quantifying larva density (Eq.8.51) using 1 st order derivatives in the penalty (Duchon[1977]) and reducing the dimension of the spatial trend function ( edf= 93.38).

3.3 Validation of spatial independency

The randomized quantile residuals of the finally selected model (Eq. 8.51 Duchon) using 1st order derivatives in the penalty and three other negative binomial regression models were compared and validated for spatial independency. The validation was conducted for the model including no spatial trend function (Eq. 8.3), the model with a spatial trend function of low dimension (Eq. 8.4) and the model with a trend function of optimal dimension concerning AIC (Eq. 8.51). Hence the finally selected model was compared to two models of lower complexity (Eqs. 8.3, 8.4) and one model of higher complexity (Eq. 8.51) of the spatial trend function. For comparison models 8.3 and 8.4 were refitted by employing approximated effects for the causal predictors DTW and CTH also (see Eq. 8.51).

The residuals from the model without any spatial trend function (Eq. 8.3) showed a significant deviation from spatial independency up to a distance of approximately 7 500 m (Figure 9A). The implementation of a spatial trend function of low dimension into the model (Eq. 8.4) removed the spatial dependencies for larger distances but significant deviations remained up to a distance of approximately 1 500 m (Figure 9B). The residuals for the models with more complex spatial trend functions (Eq. 8.51 and Eq. 8.51 Duchon) showed no spatial dependencies anymore (Figure 9C and 9D). Hence the results from the empirical semi-variograms confirmed the results from the model selection process. Moreover the slightly more parsimonious and finally selected model (8.51 Duchon) is sufficient to describe the spatial pattern in larva density adequately.

Figure 9
figure 9

Empirical semi-variograms for the randomized quantile residuals of four negative binomial regression models (Eqs. 8.3, A; 8.4, B;8.51, C;8.51Duchon, D) with spatial trend functions of different complexity (solid black dots). Additionally pointwise 95% confidence intervals for the null hypothesis of spatial independency are plotted that are based on 1 000 random permutations of the coordinates (dashed lines with blank dots).

3.4 Model computations

The model sensitivity resulting from the different model effects can be illustrated by predictions on the scale of the response variable (larva density) by varying only one of the causal predictors DWT or CTH while keeping the other constant. For this prediction the negative binomial model (Eq. 8.51) employing 1st order derivatives in the penalty was used. The geographic location was kept constant also because of the high spatial variation in larva density (Figure 7B). To provide an example, the larva density was estimated at three different locations differing significantly in spatial effects (Figure 10). The large differences between the curves in the figures resulted from these different geographic locations: the weaker the causal effects were, the higher the variation in predicted larva density for the different locations was. Hence, if conditions like a low CTH or a large DWT led to a high potential larva density, then the actual larva density depended strongly on geographic location (Figures 10A, 10D). With increasing CTH and declining DWT, the potential of a site to serve as habitat for larvae decreased, and the effect of geographic location also decreased.

Figure 10
figure 10

Predicted larva density (Eq.8.51) using 1 st order derivatives in the penalty at three geographic locations in the Hessian Ried by varying one of the predictors DWT or CTH while keeping the other constant. The value of the constant used is given in the title of each figure. The black, solid line represents predictions for the location coded by a black square, the dark gray, dashed line for the location coded by a dark gray circle, and the light gray, point dashed line for the location coded by a light gray triangle (Figure 7B).

The prediction accuracy of the model is illustrated exemplarily by calculating standard errors for the predictions along a northing gradient through the Hessian Ried and keeping the other predictors constant (Figure 11).

Figure 11
figure 11

Predicted larva density (Eq.8.51Duchon) using 1 st order derivatives in the penalty along a northing gradient through the Hessian Ried for a fixed easting coordinate = 3469000 and constant values of DWT= 4.0 m and CTH= 0. The black, solid line represents predictions for the expectation value and the dashed lines for confidence interval.

4Discussion

It is known that sites with a high groundwater table or high percentage of bed rock prevent the hibernation of forest cockchafer larvae at greater soil depths in cold winter climates (Schwerdtfeger [1981]). This is in part reflected in the model effects. A pure clay layer can be assumed to have similar negative effects on the suitability of a site as a habitat for larvae and an increasing proportion of clay layer leads to reduced numbers of predicted larvae per m2 (Figure 10). In this context the effect of CTH was found to be a weak indicator so far, since the confidence intervals were rather wide. However, the integration of CTH into the negative binomial regression model (Eq. 8.5) led to a reduction of the AIC, and hence improved the model (Table 3).

The predicted larva density within 1 m soil depth is affected by values of DWT up to 4 m (Figure 2A), which does not seem feasible since maximum depth for hibernation is thought to be approximately 1.1 m (Schwerdtfeger [1981]). However, DWT values used in this model were based on ground water data simulated for the month of October 2007. In spring, a lower DWT is usually observed, and single high water occurrences might result temporarily in even smaller DWT. Furthermore, depending on the specific soil substrates, a capillary ascension may reduce the capacity of the first meter of soil to serve as habitat for larvae as well. These interpretations are based on the assumption that single temporary high-water events affect the capacity of the larval habitat significantly. The constraint of a constant effect of DWT of more than 4 m (Eq. 8.51/Figure 3A) was imposed because it is biologically feasible that the effect of DWT is constant below a certain threshold. Higher values of DWT are at least partially the result of intensive groundwater withdrawals in the Hessian Ried (NW-FVA, Nordwestdeutsche Forstliche Versuchsanstalt Hrsg. [2013]): p. 30 ff.). Hence, the lowering of groundwater in areas where groundwater was formerly available to trees may also affect larva density indirectly. The degradation of forest is a direct result of a lowering of the groundwater table, and the partial dieback of single trees and forests may have improved habitat conditions for the forest cockchafer as well.

A comparison of models showed that the negative binomial distribution was superior to the Poisson and zero-inflated Poisson distributions, which is in accordance with other investigations of overdispersed animal count data (Gray [2005]; Sileshi [2008]; Vaudor et al. [2011]). In many investigations the negative binomial, the Poisson and the zero-inflated Poisson distributions have been compared. In some cases a zero-inflated modification of the negative binomial distribution results in minor, and somewhat doubtful improvements in the models (Gray [2005]; Vaudor et al. [2011]). Especially in cases with the occurrence and abundance of a species resulting from distinct processes, the zero-inflated modification of the negative binomial distribution can result in considerable model improvements (Wenger and Freeman [2008]).

Overdispersion is often assumed to result from a spatial or temporal heterogeneity of the habitat. However, even if conditional distributions are fitted by employing regression approaches (Sileshi [2008]) or stratification (Vaudor et al. [2011]), the negative binomial approach has been found to be superior in many investigations. In this investigation the modeling approach for determining the effect of spatial heterogeneity on overdispersion was much more flexible due to the complex spatial trend function (Table 2). However, even though the dimension of the spatial trend was increased considerably, a significant overdispersion was still evident (Table 2). Hence, at least for our investigation, it can be concluded that even a quite complex linear predictor is not sufficient to cover all sources of overdispersion.

Extended generalized regression models facilitate the estimation of the conditional mean, variance or mixture parameter as functions of covariates. The simplest structure of a linear predictor is to assume the effects of the covariates to be linear (Sileshi [2008]). Yet this assumption must be validated to guarantee unbiased predictions across the whole range of covariates (Hastie and Tibshirani [1990]). Therefore the extension of generalized additive models to overdispersed and zero-inflated count data, such as were used in this investigation, indicates a major advance in the methodology (Barry and Welsh [2002]; Rigby and Stasinopoulos [2005]; Wood [2006]).

Due to the heterogeneity of covariates that are spatially correlated, but unknown, or only insufficiently available, in many cases the response data are autocorrelated. Haining et al. ([2009]) presented a simple conditional autoregressive model to deal with autocorrelated count data which results in the estimation of spatially correlated random effects. Fahrmeir and Echavarrı ([2006]) introduce an extensive methodology using structured additive regression models STAR for overdispersed and zero-inflated count data. These models make it possible to model non-linear covariate effects, individual or cluster-specific uncorrelated random effects, spatially correlated random effects or 2-dimensional spatial trend surfaces simultaneously. The methodology employed in our investigation (Wood [2006]) offers similar technical possibilities for the Poisson and negative binomial distributions and, in combination with Rigby and Stasinopoulos’s ([2005]) methods, for the zero-inflated Poisson distribution. The inventory plots are located exactly via coordinates, and hence a 2-dimensional spatial trend function is fitted instead of spatially correlated random effects for distinct areas. The observations at the 4 subplots were aggregated due to their proximity, which makes the estimation of uncorrelated random effects at the plot level redundant. Overall, the approach adopted combines an appropriate distribution assumption for count data with non-linear effects of causal covariates and an advanced method for covering spatial autocorrelation.

Finally graphical representations of randomized quantile residuals are powerful validation tools that provide more detailed information about model characteristics than just comparing global statistics. However, so far, in most investigations of count data, validation has been limited to global statistics or simple comparisons of counts (Gray [2005]; Sileshi [2008]; Wenger and Freeman [2008]; Vaudor et al. [2011]).

5Conclusions

The developed negative binomial regression model can be used to predict the current spatial pattern of larva density in the Hessian Ried. The almost optimal fit by the model allows for the prediction of conditional expectation values but also of conditional quantiles to account for uncertainty in the process of silvicultural decision making. Based on the predictions local forest managers will be able to optimize the spatial pattern of regeneration and forest protection planning. Areas with different risks can be identified by combining the predictions with expert knowledge about critical larva density. Hence areas can be separated where forest protection measures are essential, reasonable or needless to ensure successful regeneration measures. A classification of forest stands can be conducted that will be mainly affected by their regional location within the Hessian Ried. Deviations from the overall spatial pattern will be affected by the spatial pattern of pure clay layer and distance to water table.

In the future inventories of larvae in the Hessian Ried will be conducted continually at an interval of 4 years in the year prior to the year of swarming. In the planning of control measures, it is of particular interest to know if the spatial distribution of the larvae varies over time or is more or less stationary. Hence, a major objective of the future inventories will be to gain insight into the spatio-temporal pattern of larva density. For this purpose the methodology developed can be used for time series data by estimating inventory-specific spatial trends or, in the case of a number of inventories, by integrating a space-time effect (Augustin et al. [2009]).

The overall sampling grid should be optimized in future inventories to enable more plots with a medium to high CTH and low DWT to be assessed. In this context the database could be improved considerably by recording the CTH when excavating instead of modeling it.

Some knowledge exists about the capacity of different stand structures to serve as habitat for forest cockchafer larvae. For example dense young stands provide a less suitable habitat (Schwerdtfeger [1981]). Hence additional information about stand structure and tree species composition should be recorded and its effects on larva density tested. Currently the model predictions might be biased for certain stand structures if relevant effects on larva density do exist. However it is likely that the strong spatial effect will be the most important effect in future model developments also. Finally the 2013 inventory will be used to test the effect of forest protection measures that have been implemented in a sub-area of the Hessian Ried.

References

  • Augustin NH, Musio M, Wilpert K, Kublin E, Wood SN, Schumacher M: Modeling spatiotemporal forest health monitoring data. J Am Stat Assoc 2009, 104(487):899–911. doi:10.1198/jasa.2009.ap07058 doi:10.1198/jasa.2009.ap07058 10.1198/jasa.2009.ap07058

    Article  CAS  Google Scholar 

  • Barry SC, Welsh AH: Generalized additive modelling and zero inflated count data. Ecol Model 2002, 157(2–3):179–188. 10.1016/S0304-3800(02)00194-1

    Article  Google Scholar 

  • Burnham KP, Anderson DR: Multimodel inference: understanding AIC and BIC in model selection. Sociol Method Res 2004, 33(2):261–304. doi:10.1177/0049124104268644 doi:10.1177/0049124104268644 10.1177/0049124104268644

    Article  Google Scholar 

  • Buuren S, Fredriks M: Worm plot: simple diagnostic device for modelling growth reference curves. Stat Med 2001, 20: 1259–1277. 10.1002/sim.746

    Article  PubMed  Google Scholar 

  • Cameron AC, Trivedi PK: Regression-based tests for overdispersion in the Poisson model. J Econometrics 1990, 46(3):347–364. 10.1016/0304-4076(90)90014-K

    Article  Google Scholar 

  • Duchon J: Splines minimizing rotation-invariant semi-norms Sobolev spaces. Construction Theory of Functions of Several Variables. Springer, Berlin, In; 1977.

    Book  Google Scholar 

  • Dunn PK, Smyth GK: Randomized quantile residuals. J Comput Graph Stat 1996, 5: 236–244.

    Google Scholar 

  • Fahrmeir L, Echavarrı LO: Structured additive regression for overdispersed and zero-inflated count data. Appl Stochastic Models Bus Ind 2006, 22: 351–369. 10.1002/asmb.631

    Article  Google Scholar 

  • Fahrmeir L, Kneib T, Lang S: Regression. Springer, Berlin; 2007.

    Google Scholar 

  • Gray BR: Selecting a distributional assumption for modelling relative densities of benthic macroinvertebrates. Ecol Model 2005, 185: 1–12. 10.1016/j.ecolmodel.2004.11.006

    Article  Google Scholar 

  • Haining R, Law J, Griffith D: Modelling small area counts in the presence of overdispersion and spatial autocorrelation. Comput Stat Data Anal 2009, 53: 2923–2937. 10.1016/j.csda.2008.08.014

    Article  Google Scholar 

  • Hastie HJ, Tibshirani RJ: Generalized Additive Models. Monographs on Statistics and applied Probability 43. Chapman & Hall, London, New York, Tokyo, Melbourne, Madras; 1990.

    Google Scholar 

  • Kleiber C, Zeileis A: Applied Econometrics with R. Springer-Verlag, New York; 2008.

    Book  Google Scholar 

  • McCullagh P, Nelder JA: Generalized Linear Models, 2nd edn. Monographs on Statistics and applied Probability 37. London, New York, Tokyo, Melbourne, Madras, Chapman and Hall; 1989.

    Google Scholar 

  • NW-FVA, Nordwestdeutsche Forstliche Versuchsanstalt (Hrsg.) (2013) Waldentwicklungsszenarien für das Hessische Ried, Entscheidungsunterstützung vor dem Hintergrund sich beschleunigt ändernder Wasserhaushalts- und Klimabedingungen und den Anforderungen aus dem europäischen Schutzgebietssystem Natura 2000. Beiträge aus der Nordwestdeutschen Forstlichen Versuchsanstalt, Band 10. Universitätsverlag Göttingen, p 397 , [http://webdoc.sub.gwdg.de/univerlag/2013/NWFVA10_HessischesRied.pdf]

  • Ott A, Delb H, Mattes J, Schröter H: Erfolgreiche Regulierung eines Nebenflugstammes des Waldmaikäfers. AFZ-DerWald 2006, 61(6):312–314.

    Google Scholar 

  • R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria; 2010.

  • Rigby RA, Stasinopoulos DM: Generalized additive models for location, scale and shape (with discussion). Appl Statist 2005, 54: 507–554.

    Google Scholar 

  • Rigby RA, Stasinopoulos DM (2009) A flexible regression approach using GAMLSS in R. Short course booklet, [http://book.gamlss.org/]

    Google Scholar 

  • Schwerdtfeger F: Die Waldkrankheiten. 4. neubearbeitete Auflage, Verlag Paul Parey, Hamburg und Berlin; 1981.

    Google Scholar 

  • Sileshi G: The excess-zero problem in soil animal count data and choice of appropriate models for statistical inference. Pedobiologia 2008, 52: 1–17. 10.1016/j.pedobi.2007.11.003

    Article  Google Scholar 

  • Vaudor L, Lamouroux N, Olivier J-M: Comparing distribution models for small samples of overdispersed counts of freshwater fish. Acta Oecol 2011, 37: 170–178. 10.1016/j.actao.2011.01.010

    Article  Google Scholar 

  • Venables WN, Ripley BD: Modern Applied Statistics with S. Springer, New York; 2002.

    Book  Google Scholar 

  • Wenger SJ, Freeman MC: Estimating species occurrence, abundance, and detection probability using zero-inflated distributions. Ecology 2008, 89: 2953–2959. 10.1890/07-1127.1

    Article  PubMed  Google Scholar 

  • Wood SN: Generalized Additive Models: An Introduction with R. Chapman and Hall/CRC, Boca Raton; 2006.

    Google Scholar 

  • Wood SN, Augustin NH: GAMs with integrated model selection using penalized regression splines and applications to environmental modelling. Ecol Model 2002, 157: 157–177. 10.1016/S0304-3800(02)00193-X

    Article  Google Scholar 

  • Yee TW: The VGAM package for categorical data analysis. J Statist Software 2010, 32(10):1–34. http://www.jstatsoft.org/v32/i10/ http://www.jstatsoft.org/v32/i10/

    Article  Google Scholar 

  • Yee TW, Wild CJ: Vector generalized additive models. J R Statist Soc B 1996, 58(3):481–493.

    Google Scholar 

  • Zucchini W, Schmidt M, Gadow K: A model for the diameter-height distribution in an uneven-aged beech forest and a method to assess the fit of such models. Silva Fennica 2001, 35: 169–183. 10.14214/sf.594

    Article  Google Scholar 

Download references

Acknowledgements

We are thankful to Johannes Sutmöller for providing the groundwater data and to Falko Engel for calculating the percentage clay thickness from soil substrate maps. We would like to thank Peter Gawehn and his team for their comprehensive and careful field work and data collection and Helen Desmond for her extensive language revision. We also thank two anonymous reviewers for valuable comments that improved the quality of the manuscript. We also thank the Hessian State Forest Enterprise, Hessen-Forst, the principal project partner.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Matthias Schmidt.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

RH provided the data. MS conducted the analyses. Both authors contributed to the writing of the manuscript. Both authors read and approved the final manuscript.

Authors’ original submitted files for images

Below are the links to the authors’ original submitted files for images.

Authors’ original file for figure 1

Authors’ original file for figure 2

Authors’ original file for figure 3

Authors’ original file for figure 4

Authors’ original file for figure 5

Authors’ original file for figure 6

Authors’ original file for figure 7

Authors’ original file for figure 8

Authors’ original file for figure 9

Authors’ original file for figure 10

Authors’ original file for figure 11

Authors’ original file for figure 12

Authors’ original file for figure 13

Authors’ original file for figure 14

Authors’ original file for figure 15

Authors’ original file for figure 16

Authors’ original file for figure 17

Authors’ original file for figure 18

Authors’ original file for figure 19

Authors’ original file for figure 20

Authors’ original file for figure 21

Authors’ original file for figure 22

Authors’ original file for figure 23

Authors’ original file for figure 24

Authors’ original file for figure 25

Authors’ original file for figure 26

Authors’ original file for figure 27

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (https://creativecommons.org/licenses/by/4.0), which permits use, duplication, adaptation, distribution, and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Schmidt, M., Hurling, R. A spatially-explicit count data regression for modeling the density of forest cockchafer (Melolontha hippocastani) larvae in the Hessian Ried (Germany). For. Ecosyst. 1, 19 (2014). https://doi.org/10.1186/s40663-014-0019-y

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40663-014-0019-y

Keywords