Strategies for Modeling Regeneration Density in Relation to Distance from Adult Trees

Research Highlights: We proposed new methodologies for the spatial analysis of regeneration processes and compared with existing approaches. Background and Objectives: Identifying the spatial relationship between adult trees and new cohorts is fundamental to understanding the dynamics of regeneration and therefore helps us to optimize the stand density and natural regeneration when undertaking regeneration fellings. Most of the statistical approaches analyzing the spatial dependence between adult trees and new individuals (seedlings or saplings) require a complete census and mapping of all individuals. However, approaches considering individuals grouped into sampling points or subplots (i.e., density data) are limited. In this study, we reviewed and compared approaches (intertype point pattern analyses and a generalized additive model) to describe the spatial relationship between adult trees and density regeneration in a Pinus sylvestris L. monospecific stand in Spain. We also proposed a new approach (intertype mark variance function) to disentangle the effect of the tree-size on sapling density and the effect of the spatial pattern. Materials and Methods: To this end, we used a half-hectare plot in which all the individuals of P. sylvestris have been mapped and measured. Results: Our results indicated that sapling distribution was related to distance from the adult trees, thus displaying distance-dependence patterns, but it was not related to the size of the adult trees. The intertype mark correlation function was an useful tool to distinguish the effect of the marks (sapling density and tree size) from the effect of the spatial pattern of the classes (trees cohorts in our case). Conclusions: The largest number of saplings was found with increased distance between adult trees (>11 m), and the generalized additive model may be useful to explain spatial relationships between adult trees and regenerating cohorts when other measured biotic variables (e.g., soil stoniness, etc.) and repeated measurements are available.


Introduction
The establishment of new individuals following regeneration fellings is a vital stage of stand development as the persistence of the forest cover is at stake [1]. The occurrence and survival of the new individuals is hindered by several factors, such as low seed production, post-dispersal predation, summer drought, limited nutrient availability, flooding, frost, fire, competition with grass and shrubs, lack of or excess light, or damage from browsing, pests and pathogens [1][2][3][4][5]. Furthermore, the impact of these factors on the regeneration process can be modified under climate change conditions [6,7]. Therefore, a greater understanding of the relationship between the regeneration process and environmental conditions is crucial for current and future forest management [7]. In this regard, silvicultural operations can modulate light and water availability to the new individuals, i.e., seedlings and saplings, through reducing stem density over the regeneration period [8].

Study Site
The study was conducted at Navafría forest (41 • 00 N, 3 • 48 W), in the Central Mountain Range (Spain). The forest grows on acid soils and is mainly north-facing. The altitude of the forest ranges from approximately 1200 to 2200 m above sea level. The annual rainfall is about 1000 mm and the mean temperature is around 9.8 • C. The main species of the forest is P. sylvestris but Quercus pyrenaica Willd. is also present at lower altitudes. The rotation period for P. sylvestris is 100-years. At Navafría, the regeneration of P. sylvestris is achieved through the uniform shelterwood method with three fellings (a preparatory cutting, an establishment cutting and a final cutting) over a 20-year regeneration period. At the beginning of the regeneration period, the density is around 200-300 trees per hectare [25,30] and the occurrence of seedlings is promoted when the density is around 100 trees per hectare after the establishment cutting [29]. Soil preparation operations, subsoil and blade scarification are commonly carried out to favor the establishment of the regeneration after the establishment cutting [25,31].
In 2001, a half-hectare square plot was installed in a monospecific, even-aged stand of P. sylvestris in Navafría forest in order to study the structure and dynamics of this species. The plot was installed at the beginning of the regeneration period when the preparatory felling, which usually takes place about 5-10 years before the establishment felling, had already been carried out. The silvicultural operations in the plot are carried out by the forest service following the forest management plan of Navafría forest; thereby the plot is representative of the uniform shelterwood application to P. sylvestris stands in Central Spain. In this plot, we labelled and measured the diameter at breast height and the height of all the stems with a height ≥ 1.30 m. We distinguished two cohorts, mother trees (individuals with dbh ≥ 10 cm) and saplings (individuals with a height ≥ 1.30 m and dbh < 10 cm). The mother trees were individually mapped whereas the saplings were grouped into a 2 × 2 m quadrat grid. These measurements were repeated in 2006, 2010 and 2014 (Table 1).
The second felling under the uniform shelterwood system, the establishment felling, was undertaken between the first (2001) and the second inventory (2006), and soil preparation was carried out in 2010, with seedling establishment taking place between 2010 and 2014. For the purposes of this study, we only considered the 2014 inventory, when both saplings and mother trees were present ( Figure 1).

