Braak, C.J.F., ter

Verboom, J.

Name not available

On deducing time parameters from spatial patterns: the (auto)logistic model in landscape ecological processes

Wageningen Yield 

  
 
 1
 ON DEDUCING TIME PARAMETERS FROM SPATIAL PATTERNS:  
 THE (AUTO)LOGISTIC MODEL 
  IN LANDSCAPE ECOLOGICAL PROCESSES. 
 
 
 Cajo J.F. ter Braak1,2 and Jana Verboom2 
 DLO-Agricultural Mathematics Group and DLO-Institute for Forestry and Nature Research 
Box 100, NL-6700 AC WageningenE-mail:C.J.F.ter.Braak@GLW.AGRO.NL 
 
ter Braak, C. J. F. and Verboom, J.  (1995). On deducing time parameters from spatial patterns: 
The (auto)logistic model in landscape ecological processes.  In: Spatial and Temporal Modelling 
in Agricultural Research, Proceedings of the Fourth HARMA workshop, R. W. Payne, (ed.) 
Harpenden, UK, IARC-Rothamsted. 
 
Abstract 
 
In this paper, a popular metapopulation model is critically examined by putting the model in the 
context of Markov random fields and the statistical analysis of binary lattice systems. The claim 
that the model can be used to estimate time-process parameters from spatial-pattern data is 
examined on a real data set where process information was available. 
 
 
1 Introduction 
 
In man-dominated landscapes, natural habitats occur in small, spatially separated fragments. An 
animal or plant species that is restricted to such fragments may form a metapopulation, that is, a 
set of local populations which interact via individuals that move among them (Levins, 1970). For 
conservation programs, we need to understand the dynamics of such metapopulations, in 
particular the prospect for survival of the species in the landscape (Verboom et al., 1993). For 
landscape design and management, we need to assess the effects of mitigating measures (e.g. 
adding landscape elements) or measures that fragment the habitat even further, e.g. construction 
of road and railways (Verboom et al., 1993). How does the size and configuration of the 
fragments influence the metapopulation dynamics? 
 
In one class of metapopulation models, fragments (in this context often termed patches) are 
either vacant or occupied. Dynamical processes are then the process of extinction of the local 
population at a patch and the process of colonization of an vacant patch from occupied 
neighbours. Hanski (1994a) proposed a 'practical model of metapopulation dynamics' in which 
the parameters that govern the colonization and extinction processes are estimated from 
'snapshot' data on the occupancy of the species in a set of patches. In his method, process 
parameters are thus derived from spatial pattern data. Hanski (1994a) acknowledged that some 
parameters may sometimes be aliased so that extra information is needed for the claim to hold, 
but this aspect is not emphasized in Hanski (1994b). Under slightly different formulations of the 
model, the extra information is always needed, as in Hanski's second model with 'rescue effect'. 
 
In this paper, we critically examine Hanski's models by putting the models in the context of 
Markov random fields and the statistical analysis of binary lattice systems (Besag, 1977). Related 
  
 
 2
models are discussed by Durrett and Levin (1994). 
 
 
 
3 Hanski's practical model of metapopulation dynamics 
 
 
Hanski (1994a) modelled the occupancy of each particular patch i (i=1...n) by a Markov chain 
with two states (occupied, vacant). Let Ci be the probability that patch i, when vacant, is 
colonized the next year and Ei is the probability that patch i, when occupied, becomes vacant. 
The matrix of transition probabilities at time t is thus 
 
  1-Ei Ci 
 Ti = 
  Ei 1-Ci 
 
If Ci and, perhaps, Ei are functions of the occupancy of neighbouring patches, Ti varies in time. 
To derive a stationary probability that the patch is occupied (pi), we assume a constant 
neighbourhood so that Ti is a constant in time. Then the stationary equation is such that the 
probability of being occupied at time t+1 equals that probability at time t, i.e.  
 
 pi = Ci (1-pi) + (1-Ei) pi 
 
