Forest Ecology and Management 108 Ž1998. 57–62
Modelling interactions between trees by means of field observations Aila Sarkka ¨ ¨ b
a,)
, Erkki Tomppo
b
a UniÕersity of Gothenburg, 411 80 Gothenburg, Sweden Finnish Forest Research Institute, Unioninkatu 40 A, 00170 Helsinki, Finland
Received 17 September 1997; accepted 19 November 1997
Abstract The Gibbs model and the pseudo-likelihood estimation method are applied to study interactions between trees of pine and spruce forests in the southern half of Finland. The coordinates of trees have been measured on small sample plots in the original data. Only some trees from the plots are needed in parameter estimation through field measurements. q 1998 Elsevier Science B.V. All rights reserved. Keywords: Spatial pattern of trees; Gibbs model; Pairwise interaction; Field observations; Pseudo-likelihood method
1. Introduction The relative spatial distribution of trees has an important role in many fields of forestry. Examples of purposes where it can be utilized are planning of forest inventory design, construction of growth models of trees or stands and problems related to forest regeneration and thinning methods. It is of importance when studying spreading of diseases or in producing of high-quality timber. In order to simulate artificial forests for different research purposes, knowledge of spatial pattern of trees is also needed. A suitable mathematical model for locations of trees is the spatial point process model. The simplest model is the homogeneous Poisson process, where complete spatial randomness is assumed. More complicated models and mathematical tools for parameter estimation are needed if the Poisson hypothesis is )
Corresponding author.
rejected. Interactions between trees are modelled here by using a homogeneous Gibbs point process. In this family of point processes interactions are taken into account and described by a potential function. We concentrate on pairwise interaction models, where a pair potential function describes interactions. In order to find out the level and strength of interaction, a pair potential function is parameterized and the parameters are estimated applying some estimation method. Non-parametric methods can also be applied. Most of the methods suggested in the literature are based on mapped data. In practical forestry or forest research, however, it may be difficult or too expensive to have data in this form. It is more convenient to collect data in a form of what is called field observations: measurements are carried out in arbitrarily chosen trees or randomly chosen sample points on the field. Tomppo Ž1986. suggested the application of the Takacs–Fiksel method and Sarkka ¨ ¨ Ž1995. the application of the pseudo-likeli-
0378-1127r98r$19.00 q 1998 Elsevier Science B.V. All rights reserved. PII S 0 3 7 8 - 1 1 2 7 Ž 9 7 . 0 0 3 4 5 - 9
A. Sarkka, ¨ ¨ E. Tomppor Forest Ecology and Management 108 (1998) 57–62
58
hood method in estimating a pair potential function for data collected through field observations. The model and the estimation methods are briefly recalled in Section 2. Interactions between the trees in pine and spruce forests in the southern half of Finland are modelled in this paper. In Section 3, we present the original INKA sample plot data measured by the Finnish Forest Research Institute and results of the interaction study are in Section 4. The coordinates of the trees on the sample plots are known instead of the entire maps of forest stands. Only one tree from each chosen sample plot is considered here in order to demonstrate the parameter estimation from field measurements. The parameters are estimated through the pseudo-likelihood method.
Fig. 1. A typical pair potential function of an inhibitive Žsolid. and of an attractive Ždashed. process and a typical mixture of these two Ždotted..
2. Preliminaries 2.1. Gibbs pairwise interaction model The Gibbs point process model is a natural tool to describe interactions between trees. For each point configuration w , the distribution f Ž w . measures how much more likely the configuration w is for that process than for a Poisson process ŽStoyan et al., 1987, pp. 156–159.. This comparison of configurations of interaction processes with those of the Poisson forest is natural because the trees of a Poisson process do not interact. We concentrate on pairwise interaction processes F in R 2 , where a pair potential function f :w0,`. ™ R is applied to describe interactions between the trees. A so-called local energy is related to the probability to add a point x to the configuration w and is given by: EŽ x ,w . s a q
Ý
f Ž 5 x y y 5. ,
yg w
where a is a constant chemical activity connected to the intensity of the process and 5 x y y 5 is the distance between x and y. We consider stationary and isotropic processes, where the pair potential function depends only on the distance between the points. The sign and shape of the pair potential function indicate a possible inhibition or attraction between
trees and its strength Žsee Fig. 1.. Roughly speaking, positive values indicate inhibition and negative ones attraction. There is no interaction at distance r with f Ž r . s 0. A well-known fact is that the pairwise Gibbs process is not usually the best model for highly clustered patterns: e.g., under hard-core condition parameters can always be estimated but testing is difficult Žsee e.g., Gates and Westcott, 1986; Møller, 1994.. In the case of a bounded observation window W ; R 2 , the density function of the pairwise interaction process is: f Ž wW . s
1 Z
exp Ž yU Ž w W . . ,
where w W s x 1 , x 2 , . . . , x n 4 consists of the locations of n points in W, Z is a scaling factor, and U Ž w W . s n a q Ý f Ž 5 x i y x j 5. i-j
the energy needed for the configuration w W ŽRipley, 1988.. In the present study, a special case of the pairwise interaction process, the multiscale pairwise interaction process, is applied ŽPenttinen, 1984.. Let 0 F R 0 - R 1 - . . . - R m , m G 1, be fixed real numbers. R 0 is the fixed hard-core radius Žminimum possible intertree distance. and R m s R the fixed interaction
A. Sarkka, ¨ ¨ E. Tomppor Forest Ecology and Management 108 (1998) 57–62
radius Ždistance after which the trees do not interact.. The pair potential function is given as:
°` if 0 F r F R - r F R , i s 1, . . . ,m ¢0 if r ) R
f Ž r . s~ a if R i
0
iy1
i
Ž 1.
m
where a1 , . . . ,a m are the parameters to be estimated. In order to avoid confusion with points of the process and points of R 2 , we refer throughout this paper to trees of the process or of the pattern and to points of R 2 . 2.2. Estimation through field obserÕations The pair potential function is parameterized and the parameters are estimated by means of some estimation method. Statistical methods developed for mapped data may be of limited use in practical forestry. Measurement of local information through sampling is often more practicable or even the only possibility. The classical field observation methods can be divided into two groups: areal counts and distance measurements. The sample units are usually sparsely located in order to get approximative independence needed in statistical reasoning. Only the Tacaks–Fiksel ŽTF. method ŽTomppo, 1986. and the pseudo-likelihood ŽPL. one ŽSarkka, ¨ ¨ 1995. have been shown to be applicable in pair potential estimation through nonmapped field observations. The TF estimation method is based on the equality of two expectations, i.e.,
lE 0 u Ž F . s E u Ž F . exp Ž yE Ž 0,F . . ,
Ž 2.
where l is the intensity of the process, u a non-negative measurable test function and E the local energy ŽFiksel, 1984.. E 0 is an expectation with respect to a conditional distribution given a tree at a given point and E an expectation with respect to the distribution of the process. For example, if uŽ w . is the number of trees in a certain distance, then E 0 w uŽF .x is the expected number of trees calculated from an arbitrary tree, the tree itself excluded, and Ew uŽF .x the expected number calculated from an arbitrary point on the study area. For the Poisson process these two expectations are equal but for general Gibbs processes weighting is needed.
59
The idea is to estimate both sides of Eq. Ž2. and use the method of least squares in order to minimize the squared difference. The expectations can be estimated by corresponding sample means calculated from measurements connected to chosen trees and points. Mapped data are not needed if the pair potential and test functions are chosen in a suitable way ŽTomppo, 1986.. Measurements that are needed depend on the pair potential and the test functions ŽSarkka, ¨ ¨ 1995.. The PL method was first developed for Markov random fields ŽBesag, 1974. and later for special type of Gibbs point processes ŽBesag, 1978.. The idea is to construct a PL function and maximize it with respect to the parameters. The PL function for Gibbs point processes can be derived from the one for the cell processes. The study area is divided into a grid of n c square cells of area D such that there is at most one tree in each cell. The PL function for the system of random variables n s Ž n1 ,n 2 , . . . ,n n c . is: PL Ž n;u . s P P Ž n i , < n j , j / i . , i
where n i is the number of trees in the ith cell. Therefore, the PL function is a product of conditional probabilities to have one or no trees in one cell given the number of trees in other cells. Under rather general assumptions, e.g., a hard-core distance exists, as D tends to zero the process converges to a point process ŽBesag et al., 1982.. In the case of Gibbs processes, maximization of the logarithm of the PL function leads to the maximization of log PL Ž u ; w W . s y
Ý
EŽ x ,w _ x .
xg w W
y
H exp Ž yE Ž j , w . . d j W
with respect to the parameter u . The PL estimation equations coincide with Eq. Ž2. if the test function u is chosen as a derivative of the pair potential function with respect to the parameters ŽSarkka, ¨ ¨ 1989; Grabarnik and Sarkka, ¨ ¨ 1992; Diggle et al., 1993.. Therefore, the PL estimation equations for unmapped data are discrete versions of Eq. Ž2. and can be derived on the basis of it ŽSarkka, ¨ ¨ 1995.. The chosen pair potential function determines the measurements that are needed.
60
A. Sarkka, ¨ ¨ E. Tomppor Forest Ecology and Management 108 (1998) 57–62
3. Data: INKA growth sample plots The empirical data consist of the permanent INKA sample plots established by the Finnish Forest Research Institute for growth study purposes. The plots are distributed all over Finland. Only plots located in the southern half of the country are considered in this study Žsee Fig. 2.. Forests can be divided into homogeneous stands according to tree and stand characteristics, such as tree species composition, age and size of trees, site fertility, etc. The size of stand vary between 0.2 and 5 ha, typically between 1 and 3 ha. Some type of stands have first been chosen for the study and then a cluster of three permanent sample plots are placed on each study stand. The sizes of the plots depend on the number of trees per hectare. Several characteristics have been measured describing both the forest stands and the trees on the sample plots. Examples of stand variables are dominant tree species, mean characteristics of trees, soil type and site fertility of soil. From each tree of the plots, coordinates, species, diameter and many other variables of importance in practical forestry, are measured. The INKA growth sample plots are a special case of field observation: small regions of the study area are investigated carefully but the entire map of the whole area is unknown.
4. Results The pairwise interaction process approach and INKA sample plots are applied in the study of interactions among the trees in pine and spruce forests in the southern half of Finland. A multiscale step function model Ž1. with nine steps at even intervals is applied. The minimum inter-tree distance is used as a hard-core radius R 0 . The length of interaction radius R is fixed as 3.5 m ŽTomppo, 1986.. The model is appropriate both in inhibitive and attractive cases because a1 , . . . ,a m can have both positive and negative values. The parameters are estimated through the PL method. The estimates of the pair potentials are given but goodness of fit of the models is not tested. It can be done by using Ripley’s L function ŽRipley, 1981.. The plots in pine and spruce dominated stands have been divided into three groups according to the mean diameter of the trees on the plot: Ž1. 0–10 cm, Ž2. 10–20 cm and Ž3. more than 20 cm. The mean diameter of trees and the age of stands are highly correlated and therefore this grouping allows to study the spatial pattern of trees in forests of different ages. In each group, the pattern formed by the dominant and subdominant trees are studied, the suppressed ones are not taken into account. Furthermore, some arbitrarily chosen trees and sample points are needed. Tomppo Ž1986. suggested that the minimum adequate sample size is from 60 to 100. Sample points and trees are collected as follows: from each sample plot of interest one arbitrary point and one arbitrary tree are chosen, both of them far enough from the boundary in order to avoid edge effects. Notice that all trees on the sample plots are known here and therefore the choice of arbitrary trees is straightforward.
4.1. Pine forests
Fig. 2. Locations of INKA stands in the southern half of Finland.
The estimated pair potential curves and the corresponding simulated forests show the spatial pattern of pine forests of different age. Young pine forests Žclass 1. have been regenerated by planting and are therefore regular ŽFig. 3a and b.. The trees inhibit each other up to about one meter. After that, interac-
A. Sarkka, ¨ ¨ E. Tomppor Forest Ecology and Management 108 (1998) 57–62
61
tion becomes negligible. Naturally originated trees will appear when stand gets older Žclass 2. causing attraction in small distances as can be seen from the Fig. 3c and d. Some of the trees have died or have been removed in old pine forests Žclass 3.. The pattern of trees is regular again ŽFig. 3e and f.. Pine needs much light which fact also favors a regular pattern. 4.2. Spruce forests Trees in young spruce forests Žclasses 1 and 2. attract each other in small distances ŽFig. 4a and c.. Trees form small clusters as can be seen in simulated
Fig. 4. Pair potential estimate of the classes Ža. 1, Žc. 2 and Že. 3 and a typical realization of the spruce forest of classes Žb. 1, Žd. 2 and Žf. 3 corresponding to the pair potential function in Ža., Žc. and Že., respectively.
patterns in Fig. 4b and d. Spruce does not require as much light as pine and is therefore capable to grow in small clusters. In old forests Žclass 3. trees form a regular pattern, probably caused by thinning ŽFig. 4e and f..
Acknowledgements Fig. 3. Pair potential estimate of the classes Ža. 1, Žc. 2 and Že. 3 and a typical realization of the pine forest of classes Žb. 1, Žd. 2 and Žf. 3 corresponding to the pair potential function in Ža., Žc. and Že., respectively.
We would like to thank Dr. Michel Goulard from INRA, Avignon, for his useful comments.
62
A. Sarkka, ¨ ¨ E. Tomppor Forest Ecology and Management 108 (1998) 57–62
References Besag, J., 1974. Spatial interaction and the statistical analysis of lattice systems Žwith discussion.. J. R. Stat. Soc. B 36, 192– 236. Besag, J., 1978. Some methods of statistical analysis for spatial data. Bull. Int. Stat. Inst. 47, 77–92. Besag, J., Milne, R., Zachary, S., 1982. Point process limits of lattice processes. J. Appl. Prob. 19, 210–216. Diggle, P.J., Fiksel, T., Grabarnik, P., Ogata, Y., Stoyan, D., Tanemura, M., 1993. On parameter estimation for pairwise interaction point processes. Int. Stat. Rev. 62, 99–117. Fiksel, T., 1984. Estimation of parameterized pair potentials of marked and non-marked Gibbsian point processes. Electron. Informationsverarb. U. Kybernet 20, 270–278. Gates, D.J., Westcott, M., 1986. Clustering estimates for spatial point distributions with unstable potentials. Ann. Inst. Stat. Math. A 38, 123–135. Grabarnik, P., Sarkka, ¨ ¨ A. Ž1992.. On parameter estimation of marked Gibbs point processes. Preprints from the Department of Statistics, University of Jyvaskyla, ¨ ¨ 5.
Møller, J. Ž1994.. Markov chain Monte Carlo and spatial point processes. Department of Theoretical Statistics, Institute of Mathematics, University of Aarhus, 293. Penttinen, A. Ž1984.. Modelling interaction in spatial point patterns: parameter estimation by the maximum likelihood method. Jyvaskyla ¨ ¨ Studies in Computer Science, Economics and Statistics, 7. Ripley, B.D. Ž1981.. Spatial statistics. Wiley, Chichester. Ripley, B.D. Ž1988.. Statistical inference for spatial processes. Cambridge Univ. Press, Cambridge. Sarkka, ¨ ¨ A. Ž1989.. On parameter estimation of Gibbs point processes through the pseudo-likelihood method. Reports from the Department of Statistics, University of Jyvaskyla, ¨ ¨ 4. Sarkka, ¨ ¨ A., 1995. Pseudo-likelihood approach for Gibbs point processes in connection with field observations. Statistics 26, 89–97. Stoyan, D., Kendall, W.S., Mecke, J. Ž1987.. Stochastic geometry and its applications. Wiley, Chichester. Tomppo, E. Ž1986.. Models and Methods for analysing spatial patterns of trees. Communicationes Instituti Forestalis Fenniae, 138.