Spatial Pattern of Pinus Sylvestris Saplings and Adult Trees
We used two local quadrat variance methods for the two dimensional data to tackle the spatial distribution of the P. sylvestris sapling density: the four-term local quadrat variance (4TLQV) and the nine-term local quadrat variance (9TLQV) [32]. These methods calculate the variance at different spatial scales estimating the variability among neighbor blocks. The block size is progressively increased by combining the number of quadrats. Thus, the block is the search window [33]. The squared differences in the number of saplings between each group of four and nine neighbor blocks were averaged for the whole study area for the 4TLQV and 9TLQV, respectively [18,32].
The spatial distribution of adult trees was studied using the ˆ( ) L d transformation of Ripley´s function K(d). ˆ( ) L d linearizes the K(d) function and gives 0 as the value under the random distribution of mother trees [15].
where λ is the density of mother trees per unit area, dij is the distance between the i and the j adult tree. dij(d) was replaced by ωij(d) in order to account for the edge effect. ωij(d) is 0 if adult trees i and j are farther than distance d from each other, and otherwise, ωij(d) is the inverse of the fraction of a circumference centered on i (with the radius being the distance from i to j) which falls within the plot [15,34]. We simulated 999 patterns of the complete spatial randomness null model in order to build 95 % quantile bounds to test statistically the random distribution of adult trees.

Spatial Pattern of Pinus Sylvestris Saplings and Adult Trees
We used two local quadrat variance methods for the two dimensional data to tackle the spatial distribution of the P. sylvestris sapling density: the four-term local quadrat variance (4TLQV) and the nine-term local quadrat variance (9TLQV) [32]. These methods calculate the variance at different spatial scales estimating the variability among neighbor blocks. The block size is progressively increased by combining the number of quadrats. Thus, the block is the search window [33]. The squared differences in the number of saplings between each group of four and nine neighbor blocks were averaged for the whole study area for the 4TLQV and 9TLQV, respectively [18,32].
The spatial distribution of adult trees was studied using theL(d) transformation of Ripley's function K(d).L(d) linearizes the K(d) function and gives 0 as the value under the random distribution of mother trees [15].L where λ is the density of mother trees per unit area, d ij is the distance between the i and the j adult tree. d ij (d) was replaced by ω ij (d) in order to account for the edge effect. ω ij (d) is 0 if adult trees i and j are farther than distance d from each other, and otherwise, ω ij (d) is the inverse of the fraction of a circumference centered on i (with the radius being the distance from i to j) which falls within the plot [15,34]. We simulated 999 patterns of the complete spatial randomness null model in order to build 95 % quantile bounds to test statistically the random distribution of adult trees.

Intertype K rx (d) Function
We used the intertype K rx (d) function proposed by Montes and Cañellas [18] to study the spatial dependence of sapling density on the spatial distribution of adult trees. This function modifies the interclass K 12 (d) function [32] allowing positive, negative or independent spatial associations between a point pattern of the adult trees and the sapling density measured at sampling points (or quadrat centers) and the scale at which it occurs to be identified [13]. Montes and Cañellas [18] expressed the K rx (d) function as: where N is the number of adult trees, n is the number of quadrats where the saplings were measured, x j is the number of saplings in the quadrat j and x is the mean sapling density in the study area. ω ij is calculated as in the K (d) function. Finally, we simulated 999 patterns of the antecedent condition null model and 999 patterns of the toroidal shift null model to build 95 % quantile bounds to test statistically the spatial dependence between adult trees and saplings [35].