Solving for pi yields the stationary probability: 
 
 pi = Ci / (Ci + Ei) = 1 / (1 + Ei / Ci) (2.1) 
 
Hanski (1994a) noted the following problem with (2.1). If Ci and Ei both approach 1, then pi 
approaches 0.5. The reason for this is that the model does not allow for 'simultaneous' extinctions 
and colonizations. Hanski (1994) therefore introduced a 'rescue effect' by changing the realized 
extinction probability from Ei to (1-Ci) Ei. In terms of an explicit Markov chain, suppose that, in 
the annual cycle of events, a population may get extinct and later in the year, after reproduction 
in other patches, become colonized. So extinctions happen only in the first half of the year and 
colonizations only in the second half. Then we have for the first and second halves of the year 
the matrices T1i and T2i of transition probabilities, respectively, with 
 
  1-Ei 0   1 Ci 
 T1i =   T2i =  
  Ei 1   0 1-Ci 
 
and Ei and Ci are probabilities on semi-annual basis. The Markov chain which uses T1i and T2i 
alternatingly has on an annual basis the transition probability matrix that is the product of T2i and 
Ti1, i.e. 
 
   (1-Ei) + CiEi   Ci 
 T2iT1i =  
   Ei(1-Ci)  1-Ci 
  
 
 3
 
The stationary equation is thus 
 
 pi = [(1-Ei) + CiEi] pi  + Ci (1- pi)  
 
whence, as an alternative to (2.1), the model with rescue effect is 
 
 pi = Ci / (Ci + Ei - CiEi) = 1 / (1 + Ei (1-Ci)/Ci) (2.2) 
 
 In further specification of the models, Hanski (1994a) toke into account the spatial 
arrangement of patches, their size and occupancy status. The extinction probability Ei was 
modelled by  a function of the area Ai as follows 
 
  Ei = ε Ai-β if Ai > A0 = ε1/β (2.3a) 
 
 Ei = 1  if Ai ≤ A0 (2.3b) 
 
with parameters ε and β (0 ≤ ε ≤ 1 and β ≥ 0). Thus, the extinction probability decreases with 
the area. Extinction is certain when the patch is smaller than size A0 (ε1/β). The colonization 
probability Ci of patch i was assumed be a function of the number of migrants Si arriving to patch 
i per year: 
 
 Ci = 1/(1 + γ / Si2) (2.4) 
 
with γ > 0. Function (2.4) is an s-shaped logistic function in Si to allow for the Allee effect 
(Hanski, 1991). The number of migrants Si is unknown but is modelled in turn by a weighted 
sum over all occupied patches, namely 
 
 Si = Σj yj Aj exp(-αdij) (2.5) 
 
where yj is the 0/1 indicator for the occupancy state of patch j, Aj is the area of patch j and dij is 
the distance between patch i and patch j. In (2.5), the number of migrants from an occupied patch 
j is proportional to its area (Aj), and so to its potential population size, and inversely related to its 
distance to patch i. Equation (2.5) expresses the connectedness of patch i to the other occupied 
patches.  
 
The models so obtained are fitted by maximum likelihood to 'snapshot' spatial data on the 
occupancy, i.e. to the binary data {yi, i = 1...n}, indicating whether patch i is vacant (0) or 
occupied (1), while treating the data as independent. In Hanski (1994a) the parameter α is 
determined apriori from dispersion data.  For small {Ci}, Ci and the odds Ci/(1-Ci) are close, so 
that it is difficult to distinguish between the models (2.1) and (2.2) from practical data. In the 
model (2.2) with rescue effect, the parameters ε and γ are aliased, except when the cut-off in 
(2.3b) is in force; without (2.3b) only the product εγ is estimable. This aliasing may thus also be 
a practical problem in model (2.1). If this problem occurs, external information is needed. For a 
subset of patches, one may have occupancy data for two or more consecutive years. From such 
data one may calculate the number of turnover events T (extinctions and colonizations) and 
  
 
 4