Intertype Mark CorrelationK rs mm (d) Function
The intertype mark correlationK rs mm (d) function [36] is a conditional mean function that incorporates the mark covariance as a test function, and can be defined as the expected covariance of the normalized mark values between pairs of points of different cohorts within a given distance: where n r and n s are the number of points in class r and class s (adult trees and saplings); s r and s s are the standard deviations of the marks in class r and class s respectively; m r and m s are the mean of the marks in class r and class s respectively and ω is calculated as in the K (d) function. As the function is normalized by the sum of the unweighted ω ij (d), its value for a distance d depends only on the covariance between the marks of points belonging to different classes located at a distance < d.
In our case, the first class corresponds to the adult tree pattern and the dbh of each tree is used as a mark. The second class corresponds to the center of the quadrats where sapling density has been measured, the sapling density being the variable used as the mark.
The toroidal shift null model was used to test independence between the diameter of the adult trees and the sapling density for each distance, keeping the spatial distribution of the marks of both the adult trees and the quadrats fixed but shifting one pattern with respect to the other by a random vector for each simulation. If the diameter of the adult trees and the sapling density are independently distributed, the empirical function would be within the quantile bounds, whereas values for theK rs mm (d) empirical function below the lower bound indicate negative correlation and values for theK rs mm (d) empirical function over the upper bound indicate a positive correlation.

Disentangling the Spatial Pattern from the Mark Correlation: Intertype Mark VarianceL rs v (d) Function
In this section, we propose an approach to disentangling the spatial relationship between the marks of the classes (dbh for adult trees and density for saplings) and the spatial pattern of the two cohorts. The approach proposed here expands that of Montes et al. [37] for univariate data to bivariate data: All the variables have been defined above. The value of the function depends both on the correlation of the marks and the number of individuals within a distance d around those individuals of different classes. The value of the function for the real plot distribution is compared with two different null models, one fixing the position of the points to detect correlation of the marks and the other fixing the spatial arrangement of the marks through a moving window algorithm to isolate the spatial attraction/repulsion. Values for theL rs v (d) function outside the quantile bounds from the toroidal shift null model, preserving the mark associated with each adult tree when shifting their pattern, indicate spatial dependence of the marks of the classes, analogous to the case ofK rs mm (d), whereas values for theL rs v (d) function outside the quantile bounds from the toroidal shift null model combined with the moving window algorithm to reassign the mark associated with the adult trees, maintaining the relative distribution of the marks (dbh) of the adult trees with respect the quadrat densities, indicate spatial dependence of the classes, analogous to the case of K rx (d).
The local quadrat variance, the L(d) function and the intertype analyses were done using a Microsoft ® VisualBasic ® application developed by the authors (available by contacting the authors).

Functional Predictor within a Generalized Additive Model Framework
We modeled the expected number of saplings E(Ns j ) = µ j in the quadrat j employing a functional predictor [38] in a generalized additive model (GAM) framework [39,40]. In the functional prediction, the effect of every adult tree on the number of saplings per quadrat is weighted according to the distance between adult trees and saplings. This approach has been satisfactorily used to model forest ingrowth [20]. We proposed the following model: with Ns ij following a negative binomial distribution as the variance exceeds the mean [41]. α refers to the intercept of the model, D ji is a matrix which contains the distances (in m) from the adult tree (I = 1, . . . , N) to the j sapling quadrat, and dbh i is the matrix of the dbh of the adult tree (I = 1, . . . , N) which dimensions are j files and i columns. The dbh was set to 0 when the i adult tree was more than 30 from the j quadrat.
AreaIn30 is the area (in m 2 ) of the 30 m radius circle within the plot and was entered to correct the edge effect. Thus, this implies accepting the assumption that the surrounding trees outside the plot show a similar distribution to those within the plot. More details of the edge correction can be found in Moreno-Fernández et al. [20]. We select 30 m as it is approximately the double of the mean height of adult trees. Additionally, previous simulations revealed that the estimated smooth coefficient function stabilizes around zero at 30 m. The function f 2 (X j , Y j ) is a spatial smooth term of the j quadrat to account for the spatial trend, spatial correlation of the number of saplings and any spatial trend caused by other unmeasured environmental variables. The selection of the variables was carried out using a stepwise selection process, using Akaike's Information Criterion (AIC) for the evaluation. Finally, we checked whether the assumption of spatial independence was met by plotting a semivariogram per inventory with envelopes from 99 permutations under the assumption of no spatial correlation [42]. These statistical analyses were conducted using "mgcv" [38] and "geoR" [43] in R version 3.6.0 [44].