resolve the alias by the extra estimation equation T -E(T) = 0 where E(T) is the expected number 
of turnovers in the model (Hanski, 1994a). Alternatively, one may have an idea about A0 in (2.3), 
the largest area below which the population in a patch that is just colonized will certainly go 
extinct next year (Hanski, pers. comm.). From A0 and β, one obtains ε, so that also γ becomes 
estimable. 
 
If the data are from a metapopulation at dynamic equilibrium between extinctions and 
colonizations and if all parameters are estimable, time-process parameters are obtained from 
spatial pattern data.  
 
 
5 Comments on Hanski's model 
 
 
The model is derived for a space-time process, but estimated from spatial data only. In the 
estimation, the spatial occupancy data {yi} play a role as predictors via (2.5) and as response 
variables as in auto-models (Besag, 1974). From the derivation of the model, the occupancy data 
used in (2.5) should have been from the previous year. The implicit assumption in the method is 
thus that the variation in the predictor values {Si} is small across subsequent years (Gyllenberg 
and Silvestrov, 1994). This assumption is trivially satisfied if Ci and 1-Ei are close to 0. The 
stationary probabilities (2.1) and (2.2) are derived under the assumption that Ei and Ci do not 
vary in time. However, Ci will vary in time. Again the implicit assumption is that the variation in 
{Si} is small. Gyllenberg and Silvestrov (1994) study the quasi-stationary distribution for the 
occupancy {yi}. 
 
It appears from the explicit formulae given for {pi} in Hanski (1994a,b) that equation (2.3b) is 
not used while fitting the models. Without the cut-off specified by (2.3b), the model (2.2) with 
rescue effect has the form of a logistic model, namely 
 
 logit (pi) = β0 + β1 log(Ai) + β2 log(Si) (3.1) 
 
with β0 = - log(εγ), β1 = β and β2 = 2. The method of fitting employed by Hanski (1994a) thus 
amounts to a logit regression with an offset. In a more extended version of the model β2 is a free 
parameter. 
 
In the method of fitting by logit regression, the autocorrelation among the spatial occupancy data 
is neglected. Consequently, the method of fitting does not maximizes a likelihood, but a 
pseudolikelihood (Besag, 1995, Preisler, 1993). Maximum pseudolikelihood estimators are in 
general still consistent but the standard errors must be adapted. The standard errors given in 
Hanski (1994a,b) thus cannot be trusted without further justification. 
 
Model (3.1) is not a proper auto-logistic model in the sense of Besag (1974). As a consequence 
of the Hammersley-Clifford theorem (Besag, 1974), the conditional probabilities given by model 
(3.1) do not define a joint probability distribution for the spatial data {yi}; the set of conditional 
probabilities is internally inconsistent. With pairwise interactions among patches, the only 
conditional model for binary data that yields a valid joint probability model is of the form 
  
 
 5
(Besag, 1974; Preisler, 1993) 
 
 logit pi = γi + Σj≠i δij yj (3.2) 
 
with δij = δji. An example of a proper auto-logistic model is 
 
 logit pi = γ0 + γ1 log(Ai) + γ2 Ai Si (3.3) 
 
Model (3.3) is of form (3.2) as can be seen by defining 
 
 δij = γ2 Ai Aj exp (-α dij). (3.4) 
 
The question is therefore whether logit regression can still yield sensible estimates of the 
parameters in (3.1).  
 
One may argue that the inconsistency of (3.1) in defining a joint probability distribution is not 
insurmountable because the model is meant to be valid for the spatial-temporal process, and not 
for the spatial data per se. Therefore, data were simulated by the spatial-temporal process and 
analyzed as spatial 'snapshot' data by Hanski's model (50 patches, 200 independent simulations). 
The configuration and size of the patches was as in Hanski (1994a: Fig. 7). Two sets of 
parameter values were used; one set as in Hanski (1994a, Fig. 7) and one set with much higher 
turnover rates. Hanski's estimators of the colonization and extinction turned out to be almost 
unbiased. For application of the model it is important to note that the variation in the estimates 
was huge.  
 
 
7 Application of Hanski's model to nuthatch data. 
 
 
Both models (3.1) and (3.3) were fitted to real data on a nuthatch metapopulation (cf. Verboom 
et al. 1991). Hanski's model (3.1) fitted better to the data than (3.3) as judged on the basis of the 
pseudo-deviance.  
 