Spatial Pattern of P. Sylvestris Saplings and Adult Trees over the Regeneration Fellings of the Uniform Shelterwood System
The occurrence of the P. sylvestris saplings took place between 2010 and 2014 (Table 1). In 2014, the saplings were distributed in clumps over the whole plot (except below the adult trees) (Figure 1). The spatial structure of P. sylvestris saplings in clumps, as revealed by the map (Figure 1), is statistically confirmed by the variance-block length graph ( Figure 2). This graph shows that both the 4TLQV and 9TLQV peaked at 5 m, suggesting spatial aggregation of saplings at scales of 5 m. The selection of the variables was carried out using a stepwise selection process, using Akaike's Information Criterion (AIC) for the evaluation. Finally, we checked whether the assumption of spatial independence was met by plotting a semivariogram per inventory with envelopes from 99 permutations under the assumption of no spatial correlation [42]. These statistical analyses were conducted using "mgcv" [38] and "geoR" [43] in R version 3.6.0 [44].

Spatial Pattern of P. Sylvestris Saplings and Adult Trees over the Regeneration Fellings of the Uniform Shelterwood System
The occurrence of the P. sylvestris saplings took place between 2010 and 2014 (Table 1). In 2014, the saplings were distributed in clumps over the whole plot (except below the adult trees) (Figure 1). The spatial structure of P. sylvestris saplings in clumps, as revealed by the map (Figure 1), is statistically confirmed by the variance-block length graph ( Figure 2). This graph shows that both the 4TLQV and 9TLQV peaked at 5 m, suggesting spatial aggregation of saplings at scales of 5 m.
There is evidence of spatial regularity of adult trees at distances between 3 to 9 m as the L(d) function lies below the lower quantile bound (Figure 3). Beyond this range, the function is within the quantile bounds, suggesting a random distribution of adult trees farther than 9 m.   There is evidence of spatial regularity of adult trees at distances between 3 to 9 m as the L(d) function lies below the lower quantile bound (Figure 3). Beyond this range, the function is within the quantile bounds, suggesting a random distribution of adult trees farther than 9 m. The selection of the variables was carried out using a stepwise selection process, using Akaike's Information Criterion (AIC) for the evaluation. Finally, we checked whether the assumption of spatial independence was met by plotting a semivariogram per inventory with envelopes from 99 permutations under the assumption of no spatial correlation [42]. These statistical analyses were conducted using "mgcv" [38] and "geoR" [43] in R version 3.6.0 [44].

Spatial Pattern of P. Sylvestris Saplings and Adult Trees over the Regeneration Fellings of the Uniform Shelterwood System
The occurrence of the P. sylvestris saplings took place between 2010 and 2014 (Table 1). In 2014, the saplings were distributed in clumps over the whole plot (except below the adult trees) (Figure 1). The spatial structure of P. sylvestris saplings in clumps, as revealed by the map (Figure 1), is statistically confirmed by the variance-block length graph (Figure 2). This graph shows that both the 4TLQV and 9TLQV peaked at 5 m, suggesting spatial aggregation of saplings at scales of 5 m.
There is evidence of spatial regularity of adult trees at distances between 3 to 9 m as the L(d) function lies below the lower quantile bound (Figure 3). Beyond this range, the function is within the quantile bounds, suggesting a random distribution of adult trees farther than 9 m.

Intertype Functions
The intertype K rx (d) function for the spatial structure of saplings and adult trees revealed a negative influence of the adult trees on the sapling density when the distance between the individuals of both cohorts is shorter than 11 m as the K rx (d) function is below zero (Figure 4). Beyond this distance, the function remains close to zero indicating an independent spatial pattern between the two cohorts. However, the K rx (d) function is only outside the 95% quantile bounds at distances between 1.5 and 8.5 m in the case of the antecedent null model (p = 0.032), suggesting a statistically significant negative influence of adult trees on saplings when the former are located closer than 8.5 m from the saplings in the study area (p = 0.045). In the case of the toroidal shift, the K rx (d) function is outside the 95% quantile bounds between 1.5 and 7.5 m.

Intertype Functions
The intertype Krx(d) function for the spatial structure of saplings and adult trees revealed a negative influence of the adult trees on the sapling density when the distance between the individuals of both cohorts is shorter than 11 m as the Krx(d) function is below zero (Figure 4). Beyond this distance, the function remains close to zero indicating an independent spatial pattern between the two cohorts. However, the Krx(d) function is only outside the 95% quantile bounds at distances between 1.5 and 8.5 m in the case of the antecedent null model (p = 0.032), suggesting a statistically significant negative influence of adult trees on saplings when the former are located closer than 8.5 m from the saplings in the study area (p = 0.045). In the case of the toroidal shift, the Krx(d) function is outside the 95% quantile bounds between 1.5 and 7.5 m.  The intertype mark correlationK rs mm (d) function for the dbh of the adult trees and the sapling density is within the 95% quantile bounds across the studied distances ( Figure 5, p = 0.536). This indicates that sapling density is independently spatially distributed with respect to the size of the adult trees. The intertype mark correlation ( ) rs mm K d function for the dbh of the adult trees and the sapling density is within the 95% quantile bounds across the studied distances ( Figure 5, p = 0.536). This indicates that sapling density is independently spatially distributed with respect to the size of the adult trees. Figure 5. Intertype mark correlationK rs mm (d) function (continuous line) between the diameter at breast height of the adult trees (stems with dbh ≥ 10 cm) and the sapling density (stems with height ≥ 1.30 m and dbh < 10 cm) and the 95% quantile bounds corresponding to the null model.
The intertype mark varianceL rs v (d) function is within the 95% quantile bounds of the toroidal null model without mark reassignation across the studied distances (red bounds in Figure 6, p = 0.304). This confirms the results of the intertype mark correlationK rs mm (d), the sapling density is not statistically related to the diameter of the adult trees. However, the intertype mark varianceL rs v (d) function lies below the 95% quantile bounds of the toroidal null model using the moving window algorithm with a 5 m window size to reassign the dbh of the trees at distances between 2 and 10 m (purple bounds in Figure 5, p = 0.001). All of this indicates that the sapling density is negatively related to the position of the adult trees at short distances, although this relationship is independant of the dbh of the adult trees.

Functional Predictor within a Generalized Additive Model Framework
The variable selection process revealed that the model with the lowest AIC was that containing the functional predictor N i=1 f 1 (D ji ) · dbh i /AreaIn30 j and the spatial smoother f 2 X j , Y j . The largest reduction in AIC is due to the inclusion of the spatial term in the model (Table 2). This suggests that the spatial distribution of the P. sylvestris saplings is primarily associated with the unmeasured variables such as microsite conditions (i.e., low stoniness) and the spatial pattern of previous stages of development and secondly to the size and position of adult trees under the studied conditions. The final model explained 47.4 % of the deviance and both the functional predictor and the spatial component had a significant effect on the number of saplings (p < 0.0001). In the case of the functional predictor, Figure 7 reveals a negative association between the diameter of adult trees and saplings when the distance between both cohorts is 0-6 m. Beyond this range, the coefficient turns positive, reaching the maximum at 11 m. This means that the largest number of saplings are expected to be found when the distance between adult trees and saplings is 11 m. The confidence intervals are, however, quite wide and contain the zero, indicating that the effect of adult trees on saplings is not so strong. The basis dimension of the functional predictor (10.0) is larger than the efficient degrees of freedom (4.8), ensuring that the underlying process, in this case, the sapling distribution, is represented reasonably well.
In bold, the selected model.
variables such as microsite conditions (i.e., low stoniness) and the spatial pattern of previous stages of development and secondly to the size and position of adult trees under the studied conditions. The final model explained 47.4 % of the deviance and both the functional predictor and the spatial component had a significant effect on the number of saplings (p < 0.0001). In the case of the functional predictor, Figure 7 reveals a negative association between the diameter of adult trees and saplings when the distance between both cohorts is 0-6 m. Beyond this range, the coefficient turns positive, reaching the maximum at 11 m. This means that the largest number of saplings are expected to be found when the distance between adult trees and saplings is 11 m. The confidence intervals are, however, quite wide and contain the zero, indicating that the effect of adult trees on saplings is not so strong. The basis dimension of the functional predictor (10.0) is larger than the efficient degrees of freedom (4.8), ensuring that the underlying process, in this case, the sapling distribution, is represented reasonably well. Table 2. Summary of the variable selection process according to Akaike's Information Criterion (AIC) and the percentage of deviance explained.
Variables included in the model during the selection process AIC Percentage of deviance explained (%) α 3327.3 0 In bold, the selected model.
The spatial component accounted for the remaining spatial trend of the saplings and, together with the functional predictor, removed the spatial autocorrelation of the residuals (data not shown).  The spatial component accounted for the remaining spatial trend of the saplings and, together with the functional predictor, removed the spatial autocorrelation of the residuals (data not shown).