Using α = 0.1 and β2 = 2, models (2.1) and (2.2) were fitted to data of each of six consecutive 
years in turn (1988-1993). We experienced numerical problems in all fits of model (2.1) [with 
2.ab, 2.4 and 2.5] by GENSTAT FITNONLINEAR, despite the good starting values obtained by 
logit regression. For two years no estimates could be obtained, indicated by stars in the table 
  
 
 6
 below. The fit of model (2.2)/(3.1) presented no problems. For model (2.2) an expert guess of A0 
of 0.5 ha was used to separate ε from γ. The expected number of turnovers between two 
consecutive years was estimated from each fit. The expected and observed number of turnovers 
(E(T) and T) were as follows: 
 
year  88/89 89/90 90/91 91/92 92/93 
 
E(T) 
(2.1) * 2.5 * 42.0 23.4 
(3.1) 12.5 11.3 5.8 20.1 21.6 
 
T 14 10 12 24 15 
 
By fitting (3.1) with β2 free, a deviance test on β2 = 2 was obtained. In none of the six cases, 
there was statistical evidence against β2 =2. Of course, the validity of the test is hampered by the 
fact that the fit is by pseudo-likelihood instead by regular maximum likelihood. The utility of 
Gibbs sampling of the spatial-temporal process with missing data (Augustin et al., 1994) is under 
investigation. 
 
References 
 
Augustin, N. H., Mugglestone, M. A. and Buckland, S. T. (1994). The role of the Gibbs sampler 
in modelling the spatial distribution of wildlife from incomplete survey data. 
Proceedings of the 9th International Workshop on Statistical Modelling, Exeter, 
8 pp.  
Besag, J. E. (1974). Spatial interaction and the statistical analysis of lattice systems (with 
discussion). J. R. Statist. Soc. B 36, 192-236. 
Durrett, R. and Levin, S. A. (1994). Stochastic spatial models: a user's guide to ecological 
applications. Phil. Trans. R. Soc. Lond. B 343, 329-350. 
Gyllenberg, M. and Silvestrov, D. S. (1994). Quasi-stationary distributions of a stochastic 
metapopulation model. J. Math. Biol. 33, 33-70. 
Hanski, I. (1991). Single-species metapopulation dynamics: concepts, models and observations. 
Biological Journal of the Linnean Society 42, 17-38. 
Hanski, I. (1994a). A practical model of metapopulation dynamics. Journal of Animal Ecology 
63, 151-162. 
Hanski, I. (1994b). Patch-occupancy dynamics in fragmented landscapes. TREE 9, 131-135. 
Levins, R. (1970). Extinction. Pages 75-108 in Some mathematical questions in biology, vol II, 
Rhode Island: American Mathematical Society. 
Preisler, H. K. (1993). Modelling spatial patterns of trees attacked by bark-beetles. Applied 
Statistics 42, 501-514. 
Verboom, J., Metz, J. A. J. and Meelis, E. (1993). Metapopulation models for impact assessment 
of fragmentation. In Landscape Ecology of a Stressed Environment, C. C. Vos, 
and P. Opdam (eds), 172-191. London: Chapman & Hall. 
Verboom, J., Schotman, A., Opdam, P. and Metz, A. J. (1991). European nuthatch 
metapopulations in a fragmented agricultural landscape. OIKOS 61, 149-156. 
 


On deducing time parameters from spatial patterns: the (auto)logistic model in landscape ecological processes

Abstract

Similar works

Full text

Available Versions

Name not available

Wageningen Yield