Pinus sylvestris Sapling Distribution over the Fellings of the Uniform Shelterwood Method
In this study, we presented the first attempt at describing the spatial dependence of the number of saplings of P. sylvestris (sapling density) and adult trees during regeneration fellings under the uniform shelterwood method. In previous studies, this dependence has been studied in a monospecific pine forest which is regenerated through the group shelterwood method. Moreno-Fernández et al. [20] found stronger spatial relationships between the size of adult trees and saplings than that revealed by the GAM used in the present study; the number of saplings decreases significantly when the distance from adult trees is less than 7-8 m. In contrast, Montes et al. [18] detected weaker spatial dependence between both cohorts during shelterwood fellings. The intertype mark varianceL rs v (d) and the intertype mark correlationK rs mm (d) revealed that the dbh of the adult trees does not play a determinant role in the recruitment processes in this forest. This is not surprising as the forest has an even-aged structure with homogeneous distribution of adult tree dbh. The selection of trees that are maintained during the regeneration fellings may be based both on their spatial distribution to attain openings large enough and on the stem production characteristics (large dbh, straight stem without branches) of the individuals that are able to provide the seed for the new generation. However, the density of mother trees affects significantly the recruitment establishment in their proximity, pointing to the management of the tree density as the key factor be controlled during the regeneration process. The spatial repulsion between saplings and adult trees at shorter distances suggests that adult trees must be removed when the new individuals grow into saplings, but an excessive density reduction during the preparatory cuttings may favor the establishment of a dense herb layer that may hinder seedling establishment [45].
The clustered pattern of the new individuals of P. sylvestris has already been reported for seedlings and saplings of this species as the younger individuals perform better in microsites with more favorable conditions [11,46]. Similar findings have been reported for other pine species, e.g., Pinus uncinata Ramond ex DC [12]. Carrer et al. [47] employed bivariate point pattern analyses to describe the spatial dependence between Pinus cembra L., Picea abies (L.) H. Karst. and Larix decidua Mill. recruits and the adult trees at different altitudes in the Alps. Overall, they report similar findings to ours, that is, the repulsion between adult trees and recruits at short distances except in the case of P. cembra at higher altitudes where the relationship becomes positive. Batllori et al. [12] also found spatial repulsion between P. uncinata saplings and adult trees of this species in the Pyrenees. Malone et al. [48], however, found attraction between adult trees and seedlings of conifers at very short distances (<0.5 m) and repulsion at distances beyond 1 m in a dry conifer forest located in Colorado (USA). Ledo et al. [13] used the K rx (d) function to analyze the spatial dependence between adult trees and new cohorts in multispecies forests under tropical conditions. Under the conditions present in the aforementioned studies, the spatial association between classes varies depending on the species [13] and spatial scale [48].
It is worth to note that our results come from a single plot so that they must be interpreted with caution. However, the spatial association between adult trees and sapling density of P. sylvestris over the uniform shelterwood regeneration fellings are in accordance with other studies carried out in P. sylvestris stands in Central Spain [18,20].

Method Comparison
In this work, we used two types of techniques to describe the spatial dependence of sapling density and adult trees: three intertype correlation functions and a functional predictor within a GAM framework. The intertype K rx (d) function allows a spatial relationship between sapling density and the spatial pattern of the adult trees to be detected [18]. The antecedent conditions show a greater range of spatial repulsion between the adult trees and the saplings, but under the regular pattern of the adult trees, this null model performs in a very similar way to the toroidal shift null model. The intertype mark correlationK rs mm (d) evaluates the spatial correlation between the mark of the adult trees (tree dbh in our study) and the sapling density [13], and the toroidal shift null model seems the most appropriate in this case because this model preserves the intra-class pattern of the dbh of the adult trees and the sapling density in the quadrats, removing the spatial association between both processes through random shifting. The new point pattern function proposed, the intertype mark varianceL rs v (d) function, which combines the information reported by intertype K rx (d) function and the intertype mark correlationK rs mm (d) function in a single function, disentangles the effect of the size and the spatial pattern of adult trees on the regeneration density, in a similar way to the univariate L m (d) function proposed by Montes et al. [37]. The similarity between the range of spatial repulsion and mark correlation detected by theL rs v (d) and the ranges of spatial repulsion and mark correlation detected by K rx (d) and theK rs mm (d) respectively indicates that the toroidal shift preserving the dbh associated with each adult tree is appropriate to assess the mark correlation, whereas the toroidal shift null model combined with the moving window algorithm to reassign the dbh to the adult trees is appropriate in this case for analyzing the spatial association between adult tree closeness and sapling density.
The linear functional predictor within the GAM is an analogous approach to theK rs mm (d) function, incorporating the product between the dbh and the distance between the sapling quadrat and the adult tree. In fact, both approaches detected negative, although non significant, spatial dependence between sapling density and adult trees. It is likely that the functional predictor of the GAM, N n=1 f 1 (D ji ) · dbh i /AreaIn30 j , is confounded with the spatial smoother, f 2 X j , Y j , as both terms deal with the spatial structure of the sapling density [20]. The spatial smoother, however, accounts for the spatial correlation of the residuals and allows the effect of unmeasured variables to be considered [49]. The main advantages of the GAM over the point pattern methods used here are: (i) variables defined in several forms may be incorporated in the GAM, such as smooth terms, variables with linear effects, tensor products of several variables, with varying coefficients or as functional predictors; and (ii) the GAM can incorporate several temporal and spatio-temporal structures which are of great interest where there is data with repeated measurements [40,50,51].
Therefore, it seems reasonable to use the intertype functions for data without repeated measurements, as well as when there is no information on other predictor variables (e.g., stoniness, soil humidity or shrub cover) or where the effect of other variables on the spatial dependence between cohorts is not of interest. Further studies focusing on other types of data, e.g., data dealing with spatial dependence of species or variables characterizing niche suitability are necessary to further our knowledge with regard to the performance of the approaches used here. In a future step, growth dynamics of both saplings and remaining adult trees (as in Montoro et al. [52]) after the regeneration fellings may be incorporated through space-time modeling approaches to provide a deeper insight of the regeneration process and stand response under the uniform shelterwood method.

Conclusions
In our study, we found that regeneration fellings under the uniform shelterwood system resulted in an even-aged stand (i.e., homogeneous distribution of adult trees by dbh) and sapling density is not expected to change according to adult tree size. Taking into consideration the dbh of adult trees to explain the spatial distribution of saplings is of particular interest when the size of the adult trees is more heterogeneous than we had in this study. However, as our findings resulted from a single plot, replications may provide further information on the spatial association between adult trees and saplings in P. sylvestris stands in the studied area. In addition, similar studies may be carried out across the range of the species to modulate the shelterwood system for the different areas under the forecasted climate change.
We found the intertype mark varianceL rs v (d) function to be a useful approach to disentangle the spatial effect of the marks of the classes (sapling density and adult trees size) and the effect of the spatial pattern of the classes. However, to determine whether there are other measured variables (e.g., stoniness, shrub cover) driving the underlying process (regeneration in our study case) and where repeated measurements are available, the GAM may be highly effective. Finally, comparisons among approaches must be extended to other types of data, such as multispecies